La factura de un paso de tiempo 1001 veces mayor — por qué un esquema semiimplícito de presión no es L∞ estable

Al borrar la velocidad del sonido de la condición CFL el paso de tiempo queda libre, pero la oscilación en una discontinuidad permanece en lugar de decaer.

1971, el año en que la presión se separó de la velocidad del sonido#

Fue cuando Francis Harlow y Anthony Amsden publicaron el método ICE en Los Álamos. El título llevaba la frase "all flow speeds". Anunciaba que la lentitud de un código compresible a bajo número de Mach se atacaría desde el paso de tiempo y no desde la malla. La entrada más reciente de esa línea es el esquema semiimplícito basado en presión de Boscheri y Pareschi, de 2021, que resuelve las ecuaciones de Navier–Stokes tridimensionales con tercer orden en todo el rango de Mach. Lo que se gana y lo que se pierde se contrasta con los números impresos. El paso de tiempo creció hasta un factor de 1001, y la oscilación en la discontinuidad seguía ahí después de cuatro vueltas.

Conviene empezar por lo que gana. En la figura siguiente hay que bajar el deslizador del número de Mach.

Mach 0.100CFL 0.50

dt_explicit = 1.136e-2   dt_semi = 1.250e-1   gain = 11.0x   steps so far: 0 vs 0

Both panels cover the same physical time and the same 40 cells. Pull the Mach slider down and the red acoustic characteristics tilt until they are nearly horizontal — the explicit panel has to slice that window into ever thinner layers to keep them inside one cell, while the green panel keeps taking the same eight steps. The step counter is the whole argument for treating pressure implicitly.

Ambos paneles cubren el mismo tiempo físico y las mismas cuarenta celdas. Al bajar el número de Mach las características acústicas rojas se inclinan hacia la horizontal, y el panel superior empieza a partir la ventana en capas de tiempo cada vez más delgadas para mantenerlas dentro de una sola celda. El contador de pasos del panel inferior no se mueve.

Dónde queda amarrado el paso de tiempo al número de Mach#

En las ecuaciones de Navier–Stokes reescaladas, la rigidez la carga $\varepsilon$ , el cuadrado del número de Mach, y la velocidad del sonido aparece como $c/\sqrt{\varepsilon}$ . Una condición CFL totalmente explícita la sigue arrastrando.

\Delta t \le \mathrm{CFL}\,\frac{\Delta x}{\max\left(|u| + c/\sqrt{\varepsilon}\right)}

$\Delta x$ es el espaciado de malla, $u$ la velocidad del fluido y $c$ la velocidad del sonido de referencia. Cuando $\varepsilon \to 0$ el denominador diverge y $\Delta t$ cae como $\sqrt{\varepsilon}$ . Un flujo donde no ocurre nada acústicamente queda secuestrado por ondas sonoras que a nadie le interesan.

El esquema semiimplícito traslada solo los términos de presión al lado implícito. La velocidad del sonido desaparece entonces de la condición de estabilidad.

\Delta t \le \mathrm{CFL}\,\frac{\Delta x}{\max |u|}

La razón entre ambas es $1 + 1/M$ . En $M = 0.001$ eso es $1 + 1000 = 1001$ . La tabla imprime esa razón como ganancia 1001.0, con un paso explícito de $1.249\times 10^{-3}$ junto a uno semiimplícito de $1.250$ . En $M=0.1$ la ganancia es 11.0, explícito $1.136\times 10^{-3}$ frente a semiimplícito $1.250\times 10^{-2}$ .

El paso del código es $\mathrm{CFL}\,\Delta x/|U_0|$ . Con $\Delta x = 1/400 = 0.0025$ , CFL $= 0.5$ y $U_0 = 0.1$ , $dt = 0.0125$ . Un tránsito, $1/(U_0 + c) = 1/1.1$ , dividido por 0.0125 da 72.7, y la salida informa 73 pasos. El CFL acústico es $0.5 \times 1 / 0.1 = 5$ . Separar esa misma velocidad del sonido del CFL desde el integrador temporal está en el IMEX de transporte acústico de Peluchon.

Cómo la presión se convierte en una ecuación elíptica#

La semidiscretización semiimplícita de primer orden son tres líneas. La continuidad queda totalmente explícita, la ecuación de momento trata de forma implícita solo el gradiente de presión, y la ecuación de energía toma implícitamente el momento dentro del flujo de entalpía.

\frac{\rho^{n+1} - \rho^{n}}{\Delta t} + \nabla\cdot(\rho\mathbf{u})^{n} = 0

