El número de Mach tiende a cero y el paso temporal no se inmuta — Reproduciendo un esquema IMEX TVD

El esquema AP de primer orden corrió sin problemas. Bajara el número de Mach a $10^{-2}$ o a $10^{-4}$, el paso temporal no se inmutaba. En cuanto le monté encima la discretización temporal de segundo orden ARS(2,2,2), brotaron sobreimpulsos a ambos lados del pulso. No era una divergencia. La amplitud quedaba acotada, pero por más pasos que corriera no desaparecía, y solo se esfumó al bajar el paso temporal hasta el CFL acústico (Courant–Friedrichs–Lewy: la condición sobre el paso temporal que impide que la información avance más de una celda por paso). Es justo la restricción que el esquema pretendía evitar.

Al bajar el deslizador de ε hasta $10^{-4}$ las características acústicas se tumban, y los botones Explicit e IMEX intercambian el número de pasos necesario para el mismo $t=1$.

En esa marca la lectura es 4489 pasos explícitos frente a 89 de IMEX. Conviene mirar si los rayos acústicos se pegan a la horizontal y en qué botón se encoge la diferencia.

El artículo que leí hoy dice que ese fracaso no es un error de implementación sino un teorema. De este post quedan tres cosas. Por qué el segundo orden y un paso temporal sin restricción no pueden sostenerse a la vez, cómo esquivaron los autores ese muro, y qué dicen los números reales de una reproducción en 30 líneas de Python. Adelantando el final: las oscilaciones desaparecieron exactamente como indica el artículo, y a cambio hubo un precio mayor de lo previsto.

• Título: Second order Implicit-Explicit Total Variation Diminishing schemes for the Euler system in the low Mach regime
• Autores: Giacomo Dimarco (Ferrara), Raphaël Loubère (Bordeaux), Victor Michel-Dansac, Marie-Hélène Vignal (Toulouse)
• Fuente: Preprint enviado a Elsevier, 2017-10-23
• DOI: El resumen interno no tiene campo DOI — el original hay que buscarlo por su título.
• Resumen en una línea: Combinar de forma convexa un esquema AP de primer orden con un esquema IMEX de segundo orden, para usar un paso temporal independiente del número de Mach sin perder la propiedad TVD (variación total no creciente: la variación total de la solución nunca aumenta con el tiempo) ni la propiedad $L^\infty$.

Q1. ¿Por qué muere el paso temporal cuando el número de Mach se hace pequeño?#

El punto de partida es el sistema de Euler isentrópico, adimensionalizado tomando el cuadrado del número de Mach como $\varepsilon$.

$\rho$ es la densidad, $U$ es el campo de velocidad y la primera ecuación es la conservación de masa.

$\rho U$ es la cantidad de movimiento, $p(\rho)=\rho^\gamma$ es la presión, $\gamma$ es la razón de calores específicos y $1/\varepsilon$ es el inverso del cuadrado del número de Mach.

Conviene fijarse en que solo el término de presión lleva $1/\varepsilon$. Cuando el gradiente de presión crece, la velocidad del sonido crece con él. Más exactamente, la velocidad del sonido escala como $1/\sqrt{\varepsilon}$. Un solver totalmente explícito debe seguir a la onda más rápida, así que queda atado a $\Delta t \le \Delta x / (|u| + c/\sqrt{\varepsilon})$. Con $\varepsilon = 10^{-4}$ eso es un paso unas 100 veces menor que el que pediría la escala temporal convectiva.

El problema es que nadie necesita ese factor 100. En flujo de bajo Mach lo que interesa son las estructuras convectivas, no las ondas acústicas. De hecho, en el límite $\varepsilon \to 0$ el sistema converge al Euler incompresible. La densidad queda fija en una constante $\rho_0$, se cumple $\nabla\cdot U = 0$ y la ecuación de cantidad de movimiento toma la forma $\rho_0\partial_t U + \rho_0\nabla\cdot(U\otimes U) + \nabla\pi_1 = 0$. Aquí $\pi_1$ es la perturbación de presión de primer orden y hace de multiplicador de Lagrange que sostiene la restricción de incompresibilidad.

La propiedad que hace falta entonces es AP (asymptotic-preserving, preservación asintótica). La definición es simple. Es la propiedad de que el esquema degenere, cuando $\varepsilon \to 0$, en una discretización consistente de las ecuaciones límite. Es decir, que el límite se capture bien sin encoger la malla ni el paso temporal para ajustarlos a $\varepsilon$.

Los 4489 y 89 ya leídos en el deslizador con $\varepsilon=10^{-4}$ son esta escala $1/\sqrt{\varepsilon}$ escrita como número de pasos. 4489/89 da cerca de 50, por debajo de $1/\sqrt{\varepsilon}=100$, porque el denominador explícito todavía arrastra una velocidad convectiva.

