Implementación paralela de un método numérico híbrido global en ecuaciones de convección difusión

Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería. Vol. 11,4, 695-705(1995) * Departamento de Matemática Aplicada Universidad Politécnica de Valencia Camino de Vera s/n 4 6071 Valencia, Spain ** Departamento de Sistemas Informáticos y Computación Universidad Politécnica de Valencia Camino de Vera s/n 46071 Valencia, Spain RESUMEN Presentamos un método numérico híbrido global, y su implementación en paralelo sobre el multiprocesador de memoria compartida Alliant FX/80, para el tratamiento de ecuaciones de convección-difusión en regimen transitorio en dos dimensiones. El método se basa en la discretización realizada en el caso estacionario con bases de Legendre adaptadas al operador, conjuntamente con la utilización de la transformada de Laplace y un algoritmo por bloques de reducción cíclica o par-impar que implementamos en paralelo para la resolución de los sistemas tridiagonales por bloques asociados al problema. Esto nos permite obtener soluciones aproximadas de forma rápida y sin restricciones sobre el paso de tiempo, al contrario de lo que sucede en otros métodos como diferencias finitas o elementos finitos. Se comparan las prestaciones obtenidas en el algoritmo por bloques de reducción cíclica con las de otros algoritmos por bloques basados en eliminación Gaussiana y técnicas del tipo divide y vencerás. SUMMARY We present a global hybrid method and its parallel implementation on the shared memory multiprocesor Alliant FX/80 for the treatment of two dimensional transient convection-diffusion equations. The method is based on the discretization of the steady state using Legendre basis adapted to the operation, together with the use of the Laplace transformation and a cyclic block reduction or odd-even algorithm implemented in parallel for the solution of the tridiagonal block system associated to the problem. This allows us to obtain approximate solutions very Recibido: Enero 1995 OUniversitat Politecnica de Catalunya (España) ISSN 0213-1315 696 E CHECA MARTÍNEZ Y V. HERNÁNDEZ GARCÍA fast and without restrictions on the time step comparing favourably with other methods such as finite differences and finite elements. The obtained results using the cyclic block reduction are compared to algorithms based on Gaussian elimination and divide and win techniques. INTRODUCCI~N La elección de bases para la discretización de una ecuación diferencial, como por ejemplo la de difusión o convección-difusión, es uno de los aspectos claves en la aproximación numérica de la solución. En particular, en problemas de conveccióndifusión aparecen peculiaridades y dificultades propias debidas al término de convección que destruye la simetría del operador y que complica su resolución. La elección propuesta en este trabajo conduce a estructuras matriciales sencillas, que evitan los problemas computacionales asociados a la obtención de sistemas con matriz coeficiente densa, mal condicionada etc. En la sección siguiente consideramos bases polinomiales que son combinación de polinomios de Legendre. Se consiguen así estructuras matriciales en la discretización que nos permiten obtener resultados precisos en régimen estacionario, tanto en el caso unidimensional como bidimensiona13. En la sección 2 presentamos un algoritmo numérico híbrido global para el caso transitorio, en el contexto de las ecuaciones citadas, que se apoya en el tratamiento realizado en la sección 1 y que está relacionado con los resultados expuestos en4Y5. Para llevar a cabo este nuevo enfoque hacemos uso formal de la transformada de Laplace y de un algoritmo de inversión desarrollado en8 cuya implementación en paralelo ha sido estudiada en3. Finalmente, en la sección 3 desarrollamos y analizamos, en el contexto de las ecuaciones citadas con la discretización realizada, un algoritmo de reducción cíclica por bloques desarrollado e implementado sobre el sistema multiprocesador de memoria compartida Alliant FX/80. Las prestaciones de este algoritmo se comparan con las de un algoritmo del tipo divide y vencerás por bloques propuesto también en3. 1. RÉGIMEN ESTACIONARIO. DISCRETIZACION EN BASES NO CONVENCIONALES Consideramos el subespacio vectorial de polinomios de grado menor o igual que N + 2 siguiente VN = {P(x) E PN+2(x) : P(-1) = P(1) = O) y la base {Q3(x)), Q3(x) = L3+2(x) - L3(x) , siendo Lk(x) los polinomios de Legendre de grado k. Estas bases están incluidas en la literatura científica dentro de otras más generales denominadas hieráticasl. Es bien conocido que los polinomios Q3 (x) se pueden expresar como (1 - X~)R~(X), y por tanto se cumple que se anulan en la frontera. Igualmente se verifica que Q; (x) = (23 + 3) L3+i (x), propiedad importante en los desarrollos que realizamos para los problemas tratados. Otras fórmulas explícitas deducidas de las anteriores y que son útiles en nuestros desarrollos aparecen en3. El esquema está basado en aproximaciones polinomiales en subespacios cuyas bases seleccionadas son producto tensorial de bases del caso unidimensiona13. De este modo se IMPLEMENTACIÓN DE UN MÉT. NUM. HÍBRIDO GLOBAL EN EC. DE CONV.- DIF. 697 obtienen estructuras matriciales, susceptibles de ser resueltas mediante métodos propios de la computación en paralelo. Consideramos el problema1' ( ~(5,742) = (71/2)2 - x2 donde aparece un término de difusión dado por las derivadas de segundo orden, junto con uno de convección, dado por un fluido con velocidad -2a, paralelo al eje y. Transformando este problema en el intervalo (-1,l) x (-1,l) se tendrá, en un planteamiento Galerkin, después de aplicar un teorema de Green, y considerando que la solución se desarrolla en bases que son producto tensorial de las utilizadas en el caso unidimensional, las ecuaciones siguientes con I=(-l,l)x(-l,l), k=O ,..., N yj=O ,..., M El sistema que se obtiene tiene matriz coeficiente tridiagonal por bloques, con bloques en la diagonal principal con estructura pentadiagonal y el resto con estructura tridiagonal. Tenemos que añadir, que como sucederá en el problema transformado que se considera en la sección siguiente, si en el problema anterior se considera en la ecuación correspondiente un término adicional de la forma su(x, y), como será el caso del problema transformado de conveccióndifusión, se obtiene que todos los bloques tienen estructura pentadiagonal. Cuando se hace uso de la transformada de Laplace es importante disponer de un método de inversión que sea eficiente, pues en caso contrario esta herramienta analítica pierde, en el contexto de los métodos numéricos, parte de su eficiencia. Utilizamos un algoritmo de inversión8 que se está mostrando eficiente, como puede verse en4Y5 donde 698 E CHECA MARTÍNEZ Y v HERNÁNDEZ GARCÍA se combina con el método de diferencias finitas y con el de elementos finitos, todo ello en el contexto de la ecuación del calor. Caso unidimensional: Método numérico híbrido global Consideramos la ecuación de convección-difusión en régimen transitorio au(x,y) D++V~=O a2u(x t) ~(0, t) = 1, u(1, t) = O, u(x, O) = O {TConvirtiendo la ecuación en otra equivalente con condiciones de Dirichlet N homogéneas y transformando el problema, admitiendo que h(x, S) = an(s)Qn(x) es una aproximación a la solución del problema transformado se tiene, en un planteamiento Galerkin estándar, las ecuaciones siguientes N 4~ an(s)(S1 (2j + 3)(2n + 3)~n+i(~)~j+i(x)d~- n=O -1 1 -S_l (2~ + 3) (Ln+2(~) - Ln(x))dx+ + S S_:(~n+2(x) - L~(x))(L~+~(x) - ~j (x))dx = ~c, con ko = -2((v/s) - 1/2), kl = -11 y Ic:, = O si j > 2 Estas ecuaciones conducen a un sistema lineal con estructura pentadiagonal. Los resultados obtenidos para distintos tiempos aparecen descritos en la siguiente gráfica. Los resultados han sido contrastados con la solución analítica desarrollada en forma de serie. ~{l~(,l'lll~h THANs PHO(_'EsC) I)I\IJb,I(CiION 1 ' 08 o 11 ( ' u (5 11 t C a. 04 ( 1 n o2 u -1 -0 8 -0 6 -0 4 -0 2 O U '2 O 11 O O O X 1 0 t = 10, + t = 100, o t = 2000 La situación se complica cuando el término fuente es una distribución como la S de Dirac, pues añade problemas en torno al punto donde está centrada dicha distribución. También en estos casos, en régimen transitorio, se consiguen buenos resultados, aunque con aproximaciones de orden alto7. Extensión al caso bidimensional. Soluciones numéricas por secciones y evolución temporal en puntos fijos En este caso, tomamos como referencia un problema de difusión que analiza la dispersión de calor en una placa delgada con temperatura inicial dada y temperatura en la frontera constante. Posteriormente se obtiene la estructura matricial en el caso de convección-difusión que se estudiará en la sección siguiente en un contexto de computación en paralelo. La ecuación en derivadas parciales que modeliza el fenómeno anterior, con valores particulares en la frontera y en la condición inicial, es con k = 519. Aplicando la transformada de Laplace a este problema, definimos sobre el problema transformado el funcional con f = 30 y A6 = sil-k (g + $1 siendo D t H~(O), y O el abierto (-1,l) x(-1,l). Aplicando un teorema de Green sobre la expresión del funcional y tomando como base el producto tensorial de bases polinomiales en el caso unidimensional, Qjk(x, y) = Qj(x)Qk(y) donde Qj(x) = (1 - x2)~j(x) y Qk(y) = (1 - y2)~k(y), buscamos una a proximación global de la forma Esta aproximación puede implementarse en general, aunque aparezcan pequeños problemas por las dimensiones N,M con que se desarrolla. Construimos la aproximación para el par (N, O) suponiendo por las condiciones e interpretación física del caso tratado que la solución es suave en la dirección y, tomando en esta solamente el polinomio Qo(y). La aplicación del teorema de Green conduce a la siguiente ecuación Después de sustituir derivar parcialmente respecto de ajo y simplificar se obtiene 700 E. CHECA MARTÍNEZ Y V. HERNÁNDEZ GARCÍA un sistema tridiagonal desacoplado siendo las incógnitas impares nulas a(2k+l)~, k = 0,1,2, .... Los resultados obtenidos, tomando dos funciones polinómicas en la dirección y, se describen seguidamente atendiendo a ciertas posiciones nodales3. Tpo t = 1.2 h Posición N21 N28 N41 N6 1 N80 Si el desarrollo se realiza con Qo(y),Qi(y), se obtendrán dos sistemas desacoplados análogos al anterior, donde el único que ofrece información como consecuencia de la discretización es el que antes hemos descrito. En efecto, basta observar que las incógnitas adicionales son cero ya que J'~ Ql(y)Qo(y)dy = O y la estructura de las derivadas parciales respecto a ajl son análogas a las tratadas. Si introducimos los tres primeros elementos de la base en la dirección y, se puede dividir el sistema resultante en subsistemas desacoplados, algunos de los cuales se pueden suprimir por tener la solución trivial. El único que no conduce a esta solución es el dado por las ecuaciones 8Fii aFü correspondientes a aa;M, aaJzM La estructura que aparece conduce a una matriz por bloques de la forma (C E) siendo A,B,C,D tridiagonales. El estudio realizado conduce a los resultados ofrecidos en la tabla anterior. Bajo el encabezamiento S. Aprox. aparecen los resultados más precisos, y el número de grados de libertad son 15. En realidad solo contribuyen 10 funciones de la base, cinco en la dirección x y dos en la dirección y, que son Qo(y),Q2(y). El desarrollo realizado sirve como base para establecer el siguiente resultado en problemas de convección-difusión. PROPIEDAD 1 Dado el problema + Lu = f (x, y) con condiciones de Dirichlet homogéneas y condición inicial nula, siendo L el operador de convección-difusión Lu = Au - a$, se cumple que la discretización en las bases de Legendre adaptadas en el caso bidimensional con un enfoque Galerlcin y con transformación previa, conduce a un sistema Ax = b con estructura tridiagonal por bloques, siendo, en general, las matrices que constituyen estos bloques, matrices banda. La justificación del resultado anterior se obtiene teniendo en cuenta la discretización realizada con estas bases en el caso bidimensional, para el operador Lu = Au - a-" ay Y la del caso tra.nsitorio del de difusión, en cuyo problema transformado aparece el factor su que se corresponde con E al aplicar la Transformada de Laplace. Brush (A t = 0.05) 0.185 1.039 1.286 1.938 0.922 TL-EFH 0.166 0.998 1.116 1.697 0.813 TL-DF 0.175 1.176 1.198 1.831 0.871 S. Exacta 0.173 1.065 1.186 1.812 0.862 S. Aprox. mod(4,2) 0.172 1.064 1.186 1.810 0.862 IMPLEMENTACI~N DE UN MÉT. NUM. HIBRIDO GLOBAL EN EC. DE CONV.- DIF. 701 Destacamos el hecho de que el planteamiento realizado permite el estudio del fenómeno considerado en secciones del dominio para un tiempo cualquiera, sin que ello suponga el conocer la aproximación a la solución del problema en puntos contiguos. Además conviene añadir el interés que supone fijar un punto cualquiera del dominio y estudiar aquí la evolución temporal del fenómeno, con el mismo coste computacional que para un único valor de la variable tiempo, pues los procesos son independientes. En las secciones anteriores se ha realizado un estudio de cierto tipo de ecuaciones de convección-difusión que aparecen en la práctica y que, con el enfoque realizado, conducen en el caso de coeficientes constantes a sistemas de ecuaciones lineales con matriz coeficiente tridiagonal por bloques, susceptibles de ser resueltos en paralelo de forma rápida y eficiente. Los resultados del apartado siguiente no sólo tienen aplicación en los problemas señalados sino en contextos más generales, como al aplicar esquemas de diferencias finitas o elementos finitos en ecuaciones de tipo elípticog. En ambos se obtienen sistemas de ecuaciones lineales con matriz coeficiente tridiagonal por bloques con bloques de estructura determinada. Algoritmo de reducción cíclica o par-impar por bloques Consideremos el sistema de ecuaciones lineales Ax = b, siendo A una matriz tridiagonal por bloques expresada de forma genérica como A = (O,. . . , B,, A,, C,, . . .). El desarrollo es válido igualmente para la ecuación matricial AX = B, siendo B una matriz. Suponemos que se verifican condiciones de regularidad sobre las submatrices (bloques) de forma tal que puedan llevarse a efecto las operaciones siguientes. Consideramos el estudio del caso general, donde A es una matriz tridiagonal por bloques, y donde particionamos x y b de acuerdo a los bloques de la matriz A. Suponemos que el número de bloques que tiene A en la diagonal es N = 2m - 1, sin que esto suponga ninguna restricción, pues en la implementación se ha tomado un valor de N cualquiera, completando el sistema con filas de la forma (. . . , O, 1, O, . . .) y término independiente (1,1,. . . , I)~. Una ecuación característica del sistema tendrá la forma B,X,-~ +Aix, + C,X,+~ = b, donde i = 1,2, ..., N siendo Bl = CN = O y x0 = XN+~ = O Consideramos pues la fila i-ísima de A, con i par,(. . . ,O, Bi, Ai, Ci, O, . . .). La idea básica consiste en premultiplicar, respectivamente las filas i - 1, i + 1 por matrices convenien tes tales que al restarlas ambas a la fila i se obtenga otra con la siguiente estructura (. . . , O, B:, O, A:, O, c:, O, . . .). Estas operaciones crean un subsistema tridiagonal a partir de las 2m-1 - 1 ecuaciones de índice par. Aunque las ecuaciones impares no se consideren, las incógnitas impares se pueden obtener posteriormente a partir de las pares por sustitución regresiva. Dadas ahora las 2m-1 - 1 ecuaciones con incógnitas pares solamente, se puede eliminar un nuevo grupo de incógnitas, consiguiendo así, un conjunto de 2m-2 - 1 ecuaciones que implican incógnitas cuyo índice es un múltiplo de 4. Este esquema puede repetirse hasta obtener una única ecuación para x2,-i y resolviendo el sistema correspondiente se habrá calculado xp-l. Llegados a este punto se calculan las incógnitas en orden inverso al que fueron eliminadas, mediante sustitución regresiva. Supuesto que estamos en la iteración Ic del proceso, con vistas a obtener las fórmulas generales para la posterior implementación del algoritmo se deben calcular matrices (H~)!, (H~)!, tales que El algoritmo de resolución implementado puede consultarse en3. Implementación y comparación de resultados sobre el multiprocesador de memoria compartida alliant FX/80 Consideramos en primer lugar un ejemplo test de matriz tridiagonal escalar que figura en2 y posteriormente los casos que surgen al tratar los procesos de conveccióndifusión discretizados en bases hieráticas. Estos procesos conducen a estructuras matriciales por bloques, con los bloques de la diagonal principal pentadiagonales y los que están en la subdiagonal y superdiagonal principal con estructura tridiagonal. Aquí no se tendrá en cuenta este hecho para darle a nuestro tratamiento un carácter más general. Se ha comprobado que el algoritmo de reducción cíclica con 8 procesadores ofrece, frente al algoritmo de descomposición LU con 1 procesador, un speed-up de hasta 3.5 en el caso N1 = 2016, siendo N1 el número de bloques en la diagonal principal, y cada bloque de dimensión tres en este caso. Seguidamente presentamos el speed-up conseguido en la implementación realizada del algoritmo de reducción cíclica AL(~OKI'ih10 PAKIMPAK 8 I I I I I I I BLOQUES EN LA DIAGONAL Se observa el alto speed-up de este algoritmo, aunque cabe añadir que como algoritmo secuencia1 es francamente ineficiente, al contrario de lo que sucederá con IMPLEMENTACIÓN DE UN MÉT. NUM. HÍBRIDO GLOBAL EN EC. DE CONV.- DIF. 703 el algoritmo DAC que con un procesador supera, en algunos casos, a la descomposición LU en tiempos de ejecución3. La siguiente gráfica muestra una comparación en tiempos de los dos algoritmos 0 Parimpar Global, + Dac Combinado Del estudio realizado conclumos que, en el caso general de matrices tridiagonales por bloques, el mejor algoritmo de los estudiados es el de reducción cíclica. La situación puede cambiar si el sistema es tridiagonal por bloques con bloques de estructura determinada, como sucede en el tratamiento realizado en los procesos de conveccióndifusión donde la estructura de las submatrices es tridiagonal y pentadiagonal. Tenemos que señalar también en nuestro contexto que la mejor precisión se obtiene con eli algoritmo de reducción cíclica, igualada por el algoritmo divide y vencerás, DAC, para tamaños determinados, y lo que es más importante con N1 dado, la elección que se haga de los parámetros IB e IK en la descomposición en submatrices pues Nl=IB.IK, siendo IB el número de subdivisiones en la matriz e IK el número de matrices de cada subdivisión7. Presentamos a continuación un estudio del speed-up para ambos algoritmos, en función del número de bloques en la diagonal principal y del tamaño de cada bloque. Algoritmo: REDUCCIÓN CÍCLICA. Speed up N1 = Número de bloques diagonal principal N = Tamaño de bloque