[Reseña] Sospechar de la celda después de resolverla

En el test 3 de Toro el MUSCL sin limitador da presión negativa hacia el paso 14; MOOD recalcula solo cuatro a seis celdas del choque en primer orden.

Al reimplementar este artículo, lo primero que descolocó fue el orden. Casi todo el mundo intercala el limitador (limiter, postproceso que reduce las oscilaciones) antes de que la celda se descontrole. Lo hizo Van Leer, lo hizo Sweby, lo hace toda la familia TVD. Ioriatti, Dumbser y Loubère (2020) invierten la línea de tiempo en su esquema DG semi-implícito sobre malla escalonada — primero se resuelve hasta el final y luego se sospecha. Solo después de que aparezca una densidad negativa la celda se aísla y se recalcula con un esquema de volúmenes finitos de primer orden. Este texto recorre ese algoritmo MOOD (Multi-dimensional Optimal Order Detection) y reproduce el mapa de troubled cells en el Toro Test 3, donde la presión salta de 1000 a 0.01.

Al mover la relación de presiones y $\delta_0$ cambia el número de celdas rojas junto al choque. Con overlay unlimited candidate activo, el candidato sin limitador aparece en trazo discontinuo.

1D Riemann problem · candidate (MUSCL) → MOOD detection → fallback (LF) in troubled cells

pressure ratio p_L / p_R = 1000 : 0.01relaxation delta_0 = 1e-4

overlay unlimited candidate (red dashed)

Cyan = density, yellow = pressure under MOOD. Red bands mark cells that failed PAD or NAD this step and were recomputed with the Lax–Friedrichs fallback. Toggle the overlay to see the unlimited MUSCL candidate crash near the shock once the pressure ratio gets large.

Con 1000:0.01 hay cuatro a seis celdas troubled. Con $\delta_0=10^{-6}$ también saltan las zonas suaves. Con $10^{-2}$ el choque se cuela. La banda útil va de $10^{-4}$ a $10^{-3}$ .

Datos del artículo#

Autores: Matteo Ioriatti, Michael Dumbser, Raphaël Loubère
Revista: Journal of Scientific Computing, 83(2), 2020
DOI: 10.1007/s10915-020-01209-w
Título: A Staggered Semi-implicit Discontinuous Galerkin Scheme with a Posteriori Subcell Finite Volume Limiter for the Euler Equations of Gasdynamics

En una línea — limitador a priori frente a a posteriori#

Los limitadores clásicos cobran impuesto a todas las celdas. Sospechan de cada pendiente del polinomio incluso en zonas suaves, y la recortan. MOOD solo cobra impuesto a las celdas que se portan mal. Mantiene el polinomio de alto orden donde la solución candidata parece razonable y baja a primer orden únicamente donde algo ha ido mal. La consecuencia es decisiva: un DG de alto orden conserva su precisión formal en la mayor parte de la malla aun con un choque fuerte encima.

La idea desciende de la familia MOOD a posteriori de Clain, Loubère y Diot (2011); Dumbser et al. (2014) la trasladaron a DG explícito y el artículo aquí reseñado la lleva al DG semi-implícito escalonado. La parte semi-implícita tiene además su propio beneficio: el CFL queda gobernado por la velocidad de masa, no por la velocidad del sonido. Un mismo código sirve para flujos a bajo Mach y para choques fuertes.

Dos representaciones sobre dos mallas — polinomio y promedios de subcelda#

Dentro de cada celda principal $T_i$ conviven dos representaciones de los datos.

Grados de libertad DG $\hat{q}_i^n$ — los $P+1$ coeficientes de un polinomio de grado $P$ .
Promedios de subcelda $\bar{q}_{i,s}^n$ — los $2P+1$ promedios de volumen finito de primer orden.

Dos operadores transportan información entre ambas representaciones.

\bar{q}_{i,s}^n = \mathcal{P} \cdot \hat{q}_i^n, \qquad \hat{q}_i^n = \mathcal{W} \cdot \bar{q}_i^n, \qquad \mathcal{W}\cdot\mathcal{P} = \mathcal{I}