Q2. ¿Qué hay que pasar a implícito para que el CFL se suelte?#

Hacerlo todo implícito significa resolver un sistema no lineal en cada paso, y ese costo es inasumible. El artículo parte en dos piezas el flujo de las variables conservativas $W = (\rho, \rho U)$.

$F_e(W) = (0,\ \rho U \otimes U)$ es el flujo convectivo tratado de forma explícita, y $F_i(W) = (\rho U,\ p(\rho)/\varepsilon \cdot I)$ es el flujo acústico tratado de forma implícita.

Lo esencial es que esta partición no es arbitraria. Poner implícito el gradiente de presión da consistencia asintótica. Poner implícito el flujo de masa da estabilidad uniforme, independiente de $\varepsilon$. Si solo se pasa uno de los dos, se rompe una de las dos propiedades.

Queda entonces cómo resolver el sistema implícito. El truco es tomar la divergencia. Al sustituir la divergencia de la ecuación implícita de cantidad de movimiento en la ecuación de masa, la cantidad de movimiento se cancela y queda una única ecuación elíptica no lineal para $\rho^{n+1}$.

El último término del lado izquierdo, $\Delta p(\rho)^{n+1}$, es el laplaciano de presión implícito, y los términos que lo preceden son todos fuentes explícitas evaluadas en el instante $n$.

Una vez resuelta esta ecuación elíptica para $\rho^{n+1}$, la cantidad de movimiento se actualiza con una sola sustitución explícita. Su estructura es la misma que la del paso de corrección de presión de un solver basado en presión. La única parte explícita que queda es $F_e$, así que la restricción sobre el paso temporal solo mira la velocidad convectiva. El resultado es $\Delta t \le \Delta x / \max_j(2|u_j^n|)$. El factor 2 aparece porque el mayor autovalor del flujo explícito $\rho U \otimes U$ es $2u$. No hay $\varepsilon$ por ninguna parte.

Separar la acústica del paso temporal y dejar solo la convección es el mismo movimiento que la división de Lagrange–Projection. Aquí el corte se cierra con una elíptica para la densidad, y el tope del paso queda en $\max_j(2|u_j^n|)$. $\varepsilon$ no entra ni una vez en ese tope.

Q3. ¿Por qué aparecen oscilaciones al subir a segundo orden?#

El segundo orden en tiempo llega con un Runge–Kutta IMEX (implicit-explicit: una integración temporal que mezcla tratamiento implícito y explícito término a término) ARS(2,2,2). El coeficiente $\beta = 1 - \sqrt{2}/2 \approx 0.2929$ gobierna todo el esquema. La forma semidiscreta tiene dos etapas.

$W^{(1)}$ es la solución de la etapa intermedia, y tanto la parte explícita como la implícita avanzan solo una fracción $\beta$.

$\delta = 1 - 1/(2\beta) \approx -0.7071$ sale de la condición de segundo orden de la tabla explícita, y que sea negativo se convierte más adelante en un problema.

Aquí el artículo pone sobre la mesa un resultado negativo.

“

No existe ningún esquema Runge–Kutta implícito de orden mayor que uno que sea TVD sin restricción sobre el paso temporal. (Resultado negativo de la línea de Gottlieb, Teorema 1 del artículo)

Este no es un muro que se atraviese con mejor código. El esquema AP de segundo orden es $L^2$ estable bajo el número CFL explícito $\sigma_e = c_e\Delta t/\Delta x \le 1$. Pero en cuanto $\Delta t$ supera el CFL acústico $\Delta x/(c_e + c_i/\sqrt{\varepsilon})$, pierde por igual la estabilidad $L^\infty$ y la propiedad TVD. El sobreimpulso que vi era exactamente eso. Acotado, así que nunca diverge. Estructural, así que dar más pasos no lo borra.

Para el análisis, el artículo reduce el sistema a un problema modelo escalar (ecuación 12 del artículo). Basta con imaginar un pulso acústico montado sobre un flujo lento.

$w$ es la incógnita escalar, $c_e$ es la velocidad convectiva lenta y $c_i/\sqrt{\varepsilon}$ es la velocidad acústica rápida. Es una estructura que traslada tal cual el escalado $1/\sqrt{\varepsilon}$ de las ondas de presión de Euler.

Q4. ¿Mezclar primer y segundo orden da realmente TVD?#

La solución de los autores es simple. Calcular el mismo paso con el esquema AP de primer orden y con el de segundo orden, y luego combinarlos de forma convexa.

