An IMEX-DG solver with non-conforming mesh refinement for atmospheric dynamics with rotation

Resumen Técnico: Un Solver IMEX-DG con Refinamiento de Malla No Conformes para Dinámica Atmosférica con Rotación

Enunciado del Problema
Los flujos atmosféricos se caracterizan por la interacción de múltiples escalas espaciales y temporales, que van desde circulaciones equilibradas a gran escala hasta ondas acústicas y de gravedad rápidas. En modelos atmosféricos completamente compresibles, la coexistencia de modos geostroficos lentos y modos acústicos rápidos introduce una rigidez severa, lo cual restringe fuertemente el tamaño del paso de tiempo en métodos de integración totalmente explícitos. Si bien la integración temporal Implícito-Explícita (IMEX) y la discretización espacial de Galerkin Descontinua (DG) de alto orden se han combinado con éxito para abordar la rigidez acústica en configuraciones cartesianas no rotatorias, es esencial extender estos marcos para incluir la rotación planetaria con el fin de lograr predicciones numéricas del tiempo (NWP) realistas y configuraciones de referencia estándar. La inclusión de la fuerza de Coriolis introduce complejidad adicional, particularmente en lo que respecta a la preservación del equilibrio geostrofico y la implementación eficiente de términos rotacionales dentro de solvers de alto orden libres de matrices.

Metodología
Los autores presentan una extensión de un solver IMEX-DG de alto orden existente, implementado dentro de la biblioteca de elementos finitos deal.II, para manejar las ecuaciones de Euler compresibles con efectos rotacionales.

• Ecuaciones Gobernantes: El solver aborda las ecuaciones de Euler compresibles invíscidas en un dominio acotado, incorporando la aceleración de Coriolis ($-f \mathbf{k} \times \mathbf{u}$) y la aceleración gravitatoria. El sistema se cierra utilizando una ecuación de estado de gas ideal calóricamente perfecto.
• Discretización Temporal: Se emplea un esquema de Runge-Kutta Aditivo IMEX (ARK) de segundo orden (específicamente el método TR-BDF2). Las contribuciones rígidas (gradientes de presión, gravedad y rotación) se tratan implícitamente, mientras que los términos advectivos no rígidos se tratan explícitamente. Esto permite que el paso de tiempo sea determinado por escalas de flujo meteorológicamente relevantes en lugar de la velocidad del sonido.
• Discretización Espacial: Se utiliza un método DG de alto orden en mallas adaptativas no conformes (cuadriláteras en 2D, hexaédricas en 3D). La formulación utiliza una base nodal en elementos de producto tensorial para habilitar una evaluación eficiente de operadores libres de matrices.
• Tratamiento del Término Rotacional: Se introducen y comparan dos estrategias alternativas para manejar el término rotacional dentro de la iteración de punto fijo no lineal (Picard):
  1. Estrategia R1: Los términos rotacionales y gravitatorios se retrasan al lado derecho, evaluados utilizando la iteración de Picard anterior.
  2. Estrategia R2: Los términos rotacionales y gravitatorios se tratan implícitamente en el lado izquierdo del sistema acoplado de momento-energía.
• Implementación Libre de Matrices: Una contribución técnica clave es la derivación de una formulación para la inversa del operador $(A + R)$, donde $A$ es la matriz de masa y $R$ representa el término de Coriolis. Al explotar la estructura algebraica del operador de Coriolis (específicamente que $R = \beta A J$, donde $J$ es una rotación planar), los autores derivan una factorización $(A+R)^{-1} = (I + \beta J)^{-1}A^{-1}$. Esto permite que la aplicación de la inversa se calcule mediante operaciones estándar libres de matrices seguidas de una transformación algebraica local, evitando la necesidad de resolver sistemas lineales con el operador acoplado completo.

Contribuciones Clave