\frac{(\rho\mathbf{u})^{n+1} - (\rho\mathbf{u})^{n}}{\Delta t} + \nabla\cdot\left((\rho\mathbf{u})^{n}\otimes\mathbf{u}^{n}\right) + \frac{1}{\varepsilon}\nabla p^{n+1} = 0

$n$ es el nivel temporal actual y $n+1$ el siguiente. Al sustituir la ecuación de momento en la de energía queda una sola ecuación cuya única incógnita es $p^{n+1}$ .

\varepsilon(\rho e)^{n+1} + \varepsilon\frac{\Delta t}{2}\frac{(\rho\mathbf{u})^{n}}{\rho^{n}}\nabla p^{n+1} - \Delta t^{2}\,\nabla\cdot\left(h^{n}\nabla p^{n+1}\right) = b^{n}

$h = (\rho e + p)/\rho$ es la entalpía específica y $b^n$ es un lado derecho construido solo con cantidades explícitas. El término $\Delta t^2 \nabla\cdot(h\nabla p)$ es el que vuelve elíptica la ecuación. Para un gas ideal $(\rho e)^{n+1} = p^{n+1}/(\gamma-1)$ , de modo que el sistema es lineal; con una ecuación de estado no lineal como Redlich–Kwong se le envuelve una iteración de Newton, que en el artículo converge en cuatro pasadas como máximo con tolerancia $10^{-10}$ .

Resuelta la presión se recupera el momento, y la energía total cierra con la parte implícita de la energía cinética. En resumen: masa → presión (elíptica) → momento → energía. La estructura coincide con un método de proyección incompresible, salvo que aquí la ecuación de estado sigue viva y la compresibilidad sobrevive.

¿La divergencia realmente tiende a cero en el límite de bajo Mach?#

La preservación asintótica significa que el esquema discreto reproduce por sí solo el límite incompresible cuando $\varepsilon \to 0$ . Al desarrollar cada variable en potencias de $\varepsilon$ y quedarse con los términos dominantes, la presión de orden principal resulta constante en el espacio.

\nabla p_{(0)}^{n+1} = 0

Introducirlo en el límite de la ecuación de energía da, para un gas ideal,

h^{n}\,\nabla\cdot\mathbf{u}^{n+1} = O(\Delta t)

y la constante escondida en $O(\Delta t)$ no depende de $\varepsilon$ . Por más que se baje el número de Mach, la divergencia del campo de velocidad se mantiene en el nivel de $\Delta t$ . Por eso la presión no se parte en un tablero de ajedrez al llevar la misma malla al régimen de bajo Mach.

Cuatro vueltas del mismo pulso en Python#

La sección 3.6 del artículo deja una frase incómoda. Las discretizaciones temporales implícitas de alto orden no consiguen eliminar el sobreimpulso frente a una discontinuidad fuerte, y esas oscilaciones ni crecen ni desaparecen. Decir que un esquema implícito es estable significa estable en $L_2$ , no en $L_\infty$ .

Verificar solo esa afirmación no exige el solucionador no lineal completo. Basta aplicar la misma partición a las ecuaciones de Euler unidimensionales linealizadas: transporte material explícito, acoplamiento acústico implícito, con un peso de implicitud $\theta$ barrido desde 0.5 (segundo orden, sin amortiguamiento) hasta 1.0 (primer orden, amortiguado). Colocar la presión en el centro de celda y la velocidad en las caras convierte la ecuación de presión exactamente en un sistema de Helmholtz de tres puntos. Las fronteras son periódicas, así que se usa un solucionador tridiagonal cíclico.

El estado inicial es un pulso acústico cuadrado que viaja hacia la derecha. La solución exacta es ese pulso trasladado en $u_0 + c$ , de modo que el error es analítico. Con un CFL material de 0.5, el CFL acústico queda en 5.

import numpy as np
 
GAMMA, RHO0, C0 = 1.4, 1.0, 1.0
P0 = C0 ** 2 * RHO0 / GAMMA          # so that c = sqrt(gamma p0 / rho0) = 1
N, CFL, U0 = 400, 0.5, 0.1
DX = 1.0 / N
 
 
def cyclic_thomas(a, b, c, d):
    """resolución tridiagonal periódica: Thomas + corrección de Sherman-Morrison."""
    n = len(b)
    al, be, g = a[0], c[-1], -b[0]
    bb = b.copy()
    bb[0] -= g
    bb[-1] -= al * be / g
 
    def tri(rhs):
        cp, dp = np.zeros(n), np.zeros(n)
        cp[0], dp[0] = c[0] / bb[0], rhs[0] / bb[0]
        for i in range(1, n):
            m = bb[i] - a[i] * cp[i - 1]
            cp[i], dp[i] = c[i] / m, (rhs[i] - a[i] * dp[i - 1]) / m
        x = np.zeros(n)
        x[-1] = dp[-1]
        for i in range(n - 2, -1, -1):
            x[i] = dp[i] - cp[i] * x[i + 1]
        return x
 
    v = np.zeros(n)
    v[0], v[-1] = g, be
    y, z = tri(d), tri(v)
    fac = (y[0] + al * y[-1] / g) / (1.0 + z[0] + al * z[-1] / g)
    return y - fac * z
 
 
