Tutorial discontinuo del método Galerkin: derivación, flujos e implementación de Python

El método de Galerkin discontinuo, generalmente acortado a DG, es una técnica numérica para resolver ecuaciones diferenciales parciales. Combina varias propiedades útiles de los métodos de volumen finito y de elementos finitos. Al igual que el método de volumen finito, DG puede conservar las cantidades localmente a través de flujos cuidadosamente definidos. Al igual que el método de elementos finitos, admite mallas flexibles y aproximaciones polinómicas dentro de cada elemento.

La característica definitoria de DG es que la solución numérica no necesita permanecer continua a través de los límites de los elementos. Cada elemento tiene su propia representación polinómica. Los elementos vecinos se comunican a través de flujos numéricos evaluados en sus interfaces compartidas.

Esta estructura hace que DG sea útil para la propagación de ondas, ecuaciones de transporte, flujo compresible, sistemas de aguas poco profundas, electromagnéticos y problemas que contienen choques o interfaces de materiales. También es atractivo para la computación en paralelo porque muchos cálculos se pueden realizar de forma independiente dentro de cada elemento.

¿En qué se diferencia DG de otros métodos?

Un método de volumen finito generalmente almacena un valor promedio en cada celda y calcula el flujo de una cantidad conservada a través de las caras de las celdas. Proporciona una fuerte conservación local, pero muchos esquemas básicos de volumen finito utilizan aproximaciones espaciales de orden relativamente bajo.

Un método de elementos finitos continuos representa la solución con funciones de base polinómica que se conectan continuamente a través de los límites de los elementos. Esto ofrece flexibilidad geométrica y aproximación de alto orden, pero las discontinuidades no se pueden representar directamente sin tratamiento adicional.

DG utiliza funciones de base polinomial dentro de cada elemento mientras permite valores separados en ambos lados de una interfaz. Los flujos numéricos determinan cómo interactúan esos valores. Por lo tanto, el método proporciona conservación local, aproximación de alto orden y soporte directo para soluciones discontinuas.

La ecuación de advección continua

La ecuación de advección lineal unidimensional proporciona una clara introducción al DG:

∂u/∂t + ∂f(u)/∂x = 0

Para la velocidad de transporte constante a, el flujo físico es:

f(u) = a u

La ecuación entonces se convierte en:

∂u/∂t + a ∂u/∂x = 0

Este modelo describe un perfil que se mueve a través del dominio sin cambiar de forma. Cuando a es positivo, la información viaja de izquierda a derecha. Cuando es negativo, la información viaja de derecha a izquierda.

Dividir el dominio en elementos

Suponga que el dominio se extiende de cero a L. DG lo divide en elementos no superpuestos:

Ω = K₁ ∪ K₂ ∪ ... ∪ Kₙ

Dentro de cada elemento, la solución numérica está representada por un polinomio:

uₕ(x, t) = Σ Uᵢ(t) φᵢ(x)

Las funciones φᵢ son funciones de base local, y los coeficientes Uᵢ son los grados de libertad que cambian con el tiempo.

Un polinomio de grado cero almacena un valor constante en cada elemento. Una aproximación lineal utiliza dos grados de libertad locales en una dimensión. Una aproximación cuadrática utiliza tres.

A diferencia de los elementos finitos continuos, DG no fuerza al polinomio de un elemento para que coincida con el polinomio en el siguiente. En una interfaz compartida, la solución puede tener, por lo tanto, un valor izquierdo y un valor derecho.

Derivando la forma débil en cuanto a elementos

La formulación de DG comienza multiplicando la ecuación de gobierno por una función de prueba v e integrando sobre un elemento:

∫K (∂uₕ/∂t) v dx + ∫K (∂f(uₕ)/∂x) v dx = 0

La integración por partes mueve la derivada espacial del flujo físico a la función de prueba:

∫K (∂uₕ/∂t) v dx
- ∫K f(uₕ) ∂v/∂x dx
+ f(uₕ) v |∂K = 0

El término límite es esencial. Describe el flujo a través de los lados izquierdo y derecho del elemento.

Debido a que la solución DG puede ser discontinua en una interfaz, el flujo físico no está definido de manera única allí. Un elemento proporciona un valor desde la izquierda, mientras que su vecino proporciona otro valor desde la derecha. DG reemplaza el flujo físico ambiguo con un flujo numérico:

f̂(u⁻, u⁺)

La forma débil final para un elemento es:

∫K (∂uₕ/∂t) v dx
- ∫K f(uₕ) ∂v/∂x dx
+ f̂R vR
- f̂L vL = 0

Esta ecuación contiene una contribución de volumen local y dos contribuciones de interfaz. Los flujos numéricos son los únicos términos que conectan directamente los elementos vecinos.

Por qué importan los flujos numéricos