$\mathcal{P}$ es la proyección $L_2$ del polinomio sobre los promedios de subcelda. $\mathcal{W}$ es una recuperación por mínimos cuadrados con restricción (constrained least-squares, regresión que fija como restricción la conservación de la integral sobre la celda grande). Como $2P+1 > P+1$ , la recuperación es sobre-determinada y la restricción es necesaria.

Cuando una celda se marca como troubled, su representación pasa de polinomio a promedios de subcelda. En las interfaces entre una celda problemática y otra que no lo es aparece una interfaz mixta; los operadores $\mathcal{W}_L, \mathcal{W}_R$ recuperan el polinomio a partir de la mitad relevante de la celda y mantienen consistente la integral en la interfaz. Que las dos representaciones puedan compartir una misma celda es el activo de esta familia. Al quitar un punto de cuadratura esa convivencia se rompe, la misma divergencia del piso de cuadratura DG.

PAD y NAD — dos puestos de control#

Cuando sale la solución candidata $Q^{\circ,n+1}$ , la esperan dos puestos de control.

PAD (Physical Admissibility Detector) — ¿tiene sentido físico?

\hat{\rho}_{i,l}^{\circ,n+1} > 0, \qquad \hat{p}_{i,l}^{\circ,n+1} > 0

Una densidad o presión negativa rompe la ecuación de estado (EOS). Si la celda cae aquí, queda marcada como troubled.

NAD (Numerical Admissibility Detector) — ¿está demasiado lejos de las vecinas? Aquí entra un principio del máximo discreto relajado (relaxed discrete maximum principle, DMP).

\min_{T_j \in \mathcal{V}_i} \hat{q}_{j,l}^n - \delta \le \hat{q}_{i,l}^{\circ,n+1} \le \max_{T_j \in \mathcal{V}_i} \hat{q}_{j,l}^n + \delta

con $\delta = \max[\delta_0, \epsilon (q_{\max}-q_{\min})]$ , $\delta_0 \in [10^{-4}, 10^{-3}]$ y $\epsilon = 5\cdot 10^{-4}$ . La relajación es la clave: un DMP estricto marcaría como troubled incluso los máximos locales benignos de una función suave, y el limitador se dispararía sin parar. La banda relajada le concede a esos máximos un pequeño respiro. $\delta_0$ es un margen absoluto; $\epsilon$ proporciona el margen relativo al rango local. Las celdas que fallen cualquiera de los dos controles reciben el indicador troubled $\beta_i = 1$ .

Python — un paso de MOOD#

El artículo emplea DG, pero la coreografía candidato → control → recuperación se traslada igual de bien a un MUSCL clásico. Sesenta líneas comprimen la idea.

import numpy as np
 
GAMMA = 1.4
 
def state_to_primitive(U, gamma=GAMMA):
    rho = U[0]
    u = U[1] / np.maximum(rho, 1e-12)
    e = U[2] / np.maximum(rho, 1e-12) - 0.5 * u * u
    p = (gamma - 1.0) * rho * e
    return rho, u, p
 
def physical_flux_euler(U, gamma=GAMMA):
    rho, u, p = state_to_primitive(U, gamma)
    return np.array([rho * u, rho * u * u + p, u * (U[2] + p)])
 
def lax_friedrichs_flux(UL, UR, gamma=GAMMA):
    rL, uL, pL = state_to_primitive(UL, gamma)
    rR, uR, pR = state_to_primitive(UR, gamma)
    a = max(abs(uL) + np.sqrt(gamma * pL / rL),
            abs(uR) + np.sqrt(gamma * pR / rR))
    return 0.5 * (physical_flux_euler(UL) + physical_flux_euler(UR)) \
           - 0.5 * a * (UR - UL)
 