def grad_face(p):
    return (np.roll(p, -1) - p) / DX
 
 
def div_cell(u):
    """div_cell(grad_face(.)) es exactamente el laplaciano de tres puntos."""
    return (u - np.roll(u, 1)) / DX
 
 
def upwind(f):
    return (f - np.roll(f, 1)) / DX
 
 
def si_update(u, p, dt, theta):
    """transporte material explícito, acoplamiento acústico con peso implícito theta."""
    ps = p - dt * U0 * upwind(p)
    ue = u - dt * U0 * upwind(u) - dt * (1 - theta) / RHO0 * grad_face(p)
    k = C0 ** 2 * (theta * dt) ** 2 / DX ** 2
    rhs = ps - dt * GAMMA * P0 * (theta * div_cell(ue) + (1 - theta) * div_cell(u))
    p_n = cyclic_thomas(np.full(N, -k), np.full(N, 1 + 2 * k), np.full(N, -k), rhs)
    return ue - dt * theta / RHO0 * grad_face(p_n), p_n
 
 
def chi_flattener(u, k1):
    """Ec. (104): marca las celdas que se comprimen más rápido que k1 * c_min."""
    divu = div_cell(u)
    return np.minimum(1.0, np.maximum(0.0, -(divu + k1 * C0) / (k1 * C0)))
 
 
def exact_wave(x, t):
    s = np.mod(x - (U0 + C0) * t, 1.0)
    return np.where((s >= 0.3) & (s < 0.6), 1.0, 0.0)
 
 
def solve_case(n_transit, theta, k1=None):
    xc, xf = (np.arange(N) + 0.5) * DX, (np.arange(N) + 1.0) * DX
    p, u = exact_wave(xc, 0.0), exact_wave(xf, 0.0) / (RHO0 * C0)
    t_end, dt = n_transit / (U0 + C0), CFL * DX / abs(U0)
    t, nstep, nflag = 0.0, 0, 0
    while t < t_end - 1e-14:
        h = min(dt, t_end - t)
        u2, p2 = si_update(u, p, h, theta)
        if k1 is not None:
            chi = chi_flattener(u, k1)
            nflag += int(np.sum(chi > 0))
            u1, p1 = si_update(u, p, h, 1.0)
            cf = np.maximum(chi, np.roll(chi, -1))
            u2, p2 = cf * u1 + (1 - cf) * u2, chi * p1 + (1 - chi) * p2
        u, p, t, nstep = u2, p2, t + h, nstep + 1
    ex = exact_wave(xc, t_end)
    return dict(l1=float(np.mean(np.abs(p - ex))), over=float(np.max(p) - 1.0),
                under=float(-np.min(p)), l2=float(np.sqrt(np.mean(p ** 2))),
                steps=nstep, flagged=nflag, acfl=CFL * C0 / abs(U0))
 
 
print("=== (1) time step from Eq.(20) vs Eq.(106), CFL=0.5, 400 cells ===")
print(f"{'Mach':>8}{'dt explicit':>14}{'dt semi-impl':>15}{'gain':>8}")
for M in [0.5, 0.2, 0.1, 0.02, 0.005, 0.001]:
    de, di = CFL * DX / (M * C0 + C0), CFL * DX / (M * C0)
    print(f"{M:8.3f}{de:14.3e}{di:15.3e}{di / de:8.1f}")
 
b = solve_case(1, 1.0)
print(f"\n=== (2) square acoustic pulse, M={U0}, material CFL={CFL}")
print(f"    -> acoustic CFL {b['acfl']:.1f}, {b['steps']} steps per transit ===")
print(f"{'theta':>7}{'over@1':>10}{'over@2':>10}{'over@4':>10}{'L1@4':>9}{'|p|2@4':>9}")
for th in (0.50, 0.52, 0.55, 0.60, 0.70, 1.00):
    r = [solve_case(n, th) for n in (1, 2, 4)]
    ov = [f"{x['over']:10.4f}" if abs(x['over']) < 1e3 else f"{'blow-up':>10}" for x in r]
    l1 = f"{r[2]['l1']:9.4f}" if r[2]['l1'] < 1e3 else f"{'-':>9}"
    l2 = f"{r[2]['l2']:9.4f}" if r[2]['l2'] < 1e3 else f"{'-':>9}"
    print(f"{th:7.2f}{ov[0]}{ov[1]}{ov[2]}{l1}{l2}")
 