Un flujo numérico debe ser consistente. Cuando los dos valores de la interfaz son iguales, debería reproducir el flujo físico:

f̂(u, u) = f(u)

También debe proporcionar una estabilidad adecuada. Para los problemas de transporte, el flujo debe respetar la dirección en la que se mueve la información.

La elección de flujo controla la fuerza con que interactúan los elementos vecinos y cuánta disipación numérica se introduce. Un flujo inadecuado puede crear oscilaciones, suavizado excesivo o velocidades de onda incorrectas.

Flujo contra el viento

Para la advección lineal con velocidad positiva, la información proviene del lado izquierdo de una interfaz. Por lo tanto, el flujo contra el viento es:

f̂(u⁻, u⁺) = a u⁻    when a > 0

Para la velocidad negativa, la información proviene de la derecha:

f̂(u⁻, u⁺) = a u⁺    when a < 0

Una expresión compacta es:

f̂ = a⁺u⁻ + a⁻u⁺

a⁺ = max(a, 0)
a⁻ = min(a, 0)

El flujo contra el viento es simple, estable y ampliamente utilizado. Introduce cierta disipación numérica, pero la cantidad generalmente disminuye a medida que se refina la malla o aumenta el orden polinomial.

Flujo local de Lax-Friedrichs

El flujo local de Lax-Friedrichs, también llamado Flux Rusanov, se usa con frecuencia para las leyes de conservación no lineales:

f̂(u⁻, u⁺)
= 0.5[f(u⁻) + f(u⁺)]
- 0.5 α(u⁺ - u⁻)

El parámetro α es una estimación de la mayor velocidad característica en la interfaz. Para la advección lineal constante, suele ser |a|.

La primera parte promedia los flujos físicos. La segunda parte agrega disipación que ayuda a controlar los modos de interfaz inestables. Este flujo es robusto y fácil de implementar, aunque puede ser más difusivo que un solucionador de Riemann aproximado especializado.

Solucionadores de Riemann

Los sistemas hiperbólicos no lineales pueden contener varias ondas que viajan a diferentes velocidades. Los ejemplos incluyen las ecuaciones de Euler para el flujo compresible y las ecuaciones de aguas poco profundas.

Un solucionador de Riemann examina los estados izquierdo y derecho en una interfaz y estima las ondas producidas por su interacción. Los solucionadores exactos de Riemann pueden ser costosos, por lo que los códigos DG prácticos comúnmente utilizan métodos aproximados como ROE, HLL, HLLC o flujos de Rusanov.

La elección correcta depende de la ecuación, la precisión deseada, los requisitos de robustez y la capacidad de preservar importantes propiedades físicas.

El sistema de matriz semidiscreto

Después de seleccionar las funciones de base y prueba, la forma débil se puede expresar como un sistema de ecuaciones diferenciales ordinarias:

M dU/dt = R(U)

La matriz de masa contiene integrales de productos de función base:

Mᵢⱼ = ∫K φᵢ φⱼ dx

El residuo R(U) incluye las contribuciones de flujo derivado y numérico.

Debido a que las funciones de base de DG pertenecen a elementos individuales, la matriz de masa global tiene una estructura de bloque-diagonal. Cada elemento aporta un pequeño bloque independiente. Estos bloques se pueden invertir por separado.

La matriz de masas no es automáticamente diagonal para cada base y regla de integración. Se vuelve diagonal o aproximadamente diagonal en formulaciones DG nodales comunes que utilizan puntos de interpolación y cuadratura coincidentes, como la colocación de Gauss-Lobatto. Otras formulaciones utilizan matrices de elementos densos pequeños.

Elegir funciones básicas y cuadratura

Los métodos DG comúnmente usan bases polinómicas nodales o nodales. Una base modal representa la solución a través de modos polinómicos, a menudo basado en polinomios de Legendre. Una base nodal almacena valores de solución en los nodos de interpolación dentro del elemento.

Las formulaciones nodales son convenientes porque los valores de la interfaz se pueden obtener directamente cuando los nodos se colocan en los límites de los elementos. Los puntos Gauss-Lobatto incluyen ambos puntos finales del elemento de referencia, mientras que los puntos gauss permanecen dentro de él.

La cuadratura numérica evalúa las integrales en la forma débil. La regla de cuadratura debe ser lo suficientemente precisa para el grado polinómico y cualquier término no lineal. La cuadratura insuficiente puede introducir errores de alias e inestabilidad.

Orden y precisión polinómicas

Para soluciones suficientemente suaves, un método DG bien diseñado puede lograr un error proporcional a aproximadamente h^(p+1). Por lo tanto, aumentar el grado polinomial puede mejorar la precisión sin agregar más elementos.

El orden superior no siempre es mejor. Los choques y las discontinuidades agudas pueden crear oscilaciones cerca del salto. Es posible que se necesiten limitadores, viscosidad artificial, filtrado o técnicas de captura de golpes.