1. Extensión a Flujos Rotatorios: La integración exitosa de la fuerza de Coriolis en un marco IMEX-DG completamente compresible y de alto orden dentro de la biblioteca deal.II.
2. Solver Rotacional Eficiente Libre de Matrices: El desarrollo de una factorización algebraica específica para el operador de Coriolis que preserva la eficiencia de las implementaciones libres de matrices, evitando el costo computacional de resolver grandes sistemas lineales acoplados para el término rotacional.
3. Comparación de Estrategias Iterativas: Una comparación sistemática de dos tratamientos de punto fijo no lineal (R1 y R2) para términos rotacionales, evaluándolos en términos de precisión, robustez y eficiencia computacional.
4. Validación en Mallas Adaptativas: La primera aplicación de tal marco IMEX-DG rotatorio en mallas adaptativas no conformes, demostrando la capacidad de realizar refinamiento local cerca de la orografía mientras se mantiene la estabilidad y la precisión.

Resultados Numéricos
La formulación propuesta fue validada mediante una serie de pruebas de referencia en dos y tres dimensiones:

• Ondas de Inercia-Gravedad 2D: Estudios de convergencia en una prueba de referencia estándar de ondas de inercia-gravedad mostraron que ambas estrategias R1 y R2 alcanzan tasas de convergencia cercanas a las óptimas, consistentes con el orden teórico del esquema. Los resultados confirmaron que el solver captura con precisión la propagación de ondas y la preservación de la dinámica equilibrada.
• Ondas a Escala Planetaria: Una prueba de dominio grande (extensión horizontal de 60,000 km) demostró la robustez del solver durante tiempos de integración largos (96 horas). El error de equilibrio geostrofico discreto permaneció acotado y por debajo de $10^{-7}$, indicando que la falta de una discretización exactamente bien equilibrada no condujo a oscilaciones espurias ni inestabilidad.
• Flujo 3D sobre Orografía:
  • Colina de Pendiente Media: Las simulaciones en un dominio de 60x40x16 km mostraron que la inclusión de la rotación rompe correctamente la simetría axial de la estela y los patrones de ondas, consistente con las expectativas físicas. Los resultados fueron consistentes entre mallas uniformes y no conformes, siendo estas últimas capaces de proporcionar una mejor resolución de gradientes cerca de la montaña con un costo computacional reducido.
  • Colina a Gran Escala: Una nueva prueba a mayor escala (240x240x30 km) destacó la asimetría pronunciada inducida por la rotación en los campos de velocidad, la cual fue capturada claramente sin oscilaciones espurias.
• Escalabilidad Paralela: Las pruebas de escalado fuerte en hasta 2048 núcleos (utilizando el superordenador MeluXina) demostraron que la inclusión de términos rotacionales no degrada significativamente la eficiencia paralela del solver original no rotatorio. El solver mantuvo un comportamiento de escalado cercano al ideal, y la configuración de malla no conforme mostró una escalabilidad comparable al caso de malla uniforme.

Significado y Afirmaciones
El artículo afirma que el marco IMEX-DG rotatorio propuesto extiende con éxito las capacidades de los solvers atmosféricos compresibles existentes para incluir la rotación planetaria sin comprometer la precisión, la robustez o la escalabilidad del método original. Los autores enfatizan que el tratamiento libre de matrices del operador rotacional es un habilitador crítico para la eficiencia, permitiendo que el solver maneje la complejidad adicional de la fuerza de Coriolis con una sobrecarga mínima.

Los resultados indican que el marco es adecuado para simular flujos atmosféricos a gran escala donde la rotación es un efecto físico de orden principal. La comparación entre los dos tratamientos rotacionales sugiere que son numéricamente equivalentes para los casos probados, proporcionando flexibilidad en la implementación. Los autores concluyen que el método captura correctamente la asimetría y las estructuras de ondas inducidas por la rotación, validando su uso tanto para pruebas de referencia estándar como para configuraciones atmosféricas más realistas y a gran escala. Se identifica que el trabajo futuro se centrará en configuraciones 3D a escala planetaria, comparaciones detalladas con discretizaciones exactamente bien equilibradas y extensiones a geometrías esféricas.