print("\n=== (3) flattener blend, theta=0.55 against theta=1, 4 transits ===")
print(f"{'k1':>9}{'flagged/step':>14}{'over':>9}{'under':>9}{'L1':>9}")
base = solve_case(4, 0.55)
print(f"{'off':>9}{0:14.1f}{base['over']:9.4f}{base['under']:9.4f}{base['l1']:9.4f}")
for k1 in [1e-3, 1e-1, 1e0, 1e1, 1e2]:
    r = solve_case(4, 0.55, k1)
    print(f"{k1:9.0e}{r['flagged'] / r['steps']:14.1f}{r['over']:9.4f}"
          f"{r['under']:9.4f}{r['l1']:9.4f}")
mono = solve_case(4, 1.0)
print(f"{'theta=1':>9}{'(all cells)':>14}{mono['over']:9.4f}"
      f"{mono['under']:9.4f}{mono['l1']:9.4f}")

=== (1) time step from Eq.(20) vs Eq.(106), CFL=0.5, 400 cells ===
    Mach   dt explicit   dt semi-impl    gain
   0.500     8.333e-04      2.500e-03     3.0
   0.200     1.042e-03      6.250e-03     6.0
   0.100     1.136e-03      1.250e-02    11.0
   0.020     1.225e-03      6.250e-02    51.0
   0.005     1.244e-03      2.500e-01   201.0
   0.001     1.249e-03      1.250e+00  1001.0
 
=== (2) square acoustic pulse, M=0.1, material CFL=0.5
    -> acoustic CFL 5.0, 73 steps per transit ===
  theta    over@1    over@2    over@4     L1@4   |p|2@4
   0.50  488.1288   blow-up   blow-up        -        -
   0.52    4.5456  212.3962   blow-up        -        -
   0.55    0.2285    0.2023    0.1825   0.1348   0.5169
   0.60    0.0538    0.0286   -0.0346   0.1270   0.4630
   0.70   -0.0163   -0.0900   -0.2320   0.1992   0.4064
   1.00   -0.1610   -0.3210   -0.4823   0.3098   0.3347
 
=== (3) flattener blend, theta=0.55 against theta=1, 4 transits ===
       k1  flagged/step     over    under       L1
      off           0.0   0.1825   0.1663   0.1348
    1e-03         302.9  -0.3611   0.4734   0.2870
    1e-01         280.0  -0.3544   0.4998   0.2888
    1e+00         196.0  -0.1857   0.5381   0.2448
    1e+01           8.2   0.1275   0.3135   0.1469
    1e+02           0.0   0.1801   0.1726   0.1355
  theta=1   (all cells)  -0.4823  -0.0985   0.3098

Con el deslizador de $\theta$ en 0.55, el sobreimpulso tras cuatro vueltas se queda en 0.1825. Al bajarlo a 0.50, el pulso revienta en 488 en una sola vuelta.

theta 0.55k1 = 1e1

transits 0.00   overshoot 0.0000   undershoot 0.0000   L1 0.0000   flagged 0 / 128

The acoustic CFL here is 5, so the implicit solve is doing all the work. Leave theta at 0.55 and watch the overshoot readout while the pulse laps the domain: it settles near 0.18 and stays there — bounded, but never shrinking. Push theta to 1.00 and the wiggles vanish while the pulse melts. Turning on the flattener with k1 = 1e-3 lights up the yellow ticks along almost the whole axis, which is the scheme quietly dropping to first order everywhere.

Al empujar el deslizador hasta 1.00 el sobreimpulso se vuelve negativo y $L_{1}$ crece hasta 0.3098. El 0.1348 de $\theta=0.55$ y ese extremo lejano están en el mismo deslizador.

Lo que separa a 0.55 de 0.50#

La segunda tabla es donde se parten las corridas. $\theta = 0.5$ es Crank–Nicolson, formalmente de segundo orden, y sobre la parte acústica sola conserva la energía de manera exacta. Combinado con el transporte explícito, el sobreimpulso llega a 488 en una vuelta y la corrida muere en la segunda. Cuando la mitad implícita no aporta amortiguamiento alguno, el amortiguamiento aguas arriba de la mitad explícita no alcanza a cubrir el error de partición. Por eso los esquemas de la línea de Casulli recomiendan tomar $\theta$ algo por encima de 0.5.

