Sobre la paralelización de problemas elípticos

Bol. Soc. Esp. Mat. Apl. no31(2005), 31–52 Sobre la paralelizaci´on de problemas el´ıpticos M. C. Calzada1, I. I. Albarreal2, J. L. Cruz1, E. Fern´ andez-Cara2, J. R. Galo1, M. Mar´ ın1 1Dpto. de Inform´atica y An´alisis Num´erico, Universidad de C´ordoba, Campus de Rabanales, Ed. C2-3, E-14071 C´ordoba 2Dpto. de Ecuaciones Diferenciales y An´alisis Num´erico, Universidad de Sevilla, Aptdo. 1106, 41080 Sevilla [email protected], [email protected], [email protected], [email protected], [email protected], [email protected] Resumen En este trabajo se recogen dos de los problemas que aparecen a la hora de resolver num´ericamente las ecuaciones deNavier-Stokes. La discretizaci´on en la variabletemporal en dichas ecuaciones nos conduce, en cada etapa de tiempo, a un problema de Burgers y un problema de Stokes generalizado, en las variables espaciales. La aplicaci´on de un m´etodo de punto fijo y un m´etodo de tipo gradiente conjugado, respectivamente, nos lleva a la resoluci´on de un elevado n´umero de problemas de tipo Helmholtz. As´ı, en una primera parte, recordamos un algoritmo paralelo que resuelve num´ericamente la ecuaci´on de Helmholtz. Por otro lado, la utilizaci´on de diferencias finitas centradas sobre un mismo mallado regular para la velocidad y la presi´on en la resoluci´on del problema de Stokes puede generar soluciones oscilantes, en particular presiones esp´ureas. En una segunda parte, por tanto, se presenta un esquema de discretizaci´on para el problema de Stokes, demostr´andose mediante un an´alisis asint´otico del error que ´este regulariza la soluci´on1. Fecha de recepci´on: 31/01/05 1Una versi´on previa de este trabajo aparecer´a en las Actas del I Workshop sobre “Recientes Avanaces en el An´alisis y Control de Ecuaciones Diferenciales No Lineales”, celebrado en C´ordoba en febrero de 2004. 31 32 M. C. Calzada y otros 1 Introducci´on En la resoluci´on de muchos problemas de gran complejidad, aparecen de forma frecuente ecuaciones de tipo Helmholtz: αu −β∆u=fen Ω ⊂IRN,(1) siendo α≥0, β > 0, funa funci´on conocida y Ω un dominio acotado bi o tridimensional. La ecuaci´on se completa con condiciones de contorno, que pueden ser de tipo Dirichlet, u=gsobre Γ = ∂Ω,(2) de tipo Neumann, Fourier o una mezcla de ambas. Aunque ´este es uno de los problemas m´as conocidos y estudiados, es claro que una peque˜na ganancia en el tiempo de resoluci´on de (1) puede en la pr´actica suponer una gran mejora en el tratamiento de problemas m´as complejos. Pensemos en las ecuaciones de Navier-Stokes, para las cuales hemos desarrollado varios algoritmos donde este hecho se pone de manifiesto (v´ease por ejempo [3, 4]). El m´etodo que se presenta sigue la idea de los denominados m´etodos ADI (Alternating Directions Implicit methods). Estos m´etodos fueron formulados inicialmente por Peaceman y Rachford [15] en el caso bidimensional y por Douglas y Rachford [7] en el caso tridimensional y, posteriormente, fueron generalizados al caso N-dimensional por Douglas y Gunn [6]. Para describir la idea de estos m´etodos, consideremos por ejemplo el problema de Helmholtz (1)–(2). Consideremos fijados una descomposici´on del operador L = αI−β∆ en la forma L = L1+ L2con Ln=α 2I−β∂2 ∂x2 n (n= 1,2) (3) y un esquema en diferencias finitas centradas con paso de discretizaci´on h en las dos direcciones espaciales. Un m´etodo de paso fraccionado (de tipo de Peaceman-Rachford) para esta descomposici´on es un m´etodo iterativo que comienza con U0arbitrario y, para m≥0, calcula Um+1 en dos pasos. As´ı, para Umdado, primero se obtiene Um+1 2como soluci´on de (I+τA1)Um+1 2=τf + (I−τA2)Um(4) y en segundo lugar se resuelve (I+τA2)Um+1 =τf + (I−τA1)Um+1 2.(5) Aqu´ı, τ > 0,yAnes la matriz correspondiente a la discretizaci´on del operador Ln2. As´ı, en el primer paso se resuelve un problema discreto correspondiente al operador α 2I−β∂2 ∂x2 1 , 2Por simplicidad, hemos denotado de nuevo fla funci´on discretizada. Paralelizaci´on y problemas ell´ıpticos 33 es decir, un problema en la direcci´on x1. An´alogamente, en el segundo paso se resuelve un problema en la direcci´on x2. Esto es lo que motiva el nombre de m´etodo de direcciones alternadas. Con una adecuada ordenaci´on de los nodos del mallado, se consigue que las matrices I+τAn(n= 1,2) sean diagonales por bloques, siendo cada bloque una matriz tridiagonal. Esta idea se puede aplicar, en general, a la hora de resolver un sistema algebraico Au =f(6) resultante (por ejemplo) de la discretizaci´on de un problema el´ıptico, siempre que la matriz Ase pueda descomponer en la forma3 A= p X n=1 An,(7) siendo las matrices An(1 ≤n≤p) m´as f´aciles de invertir que la propia A. Esta generalizaci´on da lugar a los llamados m´etodos de descomposici´on o m´etodos secuenciales de paso fraccionado (v´ease [13]). En este caso, para m≥0, conocido Um,el c´alculo de Um+1 se realiza resolviendo los siguientes sistemas: Um+1 p−Um τm +A1³Um+1 p−Um´=−αÃp X n=1 AnUm−f!, (8) Um+n p−Um+n−1 p τm +An³Um+n p−Um´= 0, n = 2,...,p siendo α > 0 y {τm}una sucesi´on de par´ametros reales. Como acabamos de ver, los m´etodos de paso fraccionado (8) son algoritmos iterativos secuenciales, en los que cada iteraci´on se subdivide en ppasos que se realizan sucesivamente seg´un el esquema: Um→Um+1 p→Um+2 p→ · · · → Um+p−1 p→Um+1.(9) Es decir, hasta que no se calcula Um+i pno se puede calcular Um+i+1 p. Por ello, la paralelizaci´on de estos algoritmos s´olo puede llevarse a cabo al m´as bajo nivel, a la hora de resolver los sistemas (ver por ejemplo [18] y [11]). Nosotros abordamos la resoluci´on de los problemas de Helmholtz generalizando y ampliando a dominios arbitrarios un m´etodo introducido por Lu y otros [12] formulado para problemas con condiciones de tipo Dirichlet y para dominios que resultan de la uni´on de rect´angulos. Nuestro objetivo global es desarrollar un esquema num´erico que permita la paralelizaci´on al nivel de las variables espaciales. Este m´etodo, que denominamos m´etodo impl´ıcito de direcciones simult´aneas oSDI (Simultaneous Directions Implicit method), descompone 3En general, pno tiene porqu´e coincidir con la dimensi´on Nen la que se trabaje. 34 M. C. Calzada y otros el problema original, independientemente de la dimensi´on espacial del dominio en el que est´e formulado, en un conjunto de problemas unidimensionales. Se sustituye as´ı un problema en derivadas parciales por una familia de problemas diferenciales ordinarios todos ellos con la misma estructura y por tanto con las mismas dificultades num´ericas, independientemente de la dimensi´on espacial en la que se trabaje. Adem´as, los problemas resultantes son independientes entre s´ı y pueden resolverse en paralelo. El m´etodo SDI que hemos desarrollado es aplicable a dominios Ω arbitrarios. En [17] se aplica en el caso bidimensional, efectu´andose previamente una “rectificaci´on”del abierto Ω, de tal manera que se aproxima el dominio original por una uni´on de rect´angulos. En nuestro caso no se efect´ua esta rectificaci´on y el mallado se construye a partir de una ret´ıcula ortogonal de partici´on arbitraria, admitiendo todo tipo de condiciones de contorno. Para ello se ha desarrollado un mallador (ver [9]), tanto en el caso 2D como en el 3D, que permite disponer de una estructura de datos que contribuye a una mayor eficiencia del m´etodo. El m´etodo es econ´omico en sus necesidades de almacenamiento, por lo que es compatible con problemas de gran talla. Adem´as, en [9] se demuestra su car´acter “smoother”, que explica que se adapte muy bien al uso de t´ecnicas multigrid para acelerar la convergencia sea cual sea la geometr´ıa del dominio considerado. Aqu´ı solo expondremos el algoritmo paralelo propuesto cuando la ecuaci´on de Helmholtz se completa con condiciones de contorno de tipo Dirichlet. La consideraci´on de otro tipo de condiciones de contorno rompe en general, con fronteras no necesariamente paralelas a los ejes coordenados, la independencia por direcciones al ser condiciones netamente interespaciales. En [9] se encuentra desarrollada una t´ecnica v´alida para condiciones de otro tipo. ´ Esta consiste en un m´etodo iterativo que inicialmente necesita la resoluci´on de una sucesi´on de problemas con condiciones de Dirichlet y posteriormente se reduce a la aplicaci´on del m´etodo SDI a un ´unico problema con condiciones de Dirichlet dependientes de la iteraci´on. Recientemente se ha demostrado la convergencia de este m´etodo. Hay que destacar que la resoluci´on de los problemas de Helmholtz con condiciones de Neumann tiene inter´es no s´olo por s´ı misma, sino tambi´en cuando se buscan precondicionadores que permiten acelerar la convergencia del gradiente conjugado al aplicar ´este, por ejemplo, a la resoluci´on de un problema de Stokes. Como es de preveer, disponer de un m´etodo paralelo eficiente para la resoluci´on de los problemas de Helmholtz conduce a una notable mejora en la resoluci´on de problemas m´as complejos, como puede ser el caso de los problemas de Stokes o Navier-Stokes. En nuestro caso, para resolver el problema de Stokes, reformulamos ´este como un problema de m´ınimos al que aplicamos un m´etoddo iterativo de tipo gradiente conjugado propuesto en [10] y [5]. Este m´etodo conduce a la resoluci´on de un gran n´umero de problemas de Helmholtz y/o Poisson. Por lo tanto, disponemos de una manera de resolver un problema de Stokes N-dimensional exclusivamente a partir de problemas unidimensionales independientes (problemas diferenciales ordinarios), con igual esquema de resoluci´on y dificultad num´erica independientemente de la dimensi´on del problema original. Se observa entonces que la paralelizaci´on del m´etodo es posible a tres niveles, respectivamente correspondientes a la Paralelizaci´on y problemas ell´ıpticos 35 descomposici´on de los problemas vectoriales por componentes y a los dos niveles de paralelizaci´on impl´ıcitos al m´etodo SDI. En [5], se resuelven los problemas de Helmholtz y Poisson discretizando por elementos finitos, us´andose un mallado el doble de fino para las velocidades que para la presi´on, al objeto de verificar la condici´on inf-sup. En nuestro caso, discretizamos los problemas unidimensionales resultantes utilizando diferencias finitas y el mismo mallado para todas las inc´ognitas. Esto hace que no se verifique la condici´on infsup, lo que provoca la aparici´on de presiones esp´ureas, que disminuyen la precisi´on de la aproximaci´on de la presi´on calculada sin que la aproximaci´on de la velocidad tenga que verse afectada. Las presiones esp´ureas se manifiestan en la representaci´on gr´afica de las isobaras como l´ıneas de contorno con oscilaciones. En la pr´actica, seg´un la discretizaci´on que se efect´ue, el tipo de “test” considerado y el car´acter de su soluci´on, puede que expl´ıcitamente no se observe la presencia de las mismas. En la secci´on 3.2 efectuaremos un an´alisis de esta dificultad introduciendo el desarrollo asint´otico del error cometido y siguiendo la metodolog´ıa propuesta por Wetton [19] en su estudio del m´etodo de proyecci´on. Wetton explica la causa de la aparici´on de presiones esp´ureas en base a la existencia de “errores alternantes” unidireccionales. Aqu´ı realizaremos un an´alisis similar, contemplando la existencia de estos errores pero con un car´acter multidireccional. Identificaremos los modos esp´ureos y, de esta manera, podremos buscar el filtrado adecuado para las presiones que haga que las oscilaciones se reduzcan. Este filtrado (ver [9]) se realiza siguiendo cada una de las direcciones espaciales, manteniendo el esquema conceptual establecido de reducir todo el tratamiento a problemas unidimensionales. 2 Un esquema paralelo para la ecuaci´on de Helmholtz Como acabamos de ver, los m´etodos de paso fraccionado (8) son algoritmos iterativos secuenciales, en los que cada iteraci´on se subdivide en psubetapas que se realizan sucesivamente. En [12], Lu et al. proponen un esquema de descomposici´on para la resoluci´on de (6) en el que los pasos fraccionados son independientes entre s´ı. A partir de una descomposici´on de la matriz Adada por (7) y de un par´ametro τ > 0, el paso de UmaUm+1 se realiza calculando Um+1,n (1 ≤n≤p) a partir de (I+τAn)Um+1,n = I−τ p X k=1,k6=n Ak Um+τf, n = 1,...,p, (10) y haciendo, posteriormente, Um+1 =ω p p X n=1 Um+1,n + (1 −ω)Um,(11) donde ωes un par´ametro dado. Se trata de un m´etodo de descomposici´on en el cual los ppasos fraccionados pueden resolverse simult´aneamente y la 36 M. C. Calzada y otros paralelizaci´on se realiza a un nivel m´as alto, siendo esto mucho m´as interesante y efectivo. A este m´etodo lo denominaremos m´etodo paralelo de paso fraccionado oPFS (Parallel Fractional Step method). Esquem´aticamente puede representarse como sigue: Um Um+1,1 % · · · & & · · · % Um+1,p Um+1 (12) Cuando se aplica este m´etodo al problema de Helmholtz-Dirichlet N dimensional ½Lu=αu −β∆u=fen Ω ⊂IRN, u=gsobre ∂Ω(13) considerando la siguiente descomposici´on del operador L = N X n=1 Ln= N X n=1 µα NI−∂2 ∂x2 n¶, se tiene que, en cada paso (10) se resuelven Nproblemas independientes correspondiente cada uno de ellos al operador α NI−β∂2 ∂x2 n , al igual que suced´ıa con el m´etodo ADI pero, en este caso, los c´alculos se hacen en todas las direcciones simult´aneamente. De aqu´ı que a este m´etodo lo denominemos m´etodo impl´ıcito de direcciones simult´aneas oSDI (Simultaneous Directions Implicit method)4. Como ya indicamos en el m´etodo ADI, la descomposici´on elegida es conveniente, ya que las matrices resultantes I+τAnson m´as f´aciles de invertir que la matriz Adel problema original. En este caso, la obtenci´on de cada Um+1,n puede efectuarse resolviendo en paralelo ssistemas independientes con sdel orden de h−1, lo que se corresponde con la resoluci´on simult´anea en cada una de las l´ıneas del mallado en la direcci´on xn. Se tienen por tanto dos niveles de paralelismo que conducen a N·sposibles procesos simult´aneos y, dado que ses un n´umero elevado al ser en la pr´actica h peque˜no, obtenemos un algoritmo con un alto ´ındice de paralelizaci´on. Adem´as, los problemas a resolver para la obtenci´on de las Um+1,n tienen todos la misma estructura y presentan id´enticas dificultades num´ericas. Los m´etodos paralelos de paso fraccionado en general y los SDI en particular, seg´un la descripci´on (10)–(11) efectuada, dependen de dos par´ametros (que pueden variar en cada iteraci´on en el caso no estacionario); el par´ametro τdenominado par´ametro de evoluci´on y el par´ametro ω, conocido como par´ametro de coordinaci´on. En [9] se encuentra un estudio exhaustivo de los valores de estos par´ametros, 4En general, llamaremos m´etodos SDI a los m´etodos paralelos de paso fraccionado aplicados a problemas en los que en cada Lnintervengan derivadas en una ´unica direcci´on. Paralelizaci´on y problemas ell´ıpticos 37 que hacen que el m´etodo sea convergente y adem´as su velocidad de convergencia sea lo mayor posible. En [1] aparece un estudio del error para un caso particular del algoritmo totalmente discretizado. 2.1 Resoluci´on efectiva En esta secci´on efectuamos un an´alisis de los elementos necesarios para la resoluci´on efectiva de los problemas de tipo Helmhotz-Dirichlet. Mallado Como es natural, para la resoluci´on del problema (13) efectuaremos una discretizaci´on espacial. A partir de una ret´ıcula cartesiana R(∆), obtenemos un mallado de Dirichlet M(∆) = M1(∆) ∪... ∪MN(∆), de nodos N(∆) = NI(∆) ∪NF(∆),donde NI(∆) contiene los nodos interiores, NF(∆) contiene los nodos frontera y Mn(∆) = Mn [ j=1   en(n,j) [ l=1 Ij,l n ,1≤n≤N, siendo Ij,l n= ((al),(bl)) un segmento de IRNen la direcci´on espacial enisomorfo al intervalo de IR (al, bl).En las Figuras 1 y 2, aparece un ejemplo de M(∆) para un dominio bidimensional y en la Figura 3 un dominio tridimensional an´alogo. Figura 1: Un mallado M(∆) = M1(∆)∪M2(∆),para un dominio bidimensional. 38 M. C. Calzada y otros Figura 2: Dos segmentos integrables en la direcci´on e1. Figura 3: Algunos segmentos integrables en un dominio 3D. Detalle del algoritmo ∗Etapa de resoluci´on La etapa de resoluci´on del algoritmo (10), aplicada al problema (13), se corresponde con Nsubproblemas del tipo ((I + τLn)um+1,n =³I−τPN k=1,k6=nLk´um+τf en Ω ⊂IRN, um+1,n =gsobre Γ = ∂Ω,(14) donde Ln=α NI−∂2 ∂x2 n ,1≤n≤Nyumes conocida. Para el mallado M(∆), como Lnes unidimensional, para cada uno de los problemas (14) se tiene el correspondiente problema discreto ((I + τδn)Um+1,n =³I−τPN k=1,k6=nδk´Um+τf en Mn(∆), Um+1,n =gen NF(∆),(15) donde δnes el operador discreto correspondiente a Ln.Dado que las condiciones de contorno son de tipo Dirichlet, (15) es un conjunto de problemas Paralelizaci´on y problemas ell´ıpticos 39 independientes5 ((I+τδn)Um+1,n =³I−τPN k=1,k6=nδk´Um+τf en Ij,l n= ((al),(bl)) , Um+1,n ((al)) = g((al)) , Um+1,n((bl)) = g((bl)), (16) definidos en cada uno de los segmentos integrables Ij,l nque constituyen Mn(∆). A su vez, cada problema (16) es un problema discreto correspondiente a un problema de contorno ordinario (unidimensional) en la variable xndel tipo (³1 + τα N´v−τβv00 =˜ fen (al, bl), v(al) = va, v(bl) = vb,(17) donde v00 es la segunda derivada en la variable xny˜ fes un segundo miembro gen´erico. As´ı pues, Um+1,n en (16) queda determinada mediante la resoluci´on de, al menos, Mnproblemas de contorno ordinarios del tipo (17), cada uno de ellos en el correspondiente intervalo de IR,isomorfo al segmento integrable Ij,l n∈Mn(∆). Haciendo n= 1, ..., N, vemos que (16) es una familia de, al menos, M= N X n=1 Mnproblemas del tipo (17). Todos ellos poseen la misma estructura y, por tanto, presentan igual dificultad num´erica. Adem´as, son independientes entre s´ı por lo que pueden te´oricamente resolverse de forma simult´anea. Hagamos notar que Mes un n´umero relativamente grande, ya que depende del n´umero de puntos de las particiones consideradas en cada direcci´on espacial. As´ı pues, se tiene un alto nivel de paralelismo posible. ∗Etapa de coordinaci´on Determinadas las Um+1,n, con n= 1, ..., N, en los nodos N(∆) del mallado, la etapa de coordinaci´on del algoritmo (11) es tambi´en un proceso paralelizable ya que el c´alculo Um+1(P) = ω N N X n=1 Um+1,n(P) + (1 −ω)Um(P) (18) se puede realizar de forma simult´anea para cada P∈N(∆).Esta etapa de coordinaci´on es realmente la que nos obliga a utilizar una ret´ıcula cartesiana como base para la construcci´on del mallado. ∗“Test” de parada El “test” de parada considerado en el proceso iterativo es: kUm−Um+1k∞= m´ax x∈N(∆) |Um(x)−Um+1(x)|< ε, donde εes una tolerancia dada. 5La matriz del problema discreto es una matriz diagonal por bloques. 46 M. C. Calzada y otros donde las bu(s)=u(sp)−u(si) 2, u(s)=u(sp)+u(si) 2 son regulares e independientes de h. Decimos que este tipo de desarrollo es alternante. Esto da una explicaci´on precisa de qu´e es lo que ocurre cuando se pierde la regularidad discreta. As´ı pues, aparece una soluci´on global que oscila al cambiar de submallado ya que el desarrollo asint´otico del error contiene t´erminos alternantes que, en caso de que predominen frente a los t´erminos regulares, causan la p´erdida de precisi´on, v´ease [19]. Notemos que el desacoplamiento es debido a la discretizaci´on elegida del operador y el orden en el que aparece el t´ermino alternante es el orden en el que difieren las aproximaciones frontera consideradas en cada submallado. Los tipos de desarrollo asint´otico del error indicados pueden extenderse a cualquier problema de contorno bi o tridimensional. Como conclusi´on, se tiene que mediante este an´alisis podemos averiguar si un esquema de discretizaci´on es regularizante y, en el caso de tener un desarrollo alternante, se puede determinar, mediante una adecuada ponderaci´on de los valores de nuestra inc´ognita en un nodo y sus adyacentes, un patr´on de filtrado que permita aminorar la p´erdida de precisi´on indicada. En [19] se efect´ua un an´alisis del esquema de discretizaci´on de Chorin [2] aplicado, en este caso, al problema de Stokes. Se demuestra que el error es de segundo orden para la velocidad (con errores de tipo alternante de tercer orden) y s´olo de primer orden en la presi´on debido a la presencia de errores alternantes, que predominan sobre el t´ermino regular de segundo orden. Sin embargo, este esquema es aplicable s´olo a dominios rectangulares o compuestos por uni´on de rect´angulos, dado que es preciso calcular una divergencia discreta en los nodos frontera. 3.2 El esquema regularizante propuesto Para no extender la exposici´on, haremos el desarrollo en el caso bidimensional. De forma an´aloga puede hacerse para el caso N-dimensional. Para la resoluci´on de (19), definido en un dominio arbitrario, con objeto de obtener un esquema regularizante, consideraremos el problema discreto (20), donde:  Ωhes el dominio discreto determinado a partir de una ret´ıcula cartesiana regular, que ser´a el mismo para la velocidad y la presi´on.  ∆hes el operador de Laplace discreto bidimensional, definido por la aproximaci´on de cinco puntos de orden h2, es decir: ∆hzi,j =zi−1,j +zi+1,j −4zi,j +zi,j−1+zi,j+1 h2.(25)  Gh= (Gx h, Gy h) y Dhson los operadores gradiente y divergencia discreta, obtenidos mediante aproximaciones de segundo orden basadas Paralelizaci´on y problemas ell´ıpticos 47 en diferencias centradas: Gx hPi,j =Pi+1,j −Pi−1,j 2h, (26) Gy hPi,j =Pi,j+1 −Pi,j−1 2h, (27) Dh·U=Dh·(U, V )i,j =Ui+1,j −Ui−1,j 2h+Vi,j+1 −Vi,j−1 2h.(28) Calcularemos la presi´on exclusivamente en los nodos interiores, por lo que s´olo es necesario evaluar la divergencia y el gradiente discreto en dichos nodos. De aqu´ı que este esquema puede aplicarse en problemas con dominios no necesariamente rectangulares. Adem´as, se tiene lo siguiente:  El operador divergencia discreta Dhdesacopla el mallado 2D en cuatro submallados (ocho en 3D), pero ´este s´olo se aplica en los nodos interiores. Dado que los valores de la velocidad son conocidos sobre Γh, todos los submallados “ven” igual aproximaci´on en la frontera.  El operador gradiente discreto Gh(26)–(27) respeta los submallados anteriores, pero su aplicaci´on en los nodos adyacentes a la frontera requiere la elecci´on de una adecuada aproximaci´on. Por ejemplo, en los nodos de referencia (1, j), seg´un la notaci´on habitual, la aplicaci´on de (26) conducir´ıa a Gx hP1,j =P2,j −P0,j 2h,(29) es decir, ser´ıa necesario conocer la presi´on P0,j en los nodos frontera de referencia (0, j). Pero en estos nodos hemos exclu´ıdo el c´alculo de la presi´on y, por consiguiente, debemos realizar una aproximaci´on que efectuaremos mediante una extrapolaci´on unidireccional, dada por P0,j = 3P1,j −3P2,j +P3,j .(30) El error de truncamiento asociado es e0,j =p(0, jh)−P0,j =h3pxxx +O(h4).(31) As´ı pues, de (29) y (30) se tiene que Gx hP1,j =−3P1,j + 4P2,j −P3,j 2h,(32) que no es m´as que la f´ormula de diferencias progresivas de segundo orden. F´ormulas an´alogas a ´esta (diferencias progresivas o regresivas) se obtienen en los nodos adyacentes al resto de la frontera. Seg´un (31), (32) y expresiones an´alogas, cada uno de los submallados antes citados “ve” diferentes condiciones de contorno para la presi´on, que se 48 M. C. Calzada y otros diferencian en t´erminos de orden h3. Por consiguiente, el desarrollo asint´otico de la presi´on, aproximada por este esquema, tiene un t´ermino alternante de orden h3, es decir, un t´ermino de la forma bp(3) i,j =p(3)(ih, jh) + (−1)ip(3i)(ih, jh) + (−1)jp(3j)(ih, jh) + (−1)i+jp(3ij)(ih, jh) donde las p(3),p(3i),p(3j)yp(3ij)son funciones regulares e independientes de h. Se tiene entonces el siguiente resultado, cuya demostraci´on puede encontrarse en [9]: Teorema 1 Si la soluci´on (u, p)del problema de Stokes (19) es suficientemente regular, entonces la soluci´on (U, P )del problema discreto (20), correspondiente a los operadores (25)–(28) actuando en los nodos interiores y tomando sobre la frontera los valores extrapolados para la presi´on, tiene un desarrollo asint´otico del error dado por U=u+h2u(2) +h5b u(5) +··· ,(33) P=p+h2p(2) +h3bp(3) +··· ,(34) donde u(2) yp(2) son t´erminos regulares y b u(5) ybp(3) son t´erminos alternantes. Por tanto, las aproximaciones de Uy de Pson ambas de segundo orden en h. 3.3 Experiencias num´ericas “Test” en dominio rectangular Con objeto de comprobar el comportamiento num´erico del esquema propuesto y compararlo con el esquema de discretizaci´on de Chorin (aplicable s´olo en dominios rectangulares), consideraremos inicialmente el “test” con soluci´on anal´ıtica propuesto por J.T. Oden en [14]. Consiste en la resoluci´on de (19) en el cuadrado Ω = (0,1) ×(0,1), con α= 0, ν = 1 y felegida para que la soluci´on exacta venga dada por: u1(x, y) = 8x2(1 −x)2y(y−1)(2y−1), u2(x, y) = 8y2(1 −y)2x(x−1)(1 −2x), p(x, y) = 4(x−x2). (35) La presi´on toma valores en el intervalo [0,1], siendo las isobaras segmentos verticales. La implementaci´on pr´actica se efect´ua aplicando un algoritmo tipo gradiente conjugado en el que se integra el m´etodo SDI (ver [9]), de aqu´ı que el m´etodo sea iterativo. En la Figura 9 resumimos los resultados obtenidos al aplicar los esquemas anteriormente citados: el m´etodo de Chorin sin filtrar las presiones, el m´etodo de Chorin filtrando las presiones seg´un el patr´on obtenido en [9] y el esquema regularizante propuesto. En concreto, en cada columna de la gr´afica, representamos:  El comportamiento, respecto al n´umero de iteraciones, de la norma infinito del error exacto tanto del campo de velocidades como de la presi´on, para diferentes valores del paso de discretizaci´on h. Paralelizaci´on y problemas ell´ıpticos 49  Las isobaras para el caso h= 0,05.  El error exacto de la presi´on mediante un gr´afico de superficie cuando h= 0,05. Para el primer esquema (primera fila de gr´aficos en la Figura 9), podemos observar que la existencia de un t´ermino alternante de primer orden en la presi´on conduce a la p´erdida de precisi´on, que se refleja en la situaci´on final alcanzada en el error exacto y que es del orden te´orico demostrado en [19]. La velocidad no se ve afectada. Son muy significativas las oscilaciones de las isobaras, as´ı como el predominio del error de tipo alternante frente al regular, que se aprecia en el gr´afico de superficie. Tambi´en se puede observar que el error est´a asociado a los diferentes submallados. En el segundo esquema, el filtro aplicado a la presi´on logra que el t´ermino alternante de primer orden se traslade al tercer orden, por lo que pasa a predominar el t´ermino regular de segundo orden y, por tanto, se regulariza la presi´on (este hecho est´a demostrado en [9]). Figura 9: Resultados del “test” de Oden para diferentes esquemas. Podemos observar en la segunda fila de esa misma figura que el comportamiento de la velocidad es el mismo, pero la presi´on mejora su precisi´on, alcanz´andose el orden indicado y manifest´andose peque˜nos cambios en la convergencia debido al t´ermino alternante, ya de tercer orden. Las isobaras reflejan leves oscilaciones, debido m´as al error propio de un m´etodo aproximado que a la presencia de presiones esp´ureas. En el gr´afico de superficie, se observa el predominio del t´ermino regular de segundo orden y la presencia del error alternante en un orden superior. No se aprecia ya distinci´on entre los diferentes submallados. En el tercer caso, de acuerdo con lo demostrado anteriormente, se logra alcanzar el error regular de segundo orden, tanto para la velocidad como para la presi´on, no apareciendo oscilaciones en las isobaras. El gr´afico de superficie tambi´en refleja el comportamiento regular del error y no se observa 50 M. C.˜Calzada y otros dependencia de los submallados. Hemos realizado un n´umero de iteraciones muy superior al necesario, para poder intuir el comportamiento asint´otico del algoritmo. Omitimos la representaci´on gr´afica del campo de velocidades dado que, seg´un lo descrito anteriormente, ´este no se ve afectado en la pr´actica. “Tests” 2D y 3D en dominios no rectangulares Como ejemplo de “tests” bidimensional y tridimensional definidos en dominios no rectangulares, consideraremos la determinaci´on de un fluido est´atico (“No Flow”) sobre el que act´ua una fuerza gravitacional. En la Figura 10, reflejamos las presiones obtenidas mediante l´ıneas y superficies isob´aricas, respectivamente. Se puede comprobar, de nuevo, el car´acter regularizante del esquema propuesto, no apareciendo presiones esp´ureas. En [9], pueden encontrarse resultados con otros “tests” efectuados. Figura 10: Isobaras y superficies isob´aricas en los “tests” “No Flow” 2D y 3D. Agradecimientos Trabajo parcialmente financiado por la D.G.E.S., Proyecto BFM2003–06446– C02–02 y por la Junta de Andaluc´ıa, ACC–1204–FQM124. Referencias [1] I.I. Albarreal -Paralelizaci´on en tiempo y espacio de la resoluci´on num´erica de algunas ecuaciones en derivadas parciales, Tesis, Universidad de Sevilla (2004). [2] A.J. Chorin -On the convergence of discrete approximations to the Navier-Stokes equations, Math. Comp, 22 (1968), 745–762. Paralelizaci´on y problemas ell´ıpticos 51 [3] J.L. Cruz, M.C. Calzada, E. Fern´ andez-Cara, M. Mar´ ın -A parallel algorithm for solving the incompressible Navier-Stokes problem, Comp. Math. with Appl., 25 (9) (1993), 51–28. [4] J.L. Cruz, M.C. Calzada, E. Fern´ andez-Cara, M. Mar´ ın -A convergence result for a parallel algorithm for solving the Navier-Stokes equations, Comp. Math. with Appl., 35 (4) (1998), 71–88. 1995. [5] E.J. Dean, R. Glowinski -On Some Finite Element Methods for the Numerical Simulation of Incompressible Viscous Flow, Incompressible computational fluid dynamics (M.D. Gunzburger and R.A. Nicolaides, Eds.), Cambridge U.P., USA 1993, 17–65. [6] J. Douglas Jr., J.E. Gunn -Alternating direction methods for parabolic systems in mspace variables, J. Assoc. Comput. Mach., 9(1962), 450–456. [7] J. Douglas Jr., H.H. Rachford Jr. -On the numerical solution of heat conduction problems in two and three space variables, Trans. Amer. Math. Soc., 82 (1956), 421–439. [8] C. Fletcher -Computational techniques for Fluids Dinamics. Vol 2. Specific techniques for different flow categories, Springer-Verlag, Berlin (1991). [9] J.R. Galo -Resoluci´on en paralelo de las ecuaciones de Navier-Stokes 2D y 3D mediante el m´etodo de direcciones simult´aneas, Tesis, Universidad de Sevilla (2002). [10] R. Glowinski -Supercomputing and the finite element approximation of the Navier-Stokes equations for incompressible viscous fluids, Recent advances in computational fluid dynamics (Princeton, NJ, 1988), Springer, Berlin 1989, 277–315. [11] S.L. Johnsson, Y. Saad, M.H. Schultz -Alternating direction methods on multiprocessors, SIAM J. Sci. Statist. Comput., 8 (5) (1987), 686–700. [12] T. Lu, P. Neittaanmaki, X.C. Tai -A parallel splitting-up method for partial differential equations and its applications to Navier-Stokes equations, RAIRO Mod´el. Math. Anal. Num´er., Vol. 26 (6) (1992), 673– 708. [13] G.I. Marchuk -Splitting and alternating direction methods, Handbook of numerical analysis, Vol. I, North-Holland, Amsterdam 1990, 197–462. [14] J.T. Oden, O. Jacquote -Stability of some mixed finite element methods for Stokesian flows, Comput. Methods Appl. Mech. Eng., 43(2) (1984), 231–247. [15] D.W. Peaceman, H.H. Rachford Jr. -The numerical solution of parabolic and elliptic differential equations, J. Soc. Indust. Appl. Math., 3(1955), 28–41. 52 M. C. Calzada y otros [16] G. Strang -Accurate partial difference methods. I. Non-linear problems, Numer. Math., 6(1964), 37–46. [17] X.C. Tai, P. ˜ Neittaanmaki -Parallel finite element splitting-up method for parabolic problems,7 (3) (1991), 209–225. [18] R.S. Varga -Matrix iterative analysis, Springer-Verlag, Berlin 2000. [19] B.R. Wetton -Analysis of the spatial error for a class of finite difference methods for viscous incompressible flow, SIAM J. Numer. Anal., 2(1997), 723–755.