def candidate_step_muscl(U, dx, dt, gamma=GAMMA):
    # Candidato P=1: pendiente centrada + flujo LF (sin limitador a propósito)
    N = U.shape[1]
    slope = np.zeros_like(U)
    slope[:, 1:-1] = 0.5 * (U[:, 2:] - U[:, :-2])
    UL = U + 0.5 * slope          # estado izquierdo en la cara derecha
    UR = U - 0.5 * slope          # estado derecho en la cara izquierda
    F = np.zeros((3, N + 1))
    for i in range(1, N):
        F[:, i] = lax_friedrichs_flux(UL[:, i - 1], UR[:, i], gamma)
    out = U.copy()
    out[:, 1:-1] = U[:, 1:-1] - dt / dx * (F[:, 2:N] - F[:, 1:N - 1])
    return out
 
def detect_troubled_cells(U_old, U_new, delta0=1e-4, eps=5e-4):
    rho_n, _, p_n = state_to_primitive(U_new)
    rho_o, _, p_o = state_to_primitive(U_old)
    beta = (rho_n <= 0) | (p_n <= 0) | ~np.isfinite(rho_n) | ~np.isfinite(p_n)
    for q_o, q_n in [(rho_o, rho_n), (p_o, p_n)]:
        for i in range(1, len(beta) - 1):
            qm, qi, qp = q_o[i - 1], q_o[i], q_o[i + 1]
            qmin, qmax = min(qm, qi, qp), max(qm, qi, qp)
            d = max(delta0, eps * (qmax - qmin))
            if q_n[i] < qmin - d or q_n[i] > qmax + d:
                beta[i] = True
    return beta
 
def fallback_godunov(U, dx, dt, gamma=GAMMA):
    # FV de primer orden — la red de seguridad robusta
    N = U.shape[1]
    F = np.zeros((3, N + 1))
    for i in range(1, N):
        F[:, i] = lax_friedrichs_flux(U[:, i - 1], U[:, i], gamma)
    out = U.copy()
    out[:, 1:-1] = U[:, 1:-1] - dt / dx * (F[:, 2:N] - F[:, 1:N - 1])
    return out
 
def mood_step(U, dx, dt, gamma=GAMMA, delta0=1e-4, eps=5e-4):
    U_cand = candidate_step_muscl(U, dx, dt, gamma)
    beta = detect_troubled_cells(U, U_cand, delta0, eps)
    U_safe = fallback_godunov(U, dx, dt, gamma)
    U_new = np.where(beta[None, :], U_safe, U_cand)
    return U_new, beta

El fragmento, sobre el Toro Test 3 ( $\rho_L=1$ , $p_L=1000$ , $\rho_R=1$ , $p_R=0.01$ , $N=200$ ). El MUSCL sin limitador escupe una presión negativa hacia el paso 14 y diverge. Con MOOD encima, cuatro a seis celdas alrededor del choque se encienden como troubled y se recalculan con el flujo de primer orden, mientras las restantes ~194 celdas conservan la precisión P=1 del candidato.

A la izquierda $\rho=1$ , $p=1000$ . A la derecha $\rho=1$ , $p=0.01$ . $N=200$ . El paso 14 es el momento en que el MUSCL sin limitador emite presión negativa. Solo 4–6 celdas vuelven al FV de primer orden, y unas 194 conservan P=1. Hay $2P+1$ subceldas y $P+1$ grados de libertad. Con P=1 son 3 subceldas y 2 coeficientes, así que la recuperación $\mathcal{W}$ queda sobredeterminada. Por eso la integral de la celda entra como restricción.

La banda es $\delta=\max[\delta_0,\epsilon(q_{\max}-q_{\min})]$ , con $\delta_0$ de $10^{-4}$ a $10^{-3}$ y $\epsilon=5\cdot 10^{-4}$ . Con $\delta_0=10^{-6}$ también se marcan zonas suaves. Con $10^{-2}$ el choque se escapa. En 200 celdas y CFL $=0.4$ , una o dos celdas del plateau parpadeaban en cada paso. La solución apenas se movió. El coste sí: cerca del 10% por detección en la reproducción, frente a cerca del 5% del artículo, porque $\mathcal{L},\mathcal{C},\mathcal{R}$ se rearman cada vez.