$\theta$ es el peso que carga el esquema de segundo orden. Con $\theta = 0$ queda primer orden puro; con $\theta = 1$, segundo orden puro.

La idea es la misma que la de un limitador de flujo en el espacio. La diferencia está en que la mezcla se aplica a la discretización temporal y no a la espacial. El Teorema 3 del artículo da la condición. Si $\theta = \alpha\beta/(1-\beta)$ con $\alpha \in [0,1]$, el esquema mezclado es uniformemente TVD y $L^\infty$ estable. Y eso bajo un CFL independiente del número de Mach, $\sigma_e \le \sqrt{2}$ (tomando $\alpha = 1$). De ahí queda fijado el peso máximo que puede cargar el esquema de segundo orden.

$\theta_M$ es la cota superior de la cuota de segundo orden admisible manteniendo la propiedad TVD.

Leído con honestidad, queda así. Solo se puede quemar hasta un 41% del esquema de segundo orden. Y esa cota no es una constante universal, sino un valor que viene con la discretización temporal IMEX elegida, es decir, con ARS(2,2,2). Con otra tabla IMEX, $\theta_M$ también cambia.

Con un 41% la precisión se queda corta. Por eso el artículo le monta encima un enfoque MOOD (Multi-dimensional Optimal Order Detection: un limitador a posteriori que revisa la solución una vez calculada y rehace a orden más bajo solo las celdas problemáticas). El procedimiento tiene tres pasos. Primero se calcula una solución candidata de segundo orden. Luego se buscan las celdas que violan las cotas $L^\infty$ o la condición TVD. Por último se devuelven a la solución TVD-AP únicamente esas celdas. Las regiones suaves se quedan en segundo orden pleno, y solo el entorno de las discontinuidades cae a la mezcla segura.

La curva en la que $\theta=1$ rompe la banda y $\theta_M$ deja el sobreimpulso en cero vuelve con el mismo deslizador, justo después del bloque de salida.

Q5. Reproducción en Python: las oscilaciones se van, ¿y qué se pierde?#

Con el problema modelo (ec. 12) y el esquema mezclado (ec. 17) basta para que la reproducción sea corta. Se parte de un pulso cuadrado y se miden la variación total y el sobreimpulso máximo para tres valores de $\theta$.

import numpy as np
 
BETA = 1.0 - np.sqrt(2.0) / 2.0        # ARS(2,2,2)
THETA_M = BETA / (1.0 - BETA)          # = sqrt(2) - 1
 
def solve_backward(rhs, s):
    """Resuelve (1+s)w_j - s*w_{j-1} = rhs_j con frontera periódica (upwind implícito)."""
    n = rhs.size
    A = (1.0 + s) * np.eye(n)
    A[np.arange(n), np.arange(n) - 1] -= s
    return np.linalg.solve(A, rhs)
 
def dminus(v):
    return v - np.roll(v, 1)
 
def blended_step(w, se, si, theta):
    """Ecuación (17) del artículo: theta es la cuota del esquema de segundo orden."""
    b = BETA
    star = solve_backward(w - b * se * dminus(w), b * si)          # (17a)
    rhs = (w - theta * (b - 1.0) * se * dminus(w)
             - theta * (1.0 - b) * si * dminus(star)
             - theta * (2.0 - b) * se * dminus(star)
             - (1.0 - theta) * se * dminus(w))                     # (17b)
    return solve_backward(rhs, (1.0 - theta + theta * b) * si)
 
def pulse_run(eps, theta, steps=60, n=200, cfl=0.9, ce=1.0, ci=1.0):
    dx = 1.0 / n
    x = (np.arange(n) + 0.5) * dx
    w = np.where((x > 0.25) & (x <= 0.75), 1.0, -1.0)   # pulso cuadrado, TV = 4
    se, si = ce * cfl, (ci / np.sqrt(eps)) * cfl        # dt = cfl*dx/ce
    over = 0.0
    for _ in range(steps):
        w = blended_step(w, se, si, theta)
        over = max(over, w.max() - 1.0, -1.0 - w.min())
    tv = np.abs(np.roll(w, -1) - w).sum()
    return tv, over
 
for eps in (1e-2, 1e-4):
    for theta, tag in ((0.0, "AP 1er orden"), (THETA_M, "TVD-AP      "), (1.0, "AP 2do orden")):
        tv, over = pulse_run(eps, theta)
        print(f"eps={eps:<7g} {tag} TV={tv:6.3f}  sobreimpulso máx={over:+.4f}")

La salida es la siguiente.