Un solucionador de DG lineal en numpy

El siguiente ejemplo educativo implementa un método DG lineal p = 1 para la advección unidimensional periódica. Cada elemento contiene dos grados de libertad ubicados en sus puntos finales.

La implementación utiliza el elemento forma débil, un flujo contra el viento y un método Runge-Kutta que preserva la estabilidad fuerte de tercer orden.

import numpy as np

# Domain and model parameters
length = 1.0
number_of_elements = 80
velocity = 1.0
final_time = 0.5

element_width = length / number_of_elements
jacobian = element_width / 2.0

# Reference-element mass matrix for linear basis functions
mass_reference = np.array([
    [2.0 / 3.0, 1.0 / 3.0],
    [1.0 / 3.0, 2.0 / 3.0]
])

mass_matrix = jacobian * mass_reference
inverse_mass = np.linalg.inv(mass_matrix)

# S[i, j] = integral(phi_j * derivative(phi_i)) on [-1, 1]
volume_matrix = np.array([
    [-0.5, -0.5],
    [ 0.5,  0.5]
])

left_vector = np.array([1.0, 0.0])
right_vector = np.array([0.0, 1.0])

# Physical coordinates of local DG nodes
left_edges = np.arange(number_of_elements) * element_width
right_edges = left_edges + element_width
coordinates = np.column_stack((left_edges, right_edges))

# Smooth periodic initial condition
solution = (
    0.5
    + 0.5 * np.sin(2.0 * np.pi * coordinates / length)
)

def upwind_flux(left_state, right_state, speed):
    positive_speed = max(speed, 0.0)
    negative_speed = min(speed, 0.0)

    return (
        positive_speed * left_state
        + negative_speed * right_state
    )

def spatial_residual(values):
    residual = np.zeros_like(values)

    for element in range(number_of_elements):
        left_neighbor = (element - 1) % number_of_elements
        right_neighbor = (element + 1) % number_of_elements

        # States on the left interface
        left_inside = values[element, 0]
        left_outside = values[left_neighbor, 1]

        # States on the right interface
        right_inside = values[element, 1]
        right_outside = values[right_neighbor, 0]

        flux_left = upwind_flux(
            left_outside,
            left_inside,
            velocity
        )

        flux_right = upwind_flux(
            right_inside,
            right_outside,
            velocity
        )

        local_rhs = (
            velocity * volume_matrix @ values[element]
            + flux_left * left_vector
            - flux_right * right_vector
        )

        residual[element] = inverse_mass @ local_rhs

    return residual

# Conservative time-step estimate
time_step = 0.1 * element_width / abs(velocity)
current_time = 0.0

while current_time < final_time:
    dt = min(time_step, final_time - current_time)

    # SSP-RK3 stage 1
    stage_one = solution + dt * spatial_residual(solution)

    # SSP-RK3 stage 2
    stage_two = (
        0.75 * solution
        + 0.25 * (
            stage_one
            + dt * spatial_residual(stage_one)
        )
    )

    # SSP-RK3 stage 3
    solution = (
        (1.0 / 3.0) * solution
        + (2.0 / 3.0) * (
            stage_two
            + dt * spatial_residual(stage_two)
        )
    )

    current_time += dt

print("Simulation completed")
print("Final time:", current_time)
print("Minimum value:", solution.min())
print("Maximum value:", solution.max())

Este ejemplo se limita intencionalmente a un problema lineal suave. Los solucionadores de DG de producción también necesitan tratamientos de contorno robustos, mapeos multidimensionales, cuadratura precisa, evaluación de flujo no lineal, limitadores e integración de tiempo más avanzada.

Integración del tiempo y la condición CFL

La discretización espacial DG crea un sistema de ecuaciones diferenciales ordinarias. Un método explícito como Runge-Kutta puede entonces avanzar en ese sistema a tiempo.

El paso de tiempo máximo estable depende de la velocidad de la onda, el tamaño del elemento, el grado polinomial, el flujo y el integrador de tiempo. Un escalado común es:

Δt ∝ h / [(2p + 1)|a|]

La constante de estabilidad exacta depende del método. Debe establecerse a través del análisis, la documentación o las pruebas numéricas en lugar de tratarse como universales.

El aumento del orden polinomial normalmente reduce el paso de tiempo explícito estable más grande. Por lo tanto, la DG de alto orden puede requerir más pasos incluso cuando necesita menos elementos.

DG para problemas dominados por la advección

Los métodos estándar continuos de Galerkin pueden desarrollar oscilaciones cuando la advección domina la difusión. Los métodos continuos estabilizados como SUPG modifican las funciones de prueba para agregar control a lo largo de las líneas de corriente.

