1. INTRODUCCIÓN
Es cada vez más frecuente la definición de nuevos métodos a partir de modificaciones del método clásico de Diferencias Finitas (DF). Uno de estos métodos son los esquemas conocido como métodos miméticos [1, 2, 3]. Dentro de estos esquemas están los desarrollados por Castillo et al. en [4, 5] que además de resultar conservativos logran preservar el orden de convergencia en todos los nodos del dominio discreto. Estos se conocen como esquemas de Diferencias Miméticas (DM) y han mostrado ventajas importantes sobre los esquemas clásicos de diferencias finitas en la resolución de diversos problemas que surgen en la ingeniería y las ciencias [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. No obstante, poco se ha dicho sobre cómo definir un proceso adaptativo que defina una malla que resulte óptima para encontrar una solución que cumpla con cierta precisión o tolerancia de error que imponga el problema a resolver.
Vale resaltar que los procesos adaptativos en Diferencias Finitas, y por tanto en Diferencias Miméticas, resultan poco atractivos y con poca efectividad si son comparados con los tradicionales procesos adaptativos llevados a cabo en el Método de Elementos Finitos o Volúmenes Finitos. Esto se debe principalmente al carácter estructural que deben seguir las mallas en Diferencias Finitas. No obstante, y pese a lo antes dicho, en muchos problemas (por ejemplo, frentes de choque en dinámica de fluidos) puede resultar conveniente ajustar la malla para lograr una mejor precisión con muchos menos grados de libertad que si se trabaja con una malla de paso uniforme, principalmente si se está hablando de problemas multidimensionales.
Dentro del contexto de los procesos adaptativos, un primer paso fue dado por [17] al construir los operadores gradiente y divergente sobre mallas estructuradas no uniformes, sugiriendo un proceso adaptativo usando equidistribución a partir de una función de mallado (meshsize function). Por otro lado, en [18] se define un proceso adaptativo (tipo h) para la ecuación estacionaria de convección-difusión, usando la aproximación del gradiente para la estimación del error en la variable espacial. Partiendo de estas ideas, el resultado principal de este trabajo, consiste en definir una estimación del error local (indicador de error) e implementar un procedimiento adaptativo que define una malla óptima para calcular la solución de problemas transitorios parabólicos unidimensionales, utilizando los operadores discretos no uniformes propuestos por [17].
Para definir una estimación del error cometido en la aproximación mimética (espacial), se utiliza la discretización mimética del operador gradiente sin hacer referencia a la solución analítica del problema de contorno, la cual es generalmente desconocida. En otras palabras, se hace una estimación de error en la derivada de la solución, y no en la solución. Este cálculo no amerita trabajo adicional alguno, dado que el gradiente ha sido definido previamente para calcular la solución del problema. Adicionalmente, se puede decir que el proceso propuesto para la estimación del error sigue las ideas del estimador residual SPR (superconvergent patch recovery) introducido por [19] para el suavizado de tensiones (derivadas), y por [20] en el caso de desplazamientos (incógnita del problema); ambas versiones antes implementadas para el caso del Método de los Elementos Finitos. El resto del artículo se estructura de la siguiente manera. En el siguiente apartado, se describe brevemente el método de Diferencias Miméticas utilizado y se introduce el problema de valor de frontera modelo junto a sus condiciones de frontera (tipo Robin). En el tercer apartado, se introduce el estimador de error; posteriormente, en la cuarta sección, el proceso adaptativo. Finalmente, se presenta la experimentación numérica y la discusión de los resultados obtenidos junto a las referencias del trabajo.
2. DIFERENCIAS MIMÉTICAS
Las Diferencias Miméticas se basan en la discretización de los operadores clásicos de las ecuaciones en derivadas parciales (EDP) (divergencia, gradiente y rotacional) de tal forma que ellos satisfagan una versión discreta del Teorema de Stokes o identidad de Green [5]:
En esta expresión, D, G y B son las versiones discretas de sus continuos correspondientes: gradiente (∇), divergencia (∇·) y operador de frontera ∂/∂n. Los h i representan un producto interior generalizado con pesos Q, P y I. Usando la identidad (1), se obtiene una relación para el operador de frontera B = QD + GtP.
Para la discretización espacial, se define una malla, no necesariamente uniforme, cuya geometría está dada por los nodos xi, con i = 0,1,..., N, y las celdas de la misma son los intervalos [xi−1, xi], con tamaño Hi := xi − xi−1 para i = 1,..., N. Los nodos intermedios de las celdas quedan dados por xi+1/2 = (xi + xi+1)/2, y la longitud entre dos nodos intermedios consecutivos es dada porJxi := xi+1/2 − xi−1/2, para i = 1,..., N − 1. La celda [xi−1, xi] será referida como la celda o elemento Ωi. La solución y el operador divergencia se definen en el centro de las celdas, mientras el operador gradiente en los nodos xi que definen las celdas (ver Figura 1).
El operador gradiente discreto de segundo orden, introducido por [17] para mallas no uniformes 1D, queda dado por (2) donde JG es una matriz diagonal con las inversas de las longitudes Jxi ,
y Gˆ es la parte fija del operador gradiente para mallas uniformes [5]. Por otro lado, el operador divergencia de segundo orden, para mallas no uniformes 1D, está dado por (3) donde JD es una matriz diagonal con las inversas
de las longitudes de las celdas y Dˆ la parte fija del operador divergencia para mallas uniformes. Si se implementa (2) y (3) sobre mallas uniformes, entonces D y G se reducen a los operadores presentados en [5]. El operador de frontera B se puede ver en [11] o en [6] donde se define un esquema de Diferencias Miméticas conservativas equivalente.
2.1. Problema modelo transitorio
El problema de frontera a resolver viene dado por la ecuación que modela la transferencia de calor por conveccióndifusión en su forma conservativa
donde Ω = (a,b), u(x,t) representa la temperatura (variable del problema) en un punto (x,t) del dominio Ω y un tiempo t, ν > 0 el coeficiente de transmisión térmica, y f una función escalar que describe la existencia de una fuente o sumidero en el problema. La ecuación se completa al definir la condición inicial u(x,0) = g(x) y las condiciones de frontera:
con {αi,βi,γi}, i = a,b, parámetros reales conocidos y dependiendo de su valor, determinan condiciones de contorno de Dirichlet, Neumann o Robin. El problema de contorno (4)-(5) es bien planteado en los casos de las condiciones de borde consideradas en este artículo. Para efectos de evaluar la precisión y eficiencia del proceso adaptativo, se implementarán casos en que se pueda obtener la solución analítica del problema.
2.2. Discretización mimética para el problema modelo
Usando un esquema implícito para la discretización temporal con un paso de tiempo k, más las discretizaciones de los operadores, la aproximación mimética para la ecuación de convección-difusión (4) queda dada por
representa la solución aproximada mimética en el tiempo tn = nk de la solución exacta u del problema, y
representa la restricción de f a la malla mimética. Como el operador divergencia discretizado no actúa sobre la frontera, entonces las condiciones de contorno Robin (5) se obtienen a partir de
donde el vector fb resulta de restringir el término no homogéneo de las condiciones de contorno a la malla mimética, es decir, fb = (γa,0,...,0,γb)t. Las matrices [α] y [β] son tales que α1,1 = αa, αN+2,N+2 = αb, β1,1 = βa, βN+2,N+2 = βb, y el resto de entradas son cero. A partir de (6) y (7), el esquema mimético para la ecuación de convección-difusión (4) sujeto a las condiciones de contorno Robin (5) queda dado por
3. ESTIMACIÓN DEL ERROR
Al resolver numéricamente un problema, resulta conveniente controlar la calidad de la solución aproximada o el gradiente de tal solución. Este tipo de control se puede lograr mediante un proceso de adaptatividad de la malla para aproximar eficientemente la solución del problema. Igualmente, en un proceso adaptativo resulta indispensable disponer, en todo el dominio, de la distribución local del error que se comete al usar la solución aproximada como solución del modelo. En otras palabras, es necesario estimar el error cometido en la aproximación para cada celda del dominio. En nuestro caso, la estimación del error en la variable espacial se hace sobre el gradiente de la solución mimética, es decir, se debe estimar el error que surge al usar GU para aproximar el gradiente exacto ∇u. El error numérico en la aproximación GU de ∇u será denotado por eG := ∇u − GU.
3.1. Estimación del error espacial usando el suavizado del gradiente
El proceso de estimación a posteriori del
error que se presenta proporciona una aproximación z, obtenida realizando un postproceso sobre la solución GU. La
meta es obtener el error estimado e
GU
en cada celda del dominio Ω, tal que e∗G ≈ eG .
Toda celda, Ωi, de la malla tiene asociada una parcela de celdas (del inglés patch), ωi, constituida por todas las celdas que la rodean (ver Figura 2). Para el caso de elementos de frontera, la parcela de estas celdas queda dada por el elemento de frontera y los dos elementos consecutivos (anteriores) a éste. El cálculo de z se realiza localmente para cada parcela de celdas que se forme en la malla. Los valores de GU en los nodos de la malla son usados como entrada para definir un polinomio de mayor orden (interpolación). Para problemas unidimensionales, la parcela de elementos involucra 4 nodos de la malla (ver Figura 2). En este caso, se define un polinomio cúbico usando interpolación sobre los cuatro valores de GU. En particular, para cualquier nodo xj interno en la celda Ωi se puede definir z(xj) en dicha celda. Es decir, la estimación del error restringida a la celda Ωi y evaluada en el nodo xj está dada por
donde IhGU representa la interpolación lineal de GU evaluada en el nodo xj, y z∗j el valor de z∗ en el nodo xj.
El suavizado propuesto sobre el GU es capaz de mejorar la curvatura de la solución aproximada, es decir, la solución suavizada z∗ mejora en sus derivadas. Sin embargo, esta conclusión no necesariamente ocurre para los valores sin derivar (solución U del problema). Este hecho justifica nuestra selección, para la estimación del error, en el gradiente y no en la solución U del problema.
4. CONTROL ADAPTATIVO DEL ERROR
La estrategia más simple para controlar el error cometido al usar la solución aproximada, e := u − U, consiste en un proceso iterativo que implemente la técnica de estimación dada en el apartado anterior para definir la distribución de los nodos (malla) en el paso de tiempo actual, y mantener dicha malla en los siguientes pasos de tiempo hasta que la tolerancia del error no se satisfaga y sea necesario redefinir la malla. Esta estrategia resulta muy general y no requiere información del tipo de problema a tratar. Los principales pasos que conforman el algoritmo son dados en el Cuadro 1.
Para el criterio de parada del proceso adaptativo (Cuadro 1, Paso 3), se pide que kek ≤ TOL. No obstante, la estimación del error se hace sobre el gradiente de la solución y no se tiene a la mano un valor estimado, e∗, para el error e. Para solventar esta dificultad, se usa para cada celda del dominio
donde los ηj representan los pesos y los xj dados en xˆj = 0.5 |Ωi|xj + (xj + xj−1) los puntos de la cuadratura de Gauss (en nuestro caso basta tomar npG = 2).
Continuando el proceso para todas las celdas Ωi, se logra definir el error estimado global ke∗k = Pi ke∗ki. En la práctica, se puede usar el error estimado del gradiente, ke∗k, en lugar de la estimación ke∗k del error en
la solución. Sin embargo, este cambio puede disminuir la eficiencia del proceso adaptativo, llegando a definir mallas con muchos más nodos de los necesarios, pero con un comportamiento similar a la alcanzada al usar ke∗k. Por otro lado, la cota de error local, calculada en el paso (3c), indica la forma a distribuir el error permitido (TOL) en la aproximación del problema. En este trabajo, se define una distribución uniforme del error. Sin embargo, otros criterios pueden definirse en lugar del dado acá, por ejemplo, una distribución local del error proporcional al tamaño de elemento.
5. RESULTADOS NUMÉRICOS
En este apartado, se presentan los resultados numéricos obtenidos al aplicar el algoritmo adaptativo propuesto. El carácter sintético del ejemplo analizado busca justificar la efectividad del estimador en casos donde el problema presenta fuertes variaciones en su derivadas.
La calidad de la estimación, e∗, es medida usando el co-global o índice de efectividad local, en el caso que se mida sobre una celda específica. Estos índices pueden ser usados para medir la calidad de la estimación cuando se conoce el error exacto (que es el caso del ejemplo analizado); en caso contrario, una buena aproximación del mismo es requerida. En la práctica, y de forma general, para medir Ief f se usa un error de referencia eref en lugar
G de eG, definido sobre una malla fina, de tamaño h, que puede estar dada a partir de h = H/nref. Para el análisis del error se usan las normas del máximo, k · k∞, y la norma L2 discreta, k · k2, definidas como
Ejemplo 1
Se resuelve el problema modelo (4)-(5) para un coeficiente de transmisión térmica ν = 1 y un término fuente f(x,t) definido de tal manera que la solución analítica del problema viene dada como u(x,t) = (1− x) arctan(α(x − t)) + arctan(αt). Para este ejemplo, u(0,t) = u(1,t) = 0, lo cual simplifica las condiciones de frontera:
con α = 250. El intervalo de tiempo es t ∈ [0,0.9] y un paso de tiempo dado por dt = 0.002. Además, no se avanza en el tiempo hasta que se cumpla una tolerancia de error del 0.001 para el gradiente ke∗k. La Figu-
G ra 3 muestra la solución aproximada (línea con •) junto con la solución analítica (línea continua) para el tiempo t = 0.5. En la parte inferior se muestran las mallas obtenidas (distribución espacial de los nodos) durante el proceso adaptativo, para distintos pasos de tiempo; la mayor concentración de nodos ocurre entre los valores donde el gradiente de la solución es mas pronunciado (tal como se podía esperar). La Figura 4 ilustra el número de nodos necesarios en cada paso temporal, El proceso adaptativo se inicia en una malla uniforme de 10 celdas (11 nodos) y un error relativo en la solución cerca del 30%. El error relativo en el gradiente está alrededor del 900%, y se quiere alcanzar una solución en cada paso de tiempo con un error inferior al 0.1%. Los resultados numéricos para la sucesión de mallas obtenidas en algunos pasos de tiempo, para una tolerancia de error de 10−3, se presentan en el Cuadro 2. En este cuadro, el error relativo exacto en el gradiente, keGk/kGuk,
cuarta columna, respectivamente. El error relativo en la aproximación al usar la solución mimética, kek/kuk, es dado en la quinta columna. En la última columna, se presenta el índice de efectividad global. Los resultados son mostrados solamente en norma del máximo, k · k∞, esto es debido a que los resultados obtenidos en norma k · k2 resultan, en general, análogos. En cada paso de tiempo mostrado en el Cuadro 2, se dan los valores obtenidos
índice de efectividad global.
para la ultima malla que define el proceso adaptativo. El hecho que no se alcance la tolerancia prescrita, en cada paso de tiempo, en la primera iteración del proceso iterativo, no implica que el proceso adaptativo esté funcionando de manera errática o su efectividad sea baja, esto es debido principalmente a las condiciones adicionales que se imponen al proceso, por ejemplo, porcentaje de celdas a ser divididas, número de divisiones por celda.
"Verdana" size="2">
Gráficamente, estos resultados se muestran en las Figuras 5 y 6 para un tiempo t = 0.3. Los errores relativos estimados y exactos (gradiente) en norma del m�ximo y L2 son mostrados en la Figura 5 y 6, respectivamente. La tendencia asintótica del error estimado, al ir ajustando la malla, refleja el buen comportamiento del estimador de error propuesto.
Para definir el rendimiento del proceso adaptativo, se comparan los resultados con una malla uniforme de 500 elementos para todos los pasos de tiempo. El comportamiento del error para la malla uniforme y los obtenidos del proceso adaptativo resultan equivalente (ver Figura 7). Sin embargo, cálculos hechos con mallas uniformes con menos de 500 elementos presentan un comportamiento oscilatorio en el error (la Figura 7 muestra el caso de 300 elementos), tendiendo a perder la estabilidad a medida que se avanza en tiempo (por encima de la condición de estabilidad del esquema implícito). adaptativo.
6. CONCLUSIONES Y COMENTARIOS FINALES
Se ha establecido un proceso h-adaptativo a partir de una estimación del error por suavizado en el gradiente para la variable espacial. La malla se mantiene fija por los sucesivos pasos de tiempo siempre y cuando cumpla la tolerancia prescrita del error. La experimentación numérica verifica el buen comportamiento del estimador de error por postproceso propuesto, proporcionando mallas adaptadas que aceleran la precisión (convergencia) de la solución aproximada en comparación con mallas refinadas de forma uniforme.
Es conocido que no existe literatura que defina procesos adaptativos para el control del error de soluciones numéricas dadas por las Diferencias Miméticos propuestos por [5], especialmente en el caso de problemas transitorios. La generalización de los resultados mostrados en este trabajo a problemas 2D, tanto estacionarios como no estacionarias, se puede intuir con cierta facilidad, por lo menos si pensamos en mallas estructuradas no uniformes. En el presente, los autores trabajan en dicha generalización.
7. AGRADECIMIENTOS
Giovanni Calderón and Julio C. Carrillo E. desean reconocer el apoyo financiero proporcionado por la VIE-UIS en virtud del proyecto N◦ 2415. Jorge Villamizar agradece el apoyo proporcionado por VIE-UIS bajo el Programa de Movilidad N◦ 2303.
8. REFERENCIAS
[1] M. Shashkov and S. Steinberg, “Support-operator finitedifference algorithms for general elliptic problems,” Journal of Computational Physics, vol. 118, no. 1, pp. 131–151, 1995.
[2] J. M. Hyman and M. Shashkov, “The approximation of boundary conditions for mimetic finite difference methods,” Computers Math. Applic., vol. 36, no. 5, pp. 79–99, 1998.
[3] J. M. Hyman, M. Shashkov, and S. Steinberg, “Mimetic finite difference methods for diffusion equations,” Computers Math. Applic., vol. 6, no. 3-4, pp. 333–352, 2002.
[4] J. Castillo and R. D. Grone, “A matrix analysis approach to higher-order approximations for divergence and gradients satisfying a global conservation law,” SIAM J. Matrix Anal. Appl., vol. 25, no. 1, pp. 128–142, 2003.
[5] J. Castillo and M. Yasuda, “Linear systems arising for secondorder mimetic divergence and gradient discretizations,” Journal of Mathematical Modelling and Algorithms, vol. 4, no. 1, pp. 67– 82, 2005.
[6] J. Guevara, M. Freites, and J. Castillo, “A new second order finite difference conservative scheme,” Divulgaciones Matemáticas, vol. 13, no. 1, pp. 107–122, 2005.
[7] F. Solano-Feo, J. Guevara-Jordan, O. Rojas, B. Otero, and R. Rodriguez, “A new mimetic scheme for the acoustic wave equation,” Journal of Computational and Applied Mathematics, vol. 295, no. 1, pp. 2–12, 2016.
[8] G. Sosa Jones, J. Arteaga, and O. Jiménez, “A study of mimetic and finite difference methods for the static diffusion equation,” Computers and Mathematics with Applications, vol. 76, no. 1, pp. 633–648, 2018.
[9] O. Montilla, C. Cadenas, and J. Castillo, “Matrix approach to mimetic discretizations for differential operators on non-uniform grids,” Mathematics and Computers in Simulation, vol. 73, no. 1-4, pp. 215–225, 2006.
[10] J. M. Guevara-Jordan, S. Rojas, M. Freites-Villegas, and J. E. Castillo, “Convergence of a mimetic finite difference method for static diffusion equation,” Advances in Difference Equations, vol. 2007, no. Article ID 12303, 12 pages, 2007.
[11] J. Castillo and G. Miranda, Mimetic Discretization Methods, 5th ed. Boca Raton, Florida: CRC Press, 2013.
[12] F. Hernández, J. Castillo, and G. Larrazábal, “Large sparse linear systems arising from mimetic discretization,” Computers and Mathematics with Applications, vol. 53, no. 1-11, pp. 215– 225, 2007.
[13] C. Bazan, M. Abouali, J. Castillo, and P. Blomgren, “Mimetic finite difference methods in image processing,” Computational and Applied Mathematics, vol. 30, no. 3, pp. 701–720, 2011.
[14] M. Abouali and J. Castillo, “Solving navier-stokes equation using castillo-grone’s mimetic difference operators on gpus,” APS Division of Fluid Dynamics, vol. 57, no. 17, pp. 701–720, 2012.
[15] A. Gomez-Polanco, G.-J. J. M., and B. Molina, “A mimetic iterative scheme for solving biharmonic equations,” Mathematical and Computer Modelling, vol. 57, no. 9-10, pp. 2132–2139, 2013.
[16] J. Blanco, O. Rojas, C. Chacón, J. M. Guevara-Jordan, and J. Castillo, “Tensor formulation of 3-D mimetic finite differences and applications to elliptic problems,” Electron. Trans. Numer. Anal., vol. 45, pp. 457–475, 2016.
[17] E. Batista and J. Castillo, “Mimetic schemes on non-uniform structured meshes,” Electronic Transactions on Numerical Analysis, vol. 34, no. 1, pp. 152–162, 2009.
[18] G. Calderón and A. Lugo, “Error estimation and adaptivity in mimetic schemes for boundary value problems,” Boletín Venezolano de Matemáticas, no. 2, pp. 109–124, 2015.
[19] O. Zienkiewicz and J. Zhu, “The superconvergent patch recovery (SPR) and a posteriori error estimates. Part 1: The recovery technique,” Int. J. Num. Meth. Engrg., vol. 33, pp. 1331–1364, 1992.
[20] P. Díez and G. Calderón, “Goal-oriented error estimation for transient parabolic problems,” Computational Mechanics, vol. 35, no. 5, pp. 631–646, 2007.