Aumento de la eficiencia de un método de descomposición de dominio mediante estimaciones a posteriori

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Aumento de la eficiencia de un m´etodo de descomposici´on de dominio mediante estimaciones a posteriori. C. Bernardi 1T. Chac´ on Rebollo2, E. Chac´ on Vera2, D. Franco Coronil2, 1Lab. Jacques-Louis Lions, C. N. R. S. et Univ. Pierre et Marie Curie, Boˆıte courrier 187, 75252 Paris Cedex 05 (France). E-mail: [email protected]. 2Dpto. E.D.A.N., Universidad de Sevilla, Aptdo. 1160, E-41080 Sevilla. E-mails: [email protected], [email protected], [email protected]. Palabras clave: descomposici´on de dominios, an´alisis de error a posteriori, penalizaci´on, normas fraccionarias de Sobolev Resumen En este trabajo introducimos un m´etodo de descomposici´on de dominio sin solapamiento con penalizaci´on, que viene motivado a partir de un an´alisis del error a posteriori del m´etodo estudiado por T. Chac´on y E. Chac´on en [5] y [6]. Con el objetivo de mejorar la tasa de convergencia del m´etodo de [6], en este trabajo, introducimos una nueva versi´on de este m´etodo en la cual un t´ermino de penalizaci´on H1/2 00 (Γ) reemplaza el t´ermino L2(Γ) del original de [6]. Usando este nuevo t´ermino, el n´umero de iteraciones necesarias para alcanzar una soluci´on con un error del mismo orden que el error de discretizaci´on, se reduce significativamente. Realizamos adem´as un an´alisis de error a posteriori, que nos permite desarrollar un estrategia para determinar simult´aneamente un par´ametro de penalizaci´on ´optimo y una malla optimal, para reducir el error por debajo de un valor prefijado. Varios test num´ericos muestran los buenos resultados de nuestras aproximaciones. 1. El problema penalizado Sea Ω ⊂Rd(d= 2,3) un dominio simplemente conexo y acotado con una frontera ∂Ω Lipschitz-continua. Consideramos una descomposici´on simple de Ω en dos subdominios Ω1y Ω2que no se solapan. Sean Γ = ∂Ω1∩∂Ω2, Γi=∂Ωi∩∂Ω, i= 1,2, fronteras que suponemos que todas son dominios (d−1)-dimensionales Lipschitz-continuos, con medida (d−1)-dimensional positiva. Denotamos por nij el vector normal exterior a Γ que apunta desde Ωihacia Ωjy sea n=n12. 1 C. Bernardi, T. Chac´on, E. Chac´on, D. Franco Consideramos adem´as los espacios de Sobolev Xi=H1(Ωi; Γi) = {v∈H1(Ωi) tal que v|Γi= 0}, i = 1,2; X=X1×X2. Para u= (u1, u2),v= (v1, v2)∈Xdefinimos el producto escalar y la norma sobre X ((u,v))X= 2 X i=1 (∇ui,∇vi)Ωi,kuk2 X= ((u,u))X, Consideramos el problema de Poisson en Ω con condiciones de contorno homog´eneas: Dada f∈L2(Ω), hallar u∈H1 0(Ω) tal que (∇u, ∇v)Ω= (f, v)Ω,∀v∈H1 0(Ω) (1) Para introducir la nueva versi´on del m´etodo estudiado por T. Chac´on y E. Chac´on en [6], recordemos que H1/2 00 (Γ) es el subespacio de funciones de H1/2(Γ) cuya extensi´on por cero a, por ejemplo, ∂Ω1pertenece H1/2(∂Ω1). Un producto escalar intr´ınseco en H1/2 00 (Γ) viene definido por [[w, v]]Γ=ZΓ w(x)v(x)dx +ZΓZΓ (w(x)−w(y)) (v(x)−v(y)) |x−y|ddx dy (2) +ZΓ w(x)v(x) d(x, ∂Γ) dx , donde el primer usando es el producto escalar usual en L2(Γ) (Cf. [1]). Por simplificar, denotamos tambi´en por [[·,·]]Γel producto escalar L2(Γ), y estudiamos a la vez tanto la penalizaci´on L2(Γ) como la H1/2 00 (Γ) utilizando la misma notaci´on. Distinguiremos cada caso cuando sea preciso. Introducimos entonces nuestro problema penalizado con penalizaci´on H1/2 00 (Γ) ´o L2(Γ), como (P²)     Hallar u²∈Xtal que ((u²,v))X+1 ²[[u² 1−u² 2, v1−v2]]Γ= 2 X i=1 (f, vi)Ωi,para todo v∈X. con ²un par´ametro destinado a tender a cero. Este problema tiene una ´unica para cada ² > 0 debido al Lema de Lax-Milgram. Observaci´on 1 El m´etodo original introducido en [6] consideraba s´olo un t´ermino de penalizaci´on L2(Γ) para asegurar la continuidad de u²= (u² 1, u² 2)a trav´es de Γ. Nosotros aqu´ı consideramos adem´as un t´ermino de penalizaci´on H1/2 00 (Γ) para garantizar la continuidad en un sentido m´as fuerte. En ambos casos, el m´etodo puede ser interpretado como la formulaci´on variacional de un sistema acoplado de PDEs con la siguiente estructura        −∆u1=fen Ω1, u1= 0 sobre Γ1, ∂n12 u1=1 ²b(u1−u2)sobre Γ,        −∆u2=fen Ω2, u2= 0 sobre Γ2, ∂n21 u2=1 ²b(u2−u1)sobre Γ, donde bes un operador lineal, acotado e inyectivo definido sobre Γ, que se reduce a la identidad cuando se usa penalizaci´on L2(Γ). En consecuencia, el m´etodo asegura la continuidad de los flujos normales a trav´es de Γy fuerza por penalizaci´on la continuidad de u². 2 Aumento de la eficiencia de un DDM 2. Discretizaci´on del problema penalizado Para discretizar el problema (P²), consideramos una familia regular de triangulaciones {Tih}{h>0}de cada Ωitales que Γ sea la uni´on de caras o lados completos de elementos K en cada Tih. Sea Th=T1h∪T2h. Suponemos que los conjuntos de trazas de T1hyT2hsobre Γ son iguales a un mismo conjunto, que denotamos por EΓ h. Entonces, {Th}{h>0}constituye una familia de triangulaciones regulares de Ω. Ahora sobre las triangulaciones Tih construimos una familia de subespacios de elementos finitos Xih de Xi(i= 1,2). Sea Xh=X1h×X2h. Entonces, nuestro problema penalizado discretizado es la aproximaci´on por elementos finitos Galerkin estandar u² h= (u² 1h, u² 2h)∈Xhde u², soluci´on de (P²,h)½Hallar u² h∈Xhtal que ((u² h,vh))²=F(vh) para todo vh∈Xh. Este problema admite una ´unica soluci´on para cada ² > 0. Para resolver num´ericamente (P²,h), consideramos ahora la t´ecnica paralela mediante aproximaciones sucesivas introducida en [6] siguiente: Para n= 0,1,2, ..., los t´erminos un+1 1=un+1,² 1h∈X1hyun+1 2=un+1,² 2h∈X2hse computan a partir de un 1hyun 2hresolviendo (quitamos los ´ındices ²yhpor simplificar), (Pn ²,h)            (∇un+1 1,∇v1h)Ω1+1 ²[[un+1 1−un 2, v1h]]Γ= (f, v1h)Ω1,∀v1h∈X1h, (∇un+1 2,∇v2h)Ω2+1 ²[[un+1 2−un 1, v2h]]Γ= (f, v2h)Ω2,∀v2h∈X2h. (3) Para el an´alisis que realizaremos en las secciones siguientes, supondremos que cada espacio de discretizaci´on Xih contiene el espacio X∗ ih definido por X∗ ih ={wh∈H1(Ωi) ; wh|K∈Pki(K),∀K∈ Tih, wh|Γi = 0 },(4) donde Pki(K) denota el espacio de las restricciones a Kde los polinomios con dvariables y grado total ≤kipara kientero positivo. Adem´as asumiremos que su espacio de trazas sobre Γ es Wih ={vh∈H1 0(Γ) : vh|e∈Pki(e),∀e∈ EΓ h}.(5) 3. An´alisis de error a posteriori En esta secci´on realizamos un an´alisis de error a posteriori de ambos m´etodos, con el objetivo de estudiar la optimalidad tanto de un indicador del error de penalizaci´on como de un indicador del error de discretizaci´on, que definimos como sigue: Una familia de indicadores locales del error de discretizaci´on, es: Para cada i= 1,2 y K∈ Tih, ηK i=hKkfh+ ∆u² ihkL2(K)+X e∈EK h1/2 ek[∂nKu² ih]kL2(e), 3 C. Bernardi, T. Chac´on, E. Chac´on, D. Franco donde •hKyheson el di´ametro de Kye, respectivamente, •nKes el vector normal unitario exterior a ∂K, y •[∂nKu² ih] es el salto de ∂nKu² ih a trav´es esi eno est´a incluido en Γ, y ∂nKu² 1h−∂nKu² 2hsi eest´a incluido en Γ. • EKes el conjunto de los lados (d= 2) o caras (d= 3) de Kque no est´an contenidas en Γ. Entonces, un indicador del error de discretizaci´on es la suma hilbertiana de los ηK i: ηD h=  2 X i=1 X K∈Tih ¡ηK i¢2  1/2 . Un indicador del error de penalizaci´on, es: ηP h=ku² 1h−u² 2hkH1/2 00 (Γ). Para analizar la optimalidad de ambos indicadores de error, consideramos el error de penalizaci´on E²y el error de discretizaci´on Eh, definidos respectivamente por E²= 2 X i=1 |u−u² i|1,Ωi, Eh= 2 X i=1 |u² i−u² ih|1,Ωi. Entonces, en primer lugar, establecemos las siguientes acotaciones: Teorema 1 Los indicadores de error ηP hyηD h, verifican las siguientes estimaciones ηD h≤C     Eh+  2 X i=1 X K∈Tih h2 K||f−fh||2 L2(K)  1/2     ,(6) ηP h≤C0 2 X i=1 |u−u² ih|1,Ωi≤C0(Eh+E²),(7) para ciertas constantes C > 0yC0>0que son independientes de ²yh. La demostraci´on puede verse en [2]. Nuestro resultado principal de estimaci´on de error a posteriori, que prueba la cuasi optimalidad de ambos indicadores, es el siguiente: Teorema 2 Supongamos que los espacios de trazas W1hyW2hdefinidos en (5) coinciden. Entonces, la soluci´on del problema penalizado discreto (Ph,²)satisface la siguiente estimaci´on de error a posteriori 2 X i=1 |u−u² ih|1,Ωi≤C     ηP h+  2 X i=1 X K∈Tih ¡µKηK i¢2+h2 K||f−fh||2 L2(K)  1/2     , 4 Aumento de la eficiencia de un DDM con µK=½(1 + λh)si K∩Γ6=∅, 1si K∩Γ = ∅., donde i) En el caso de penalizaci´on L2(Γ) es, λh= (hM/hm)1/2,con hM= m´ax K∩Γ6=∅hK, hm= m´ın K∩Γ6=∅hK; ii) En el caso de penalizaci´on H1/2 00 (Γ) es, λh= (hM/hm+h−1 m)1/2. La demostraci´on tambi´en puede verse en [2]. Observaci´on 2 Los teoremas (1) y (2) prueban la cuasi optimalidad de nuestros indicadores de error. La optimalidad se tendr´ıa si µK= 1, pero no la hemos podido probar por dificultades t´ecnicas. En el caso de penalizaci´on H1/2 00 (Γ), hay adem´as una perdida de optimalidad de h−1/2 m. No obstante, esta perdida se limita tambi´en a los elementos, lados o caras que intersectan a la interface Γ. 4. Experimentos num´ericos Para analizar el rendimiento pr´actico de nuestros m´etodos, hemos realizado varios ensayos num´ericos para el problema de Poisson como problema modelo, en dimensi´on d= 2. Estos ensayos han sido desarrollados sobre el c´odigo de elementos finitos FreeFEM++, usando elementos finitos P1-Lagrange. Cada experimento num´erico ha sido validado en el caso del dominio Ω la L-shape ]0,1[2\[1/2,1[2, descompuesto en dos subdominios, con interface Γ su intersecci´on con la recta y=x. Adem´as, hemos considerado un problema de Poisson con condiciones de contorno homog´eneas con soluci´on anal´ıtica conocida, definido sobre la L-shape. Finalmente, debido a la dificultad para computar el producto escalar H1/2 00 (Γ) definido en (2), en los experimentos num´ericos hemos utilizado un producto escalar que lo discretiza, construido usando f´ormulas de cuadratura (Cf. [4]). 4.1. Mejora de la tasa de convergencia: Observaci´on 3 Cuando usamos penalizaci´on H1/2 00 (Γ), es necesario un n´umero de iteraciones de orden O(|log h|h−k)para obtener un error O(hk). Por contra, en [6] se demuestra que si se usa penalizaci´on L2(Γ) son necesarias del orden de O(|log h|h−2k)iteraciones para obtener un error del mismo orden (Cf. [2]). Para verificar esta mejora, hemos tenido en cuenta que, seg´un los resultados te´oricos, mientras que en el caso de penalizaci´on H1/2 00 (Γ) la estimaci´on de orden ´optimo corresponde a²=O(hk), en el caso L2(Γ) este orden corresponde a ²=O(h2k). Entonces, hemos 5 C. Bernardi, T. Chac´on, E. Chac´on, D. Franco 0.025 0.05 0.1 0.2 100 500 1000 2500 LShape test: Influence of the penalty parameter in the number of iterations (h=1/64) Epsilon No Iter (Pen. H0012) No Iter (Pen. L2) Figura 1: Influencia del par´ametro de penalizaci´on sobre el n´umero de iteraciones. considerado una malla fija uniforme y hemos hecho varios ensayos variando el par´ametro de penalizaci´on entre ²= 0,2 y ²= 0,025 en el caso H1/2 00 (Γ), y entre ²= 0,22y²= 0,0252 en el caso L2(Γ). Hemos calculado para cada ensayo el n´umero de iteraciones necesarias para obtener el error relativo entre la soluci´on exacta y la obtenida por nuestro c´odigo en seminorma H1(Ω). Dibujamos, en el eje de abcisas, los valores de ²para el caso H1/2 00 (Γ), y los de √²en el caso L2(Γ). La figura 1 muestra el test para h= 1/64. Observamos que para valores de ²menores que un valor ²0≃0,04, el n´umero de iteraciones necesarias para alcanzar un determinado error es notablemente inferior cuando se usa penalizaci´on H1/2 00 (Γ) (linea azul) que cuando se usa penalizaci´on L2(Γ) (linea roja). Resultados similares se obtienen para mallas de otras tallas h. 4.2. Eficiencia de los indicadores de error: Para ello comparamos los indicadores de error ηP hyηD h, con los errores 2 X i=1 |Ihu−u² ih|H1(Ωi)y 2 X i=1 kIhu−u² ihkL2(Ωi), siendo Ihuel interpolado P1-Lagrange de la soluci´on exacta u. Por un lado, para verificar la eficiencia del indicador del error de penalizaci´on, hemos realizado varios tests mantenido fija la malla y haciendo variar el par´ametro de penalizaci´on entre ²= 0,005 y ²= 15. La figura 2 muestra el test para h= 1/64, en el caso de penalizaci´on H1/2 00 (Γ). Observamos que el indicador de error ηP h(trazo rojo continuo) decrece con ²hasta que el error 6 Aumento de la eficiencia de un DDM 10−2 10−1 100101 10−7 10−6 10−5 10−4 10−3 10−2 LShape test: Influence of the H00 1/2 penalty parameter (h=1/64) NormL2 EtaP SemiNormH1 EtaD Epsilon Figura 2: Eficiencia del indicador del error de penalizaci´on. 10−2 10−1 10−7 10−6 10−5 10−4 10−3 10−2 10−1 LShape test: Influence H00 1/2 of the mesh size (eps=0.01) NormL2 EtaD SemiNormH1 EtaP h Figura 3: Eficiencia del indicador del error de discretizaci´on. debido a la penalizaci´on es comparable con el error debido a la discretizaci´on. Hasta este valor, las curvas correspondientes a este indicador y a los errores (trazo verde seminorma H1y trazo azul norma L2) son paralelas. Por contra, el indicador del error de discretizaci´on ηD h(trazo rojo discontinuo) se mantiene pr´acticamente constante y totalmente independiente de ². Por otro lado, para comprobar la eficiencia del indicador del error de discretizaci´on, cuando se usa cada tipo de penalizaci´on, hemos fijado el valor del par´ametro de penalizaci´on ²y hemos hecho variar la talla de la malla entre h= 0,02 y h= 0,5, para mallas cuasi uniformes. La figura 3 muestra el test para ²= 0,01 en el caso de penalizaci´on H1/2 00 (Γ). Observamos el mismo comportamiento cualitativo que en la figura (2), intercambiando los papeles de ηD hyηP h. Los resultados obtenidos en el caso de penalizaci´on L2(Γ) son cualitativamente similares. 4.3. Optimizaci´on del par´ametro de penalizaci´on con respecto a la adaptatividad de la malla: Con el objetivo de equilibrar los errores provenientes de discretizaci´on y de penalizaci´on, consideramos la siguiente estrategia computacional, que fue propuesta por C. Bernardi et al. en [3] y consta de tres etapas. Primero elegimos una tolerancia η∗, hacemos una primera computaci´on sobre una malla cuasi uniforme y computamos ηP h. Etapa 1 Inicializaci´on: Si ηP h≤η∗entonces vamos a la Etapa 2. Si NO dividimos ²por el cociente ηP h/η∗y realizamos un nuevo c´alculo. Etapa 2 Adaptaci´on de malla: Computamos los ηK iy su valor medio ¯ηD h. Entonces, para cada Ktal que ηK ies mayor que ¯ηD h, dividimos Ken tri´angulos (o tetraedros en 7 C. Bernardi, T. Chac´on, E. Chac´on, D. Franco dimension d= 3) m´as peque˜nos, de modo que el di´ametro de estos nuevos elementos se comporta como hKmultiplicado por el cociente ¯ηD h/ηK i. Etapa 3 Penalizaci´on: Calculamos ηP hyηD h. Si Si ηP h≤ηD hentonces volvemos a la Etapa 2. Si NO dividimos ²por un n´umero constante de veces el cociente ηP h/ηD hy volvemos a la etapa 2. Hemos desarrollado varios test num´ericos, en las que dado un error prefijado, usando esta estrategia, siempre conseguimos alcanzar este error, para un valor ´optimo del par´ametro de penalizaci´on y una malla final adaptada. Observamos que, mediante esta estrategia, el tiempo de c´alculo necesario para alcanzar este error se reduce dr´asticamente con respecto al tiempo de c´alculo que necesitar´ıamos, si realiz´aramos el c´alculo partiendo de la malla final, del valor ´optimo de ²y de una soluci´on inicial u0,² h= 0. En concreto, para alcanzar un error relativo de 0.003 entre la soluci´on exacta y la obtenida por nuestro c´odigo en seminorma H1(Ω), el tiempo de c´alculo se divide aproximadamente por 4. Observaci´on 4 No obstante, el tiempo de c´alculo necesario para obtener un determinado error, es netamente mayor en el caso de usar penalizaci´on H1/2 00 (Γ) que cuando usamos penalizaci´on L2(Γ). Posiblemente, esto es debido al uso de FreeFEM++ que no permite adaptar la programaci´on. Esto conduce a que, en especial, empleemos mucho tiempo de c´alculo al calcular las normas H1/2 00 (Γ) si la malla est´a muy adaptada. Aunque la figura 1 muestra un comportamiento asint´otico mejor para el caso de penalizaci´on H1/2 00 (Γ), en cuanto al n´umero de iteraciones, los resultados anteriores indican que esta mejora s´olo se va a traducir en una ganancia en tiempo de c´alculo para valores muy peque˜nos de hy². Agradecimientos Este trabajo ha sido financiado por el Ministerio de Educaci´on y Ciencia, en el marco del proyecto MTM2006-01275 del Programa Nacional de Matem´aticas del Plan Nacional de I+D+I (2004-2007). Referencias [1] Adams, R. A., Sobolev Spaces. Pure and Applied Mathematics, Vol. 65. Academic Press, New YorkLondon, 1975. [2] Bernardi, C., Chac´on Rebollo, T., Chac´on Vera, E., Franco Coronil, D., A non-overlapping domaindecomposition method motivated by a postriori error analysis. Sometido a Math. Models and Methods in Applied Sciences. [3] Bernardi, C., Girault, V., Hecht, F., A posteriori analysis of a penalty method and application to the Stokes problem. Math. Models and Methods in Applied Sciences, 13 (2003), 1599-1628. [4] Casas, E., Raymond, J.-P., The stability in Ws,p(Γ) spaces of L2-projections on some convex sets. Numer. Funct. Anal. and Optimization, 27 (2006), 117-137. [5] Chac´on Rebollo, T., Chac´on Vera, E., A non-overlapping domain decomposition method for the Stokes equations via a penalty term on the interface. C.R. Acad. Sci. Paris, t. 334, S´erie I (2002), 1–16. [6] Chac´on Rebollo, T., Chac´on Vera, E., Study of a non-overlapping domain decomposition method: Poisson and Stokes problems. Appl. Numer. Math., 48 (2004), 169–194. 8