DG maneja el transporte a través de flujos de interfaz. Los flujos de Riemann o aproximados de Riemann introducen la estabilización consciente de la dirección y preservan la conservación local.

SUPG sigue siendo efectivo para muchos problemas suaves dominados por convección. La DG es atractiva cuando la conservación local, las soluciones discontinuas, las mallas complejas o la adaptabilidad en cuanto a elementos son requisitos centrales.

DG frente a FVM y FEM continua

Ventajas del método DG

• La conservación local está integrada en el equilibrio del flujo de la interfaz.
• El grado polinomial puede variar entre elementos.
• Las discontinuidades no violan el espacio de aproximación.
• Los cálculos de elementos se adaptan bien al hardware paralelo.
• Se pueden soportar mallas complejas y no estructuradas.
• Se puede combinar el refinamiento de malla y el enriquecimiento polinómico.
• Se pueden seleccionar diferentes flujos numéricos para diferentes ecuaciones.

Limitaciones de DG

DG normalmente utiliza más grados de libertad que un método de elementos finitos continuos del mismo grado polinomial porque los elementos vecinos no comparten valores de interfaz.

El método también requiere un cuidadoso diseño de flujo. Las ecuaciones hiperbólicas, elípticas y mixtas necesitan diferentes tratamientos de interfaz. Los operadores de difusión requieren formulaciones como penalización interior, DG local o técnicas relacionadas.

Los métodos DG explícitos de orden superior pueden tener límites de paso de tiempo restrictivos. Los problemas no lineales pueden requerir limitadores, flujos estables de entropía, preservación de la positividad o viscosidad artificial.

Estas características hacen que DG sea potente pero más difícil de implementar correctamente que la diferencia finita básica o los esquemas de volumen finito.

Errores de implementación comunes

Un error común es usar el flujo físico directamente en una interfaz discontinua sin definir cómo se deben combinar los dos valores vecinos.

Otros errores frecuentes incluyen:

• Usar el signo erróneo para la contribución del límite izquierdo o derecho
• Aplicación de orientaciones vectoriales normales inconsistentes
• Asumiendo que la matriz de masa es siempre diagonal
• Usando una cuadratura que es demasiado débil para términos no lineales
• Ignorar la dependencia polinomial-grado del límite CFL
• Aplicar incorrectamente las condiciones de contorno periódicas o físicas
• Uso de polinomios de alto orden cerca de amortiguadores sin limitador
• Probando solo la salida visual en lugar de la convergencia y la conservación

Cómo validar un solucionador de DG

Comience con un problema que tiene una solución analítica conocida. La advección lineal periódica es útil porque el perfil exacto simplemente se desplaza a t.

Ejecute el solucionador con varias resoluciones de malla y mida una norma de error. Para una solución suave, la tasa de convergencia observada debe abordar el orden teórico del método.

Compruebe la conservación integrando la solución en todo el dominio. Para la advección periódica, la masa total debe permanecer casi constante.

Pruebe las velocidades positivas y negativas para verificar la dirección contra el viento. Los datos iniciales constantes deben permanecer constantes. Los estados de contorno también deben probarse por separado antes de pasar a ecuaciones no lineales.

Cuándo elegir DG

DG es una fuerte opción cuando el problema contiene ondas, choques, interfaces de materiales o comportamiento dominado por el transporte. También es útil cuando la conservación local es esencial o cuando la simulación se beneficia de la precisión de alto orden en una malla no estructurada.

Un método básico de volumen finito puede seguir siendo más simple para problemas de conservación de bajo orden. Los elementos finitos continuos pueden ser más eficientes para problemas elípticos o estructurales suaves donde no se esperan discontinuidades.

El método numérico debe seguir la estructura matemática de la PDE en lugar de la popularidad actual o la disponibilidad de software.

Conclusión

El método Galerkin discontinuo representa la solución con polinomios independientes dentro de cada elemento. Su formulación débil produce términos de contorno naturales que son reemplazados por flujos numéricos. Estos flujos controlan cómo se mueve la información entre los elementos y permiten el método completo para conservar las cantidades localmente.

DG combina una aproximación polinomial de alto orden, flexibilidad geométrica, manejo de discontinuidad y una fuerte localidad de elementos. Estas ventajas lo hacen valioso para la propagación de olas, las leyes de conservación hiperbólica, los sistemas dominados por convección y las grandes simulaciones paralelas.

El método también introduce complejidad adicional. Los flujos, cuadratura, restricciones de paso de tiempo, limitadores y condiciones de contorno deben seleccionarse cuidadosamente. Se debe desarrollar gradualmente una implementación confiable, comenzando con una ecuación lineal simple y verificada mediante pruebas de convergencia y conservación.

Para los investigadores que ya entienden los métodos de volumen finito o elementos finitos, DG proporciona un siguiente paso natural hacia los solucionadores de PDE de alto orden localmente conservadores en Python científico.