De 200 celdas, 4–6 quedan troubled y unas 194 conservan P=1. PAD solo pregunta si la densidad y la presión siguen positivas. La presión negativa del paso 14 es un candidato que no cruza esa puerta. En la banda NAD, $\epsilon=5\cdot 10^{-4}$ es el ancho relativo y $\delta_0$ de $10^{-4}$ a $10^{-3}$ es el ancho absoluto. Cuando la presión cae de 1000 a 0.01, $q_{\max}-q_{\min}$ es grande y el término con $\epsilon$ puede superar a $\delta_0$ . La densidad vale 1 a ambos lados, así que el mismo $\delta_0$ resulta más estricto ahí. El paso de cerca del 5% a cerca del 10% es el bloque 4×4 al cambiar de DG a FV: pocas celdas troubled, pero el patrón de la matriz igual se sustituye. $\beta_i$ se reinicia en cada paso. Un choque lento vuelve a hacer sospechosas las mismas 4–6 celdas.

Visualización — dónde se encienden las troubled cells#

Al ensanchar la banda $\delta$ alrededor del mínimo y el máximo de las tres vecinas, se estrecha el rango en el que un candidato queda marcado.

Relaxed DMP — admissible band around three neighbour values

q_(i-1) = 1.00q_i = 1.00q_(i+1) = 0.10candidate q* = 0.55delta_0 = 1e-4epsilon = 1e-3

Dark band = strict DMP interval [min, max] of three neighbours. Lighter band = relaxed by delta = max(delta_0, epsilon · (max − min)). Drag q* to feel where the candidate is flagged TROUBLED. Push delta_0 down to 1e-6 and even a small local peak triggers; push it up to 1e-2 and a real shock can sneak past.

Esa banda es el intervalo admisible. Con $\delta$ en cero, cualquier extremo local queda troubled. La relajación es el margen que deja vivo un pico suave.

Lo que asoma en la reproducción#

Falsos positivos del limitador. El propio artículo admite que en RP4 aparecen "troubled cells extra en la zona del plateau, probablemente por ruido numérico". Al reproducirlo en una malla de 200 celdas con CFL = 0.4 hay una o dos celdas del plateau que parpadean en casi cada paso. No afecta al resultado, pero el coste medido del limitador queda un poco por encima del reportado.

Ajuste de $\delta_0$ por variable. La presión tiene una amplitud mucho mayor que la densidad en este test, así que un único $\delta_0$ resulta generoso con la presión y estricto con la densidad. El artículo usa un valor común; en la configuración 2 de Riemann 2D, normalizar variable a variable nos dio un detector más limpio.

Coste de recomposición del sistema implícito. Cuando se detecta una troubled cell hay que reconstruir los bloques $\mathcal{L}, \mathcal{C}, \mathcal{R}$ del sistema lineal de presión. Los bloques 4×4 cambian de forma al pasar de DG a FV, así que también hay que cambiar el patrón de sparsidad. En nuestra reproducción añadimos un ~10% por detección — casi el doble del ~5% que reporta el artículo.

Lo que MOOD deja abierto#

El algoritmo estandariza explícitamente el bucle detección → recálculo. Pero el indicador troubled $\beta_i$ se reinicia en cada paso. Un choque que se mueve despacio fuerza a marcar las mismas celdas una y otra vez. Un indicador persistente podría reducir el propio coste del limitador. La continuación del grupo de Dumbser (Boscheri 2024, ALE-DG semi-implícito) da el primer paso concreto en esa dirección.

Hoy tocamos un punto muy concreto donde dos representaciones conviven dentro de la misma celda: polinomio donde todo es suave, promedios de subcelda donde no lo es. El esqueleto de los próximos resolutores all-Mach se sostiene exactamente sobre este pequeño conmutador. Cortar una densidad negativa dentro del propio flujo es el limitador de positividad. MOOD duda primero del candidato y devuelve a primer orden solo las celdas rotas.

[Reseña] Sospechar de la celda después de resolverla — limitador MOOD subcell para DG