En $\theta = 0.55$ el panorama cambia. El sobreimpulso recorre 0.2285 → 0.2023 → 0.1825. De 0.2285 a 0.2023 hay 0.0262. El siguiente tramo, hasta 0.1825, es 0.0198. Tras cuatro vueltas, 0.1825 es el 18% de un pulso de altura 1. No crece, así que el esquema es estable en $L_{2}$ . Tampoco se encoge, así que no es estable en $L_{\infty}$ . La frase del artículo se reproduce sin retoques. En el otro extremo, $\theta = 1.0$ tiene sobreimpulso negativo: el pulso ni siquiera alcanza su propia altura. $L_{1}$ es 0.1348 frente a 0.3098. $0.3098/0.1348 \approx 2.30$ , el factor 2.3 ya dicho. $|p|_{2}$ baja, de 0.5169 a 0.3347. En $L_{2}$ , 1.0 es la corrida más quieta. En el pico solo, 0.55 se derrite menos.

El $k_1$ del flattener nunca fue un umbral#

El remedio del artículo es detectar y luego mezclar. Las celdas con choque se encuentran comparando la divergencia del campo de velocidad con la menor velocidad del sonido entre los vecinos.

\chi^{n} = \min\left[1,\ \max\left(0,\ -\frac{\nabla\cdot\mathbf{u}^{n} + k_1 c^{n}_{\min}}{k_1 c^{n}_{\min}}\right)\right]

Donde $\chi$ vale cero se conserva el candidato de alto orden, donde vale uno lo reemplaza la solución de primer orden, y en medio ambas se combinan de forma convexa. El artículo usó $k_1 = 10^{-3}$ en todos los cálculos e informó que se marcan muy pocas celdas.

La tercera tabla muestra qué ocurre al trasladar esa constante. Con $k_1 = 10^{-3}$ se marcan 302.9 de 400 celdas por paso. El sobreimpulso sí desaparece, de 0.1825 a $-0.3611$ , pero el error $L_1$ pasa de 0.1348 a 0.2870, un factor de 2.1. Eso no es detección: es todo el esquema bajando en silencio a primer orden, y aterriza casi donde aterriza $\theta = 1$ (0.3098).

La razón está en la fórmula. $\chi$ sube de 0 a 1 mientras $\nabla\cdot\mathbf{u}$ va de $-k_1 c$ a $-2k_1 c$ . Esa rampa mide un solo $k_1 c$ . Con $k_1 c = 10^{-3}$ el ancho es 0.001. Con $k_1 = 10^{-3}$ , $302.9/400 = 0.757$ , cerca del 76% de las celdas. $L_{1}$ es $0.2870/0.1348 \approx 2.13$ , el factor 2.1 ya dicho. En este problema la divergencia fluctúa a nivel $O(1)$ , de modo que cualquier celda que esté comprimiendo aunque sea un poco cae en $\chi = 1$ . Actúa como una prueba de signo, no como un umbral.

Seleccionar celdas de verdad exigió cuatro órdenes de magnitud más. Con $k_1 = 10$ se marcan apenas 8.2 celdas por paso, y el sobreimpulso pasa de 0.1825 a 0.1275. $0.1275/0.1825 \approx 0.70$ , cerca de un 30 por ciento menos, y $0.1469/0.1348 \approx 1.09$ , así que $L_{1}$ crece solo un 9 por ciento. Con $k_1 = 10^2$ no se marca nada. Nada de esto vuelve incorrecto el valor publicado: está ajustado a tubos de choque no lineales donde la divergencia es prácticamente nula en las zonas suaves. Al llevarlo a un problema que ya carga oscilaciones, la referencia pierde sentido.

Qué queda donde antes estaba la velocidad del sonido#

Borrar la velocidad del sonido de la condición CFL deja a la integración temporal sin poder seguir ese sonido. A bajo Mach esa es la decisión correcta. Al subir el número de Mach, esas ondas son choques. $\theta$ es una perilla que cambia precisión por monotonía. Bajo esta partición, 0.5 ni siquiera resultó estable: en una vuelta el sobreimpulso ya está en 488.

La constante del detector está atada a la escala de divergencia del problema. Copiar el $10^{-3}$ publicado marca 302.9 de 400 celdas. El detector se vuelve un interruptor. Otra partición que muda la velocidad del sonido al subsistema de presión, en lugar de borrarla, está en el lugar donde Toro–Vázquez traslada la rigidez. Allí la velocidad del sonido no desaparece. Cambia de ecuación.