eps=0.01    AP 1er orden TV= 0.395  sobreimpulso máx=+0.0000
eps=0.01    TVD-AP       TV= 0.977  sobreimpulso máx=+0.0000
eps=0.01    AP 2do orden TV= 3.784  sobreimpulso máx=+0.4226
eps=0.0001  AP 1er orden TV= 0.000  sobreimpulso máx=+0.0000
eps=0.0001  TVD-AP       TV= 0.000  sobreimpulso máx=+0.0000
eps=0.0001  AP 2do orden TV= 0.006  sobreimpulso máx=+0.3776

Al mover el deslizador de θ entre $\sqrt{2}-1$ y 1, el sobreimpulso y la variación total recién impresos se recalculan sobre la curva.

El sobreimpulso desaparece exactamente como dice el artículo. Solo el AP de segundo orden atraviesa la banda, en $\pm 0.38 \sim 0.42$, mientras que el AP de primer orden y el TVD-AP quedan en cero para ambos valores de $\varepsilon$. Hasta aquí la teoría y la implementación encajan sin fisuras.

El precio es la difusión. Con $\varepsilon = 10^{-2}$, tras 60 pasos la variación total se recupera desde el 0.395 del AP de primer orden hasta el 0.977 del TVD-AP. Más del doble de mejora, pero muy lejos de la variación total inicial de 4.0. Es decir, lo que vende el TVD-AP es "emborrona menos que el primer orden", no que conserve la nitidez del segundo orden.

Con $\varepsilon = 10^{-4}$ la diferencia asoma con un solo paso. El AP de segundo orden llega a una variación total de 5.510, por encima del valor inicial 4.0, y genera un sobreimpulso de $+0.378$. El TVD-AP tiene sobreimpulso exactamente cero, pero su variación total baja hasta 1.860. El número de Courant acústico es 90, así que la difusión del upwind implícito se come todo eso en un solo paso. A los 4 pasos, los valores caen a 0.062 para el AP de primer orden y 0.068 para el TVD-AP, y a los 60 pasos el pulso se ha extinguido en la práctica para los tres esquemas.

Esta difusión es justamente la razón por la que el artículo monta el limitador MOOD. Los autores sabían que la mezcla por sí sola no puede proteger la precisión.

Al dividir otra vez los números impresos aparece el precio de la difusión. Con $\varepsilon=10^{-2}$, tras 60 pasos, la variación total 0,977 del TVD-AP es unas 2,47 veces el 0,395 del primer orden, y aun así solo 0,24 de la variación total inicial 4. Un solo paso con $\varepsilon=10^{-4}$ es más empinado. El AP de segundo orden sube de 4,0 a 5,510, un aumento de 1,51, y el sobreimpulso $+0{,}378$ es casi el $+0{,}3776$ de la salida a 60 pasos. El sobreimpulso nace en el primer paso y después no crece. El TVD-AP recorta la variación total de 4 a 1,860 en ese mismo paso, una caída de 2,14. Con un número de Courant acústico de 90, el upwind implícito se come una parte grande del ancho del pulso en un paso. A los 4 pasos, 0,062 y 0,068 son cerca del 1,6% del 4 inicial. Como $\theta_M=\beta/(1-\beta)$ y $\beta=1-\sqrt{2}/2\approx 0{,}2929$, el techo $\theta_M\approx 0{,}4142$ lo fija la tabla ARS. No es un número ajustado en una malla.

¿Se corta la cuota de segundo orden en el 41%, o solo se devuelven las celdas malas?#

Llamarlo segundo orden deja una cuota utilizable solo hasta $\theta_M=\sqrt{2}-1$. El artículo escribe que un $\theta$ por celda mejoraría el resultado y deja abierta la demostración TVD de ese caso. En los ensayos de la sección 6 el CFL es $C=0{,}9$ solo para el esquema de primer orden y $C=0{,}45$ para los otros tres, atribuido a la reconstrucción espacial de segundo orden. Leído junto a la frase de que el paso temporal se liberó del número de Mach, la ganancia práctica queda en la mitad.

El problema modelo (ecs. 12 y 17) se reproduce en 30 líneas. El Euler completo es una elíptica no lineal para $\rho^{n+1}$ en cada paso, y el artículo no da ni tolerancia de convergencia ni número de iteraciones. La conclusión admite que quedan oscilaciones pequeñas en algunos casos aunque el limitador esté activo. Al lado de PISO/SIMPLE, que ya trata la presión de forma implícita, el sitio de este artículo es donde $L^\infty$ se rompe dentro de un marco conservativo basado en densidad. La factura tiene la misma forma que el $L^\infty$ de un paso de presión semiimplícito. Cuando una tobera y una región de estancamiento comparten un dominio, el paso temporal de esa mitad es lo primero que conviene poner en duda.