Local cubic spline interpolation for Vlasov-type equations on a multi-patch geometry

Resumen Técnico: Interpolación de Spline Cúbico Local para Ecuaciones de tipo Vlasov en una Geometría de Múltiples Parches (Multi-patch)

Planteamiento del Problema
El artículo aborda la resolución numérica de ecuaciones de tipo Vlasov, específicamente modelos de centro de guía (guiding-center), en geometrías complejas de múltiples parches típicas de las simulaciones de tokamaks. En estas simulaciones, las líneas de malla se alinean con las superficies de flujo magnético para minimizar las variaciones a lo lo largo de las líneas de campo. Esta geometría crea desafíos específicos: el núcleo presenta un punto O donde las líneas de malla convergen, y el borde presenta un punto X (separatriz) donde las líneas divergen. Si bien la división del dominio global en subdominios (parches) permite el refinamiento local y la optimización de mallas uniformes dentro de cada parche, introduce dificultades para mantener la regularidad de la solución en las interfaces de los parches. Estudios previos que utilizaban continuidad $C^0$ (emparejando solo los valores de la función) en las interfaces fueron identificados como potencialmente inestables para esquemas de advección semi-Lagrangianos. El problema central es cómo imponer una continuidad de orden superior ($C^1$) y calcular con precisión las derivadas de la interfaz para interpolaciones de B-spline cúbicas locales en mallas de múltiples parches arbitrarias (uniformes y no uniformes) sin sacrificar la eficiencia de los cálculos locales.

Metodología
Los autores proponen un método semi-Lagrangiano hacia atrás (BSL, backward semi-Lagrangian) utilizando B-splines cúbicas locales definidas en cada parche. Para asegurar la estabilidad y la precisión, el método impone regularidad $C^1$ en las interfaces de los parches mediante el emparejamiento tanto de los valores de la función como de sus primeras derivadas.

1. Reconstrucción de la Derivada: La innovación central es una fórmula generalizada para calcular las derivadas de la interfaz. Basándose en una fórmula de Crouseilles et al. [11] diseñada para mallas uniformes, los autores extienden el enfoque a mallas no uniformes arbitrarias utilizando polinomios de Hermite.

   • Formulación Exacta: Al imponer continuidad $C^2$ entre polinomios de Hermite locales definidos en celdas adyacentes, los autores derivan una relación recursiva que vincula la derivada en una interfaz con las derivadas en las interfaces vecinas y los valores de la función. Esto resulta en un sistema lineal global pequeño (de tamaño $(N_p-1) \times (N_p-1)$ para $N_p$ parches) que debe resolverse para obtener derivadas equivalentes a las de un spline global.
   • Formulación Aproximada: Los autores analizan el comportamiento asintótico de los coeficientes en el sistema global. Demuestran que la influencia de las interfaces distantes decae exponencialmente a medida que aumenta el número de celdas por parche. Por consiguiente, para parches suficientemente grandes (por ejemplo, $\ge 30$ celdas), la dependencia de las derivadas distantes se vuelve insignificante (alcanzando la precisión de máquina). Esto permite una aproximación local truncada, donde las derivadas se calculan utilizando solo un número limitado de celdas vecinas, evitando la resolución global pero manteniendo una alta precisión.
   • Fórmulas Explícitas: Para el caso específico de parches uniformes, los autores proporcionan expresiones cerradas explícitas para los coeficientes, optimizando aún más el cálculo.

2. Extensión 2D: El método se extiende a geometrías de múltiples parches en 2D, manejando tres casos:

   • Mallas Conformes: Una aplicación directa de producto tensorial de las fórmulas 1D para calcular las derivadas en $r$, las derivadas en $\theta$ y las derivadas cruzadas.
   • Mallas No Conformes: Donde los puntos de ruptura no se alinean a través de los parches. Los autores proponen un algoritmo que calcula progresivamente las derivadas en interfaces más gruesas e interpola estas sobre los parches refinados para imponer la continuidad $C^1$.
   • T-Joints (Uniones en T): Conexiones que involucran tres o más parches donde un borde está solo parcialmente conectado. El método adapta el orden de computación de la derivada basándose en la geometría para asegurar que todas las derivadas necesarias estén disponibles.

3. Implementación: El algoritmo se implementa en la biblioteca Gyselalib++. La ecuación de Poisson se resuelve mediante un Método de Elementos Finitos (FEM) con B-splines cúbicas sobre un spline global equivalente, mientras que la ecuación de advección se resuelve mediante el método BSL propuesto sobre splines locales.

Contribuciones Clave

• Generalización de las Fórmulas de Derivada de Interfaz: El artículo generaliza la fórmula de cálculo de derivadas de mallas uniformes a mallas no uniformes arbitrarias utilizando polinomios de Hermite, permitiendo el manejo de situaciones no conformes y el refinamiento local.
• Análisis de Estabilidad: Los autores demuestran que la regularidad $C^0$ es insuficiente para la estabilidad en esquemas BSL en dominios de múltiples parches y que se requiere regularidad $C^1$.
• Compromiso entre Local y Global: El trabajo establece que, si bien un sistema global exacto garantiza la continuidad $C^2$, una aproximación local utilizando un stencil truncado (por ejemplo, 30 celdas) logra una precisión de máquina equivalente a la solución global, ofreciendo una alternativa computacionalmente eficiente.
• Manejo de Geometrías Complejas: El método maneja con éxito las uniones en T (T-joints) y las mallas no conformes, las cuales son necesarias para modelar secciones transversales de tokamaks complejos que involucran tanto puntos O como puntos X.

Resultados Numéricos
El método fue validado utilizando la biblioteca Gyselalib++ con las siguientes pruebas:

1. Interpolación de Spline: Las pruebas en parches concéntricos (uniformes y no uniformes) confirmaron que las representaciones de spline locales convergen al spline global equivalente con orden 4. El error en las derivadas de la interfaz usando la aproximación truncada coincidió con la fórmula exacta una vez que el número de celdas excedió $\approx 30$, alcanzando la precisión de máquina ($10^{-13}$).
2. Pruebas de Advección:
   • Parches no conformes: Los splines locales en una malla no conforme lograron una precisión comparable a la de una malla global refinada, a pesar de utilizar menos celdas totales.
   • T-joints: En una configuración de T-joint en 2D, los resultados de los splines locales concordaron con los resultados del spline global hasta la precisión de máquina ($\approx 6 \times 10^{-13}$).
3. Modelo de Centro de Guía (Inestabilidad Diocotrona):
   • Las simulaciones en mallas de múltiples parches conformes reprodujeron los resultados de las simulaciones de spline global (tasas de crecimiento y conservación de masa) con errores consistentes con la precisión de la interpolación.
   • Una simulación con una malla no óptima (parche central despoblado e interfaces desplazadas) mostró que, si bien la elección de la malla afectó la magnitud del error absoluto, el método de spline local aún rastreó correctamente el comportamiento de la solución global, validando la robustez del enfoque.

Significado y Reivindicaciones
El artículo afirma proporcionar un marco robusto para resolver ecuaciones de tipo Vlasov en geometrías de múltiples parches, las cuales son esenciales para simulaciones realistas de tokamaks. La importancia principal radica en permitir el uso de interpolación de spline cúbico local con condiciones de contorno de Hermite para mantener la estabilidad y la alta precisión a través de las interfaces de los parches.

Los autores enfatizan que su método elimina la restricción encontrada en estudios previos respecto al número de puntos por parche en casos conformes, reemplazándola con el requisito de resolver un sistema global pequeño (o una aproximación local). Señalan que el método está actualmente validado para geometrías con puntos O. Si bien la introducción de las uniones en T es un paso hacia la conexión de geometrías de punto O y punto X, los autores declaran explícitamente que el manejo de las singularidades del punto X requiere un estudio adicional y aún no se ha resuelto completamente en este trabajo. El método se presenta como un paso significativo hacia solvers semi-lagrangianos de alto orden y eficientes para geometrías de confinamiento magnético complejas.