El mismo código dio 0.79 y 2.00 a la vez — qué miden realmente los casos de verificación

Un orden medido sobre una solución que contiene una discontinuidad no es el orden del esquema. El orden se mide solo en soluciones suaves; el tubo de Sod sirve para revisar la posición del choque y la monotonía.

Un esquema de segundo orden reportó 0.79#

La lista de verificación del código compresible tenía tres entradas: tubo de choque de Sod, convección de vórtice isentrópico, interacción choque–vórtice. Empecé por Sod, en ese orden.

Refiné la malla de 100 a 800 celdas en cuatro niveles y medí el error $L_1$ en densidad. El orden de convergencia salió 0.79. Reconstrucción MUSCL, flujo HLLC, Runge–Kutta SSP de dos etapas en el tiempo. El orden de diseño es 2.

Pasé dos días revisando la reconstrucción y las condiciones de frontera. No había ningún error. El 0.79 es el comportamiento correcto.

Este artículo muestra por qué, con números. También observa al mismo código reportar 2.001 en otra prueba y ordena qué puede medir cada uno de los tres casos de verificación.

El orden de convergencia es una propiedad de la solución, no del esquema#

Las demostraciones de orden salen de un desarrollo de Taylor. Usarlo exige que la solución sea diferenciable esa cantidad de veces. En una discontinuidad, $\rho'$ no existe. Si la premisa se rompe, la conclusión también.

Así que "esquema de segundo orden" es una abreviatura de "segundo orden en soluciones suaves". Al quitar la aclaración, la frase se vuelve falsa. Al subir N solo se estrecha el ancho del emborronado; con show Linf la misma solución se lee con otra pendiente.

N = 100 smear = 3 cells t = 0.200 
at N = 100: L1 0.00e+0  ·  Linf 0.00e+0  ·  max error sits at x = 0.000

El ancho de la banda difuminada se estrecha con la malla. La altura del salto no se mueve. La pendiente en $L_1$ y la pendiente en $L_\infty$ se separan en la misma pantalla. La solución numérica se modela como el perfil exacto difuminado sobre un ancho proporcional al tamaño de celda.

El área que una discontinuidad difuminada deja en $L_1$ #

Un esquema de captura de choques resuelve una discontinuidad repartiéndola en unas pocas celdas. El número de celdas $\alpha$ se mantiene casi fijo al refinar. El ancho físico de difusión es $\alpha h$ .

\|e\|_{1} = \int_0^1 |\rho_h - \rho|\,\mathrm{d}x \approx |\Delta\rho|\,\alpha h

Aquí $\Delta\rho$ es el salto de densidad a través de la discontinuidad, $\alpha$ el número de celdas que cubre la difusión y $h$ el tamaño de celda. La altura queda fija y solo el ancho se reduce. El área escala con $h$ . Eso es primer orden.

La solución de Sod lleva dos discontinuidades: el choque y el contacto. Se difuminan de manera distinta.

Un choque tiene un mecanismo compresivo que lo vuelve a afilar. Su ancho de difusión queda sujeto a unas pocas celdas y da $O(h)$ . Un contacto no tiene ese mecanismo. La difusión numérica se acumula con el tiempo y el ancho se reduce más despacio que $h$ . La estimación conocida para esquemas monótonos ronda $O(h^{2/3})$ .

El 0.79 medido es la mezcla de ambos. Queda entre 2/3 y 1.

Por qué $L_\infty$ se detuvo en 0.083#

Al leer el mismísimo cálculo en norma máxima, el número deja de moverse por completo.

\|e\|_{\infty} = \max_i |\rho_{h,i} - \rho(x_i)| \approx \beta\,|\Delta\rho|

En algún punto de la región difuminada siempre hay una celda donde el valor numérico y el exacto difieren en una fracción fija de $|\Delta\rho|$ . Esa fracción $\beta$ no depende de la malla.

Los valores medidos fueron 0.0678 → 0.0830 → 0.0834 → 0.0857. La columna impresa es $6.7837\times 10^{-2}$ , $8.3037\times 10^{-2}$ , $8.3428\times 10^{-2}$ , $8.5667\times 10^{-2}$ , con órdenes −0.292, −0.007, −0.038. En N=200 el error ya está pegado cerca de 0.083, y refinar hasta 800 no lo pasa de 0.0857. Un ancho de difusión clavado en unas pocas celdas ocupa el mismo sitio que la viscosidad artificial escrita para emborronar un choque. Allí un término viscoso fijaba el ancho. Aquí la norma fija la altura.

La columna $L_1$ se mueve de otra manera. En N=100 vale $6.5089\times 10^{-3}$ , en 200 vale $3.7579\times 10^{-3}$ con orden 0.792, en 400 vale $2.1461\times 10^{-3}$ con 0.808 y en 800 vale $1.2420\times 10^{-3}$ con 0.789. Oscila cerca de 0.79 y no supera 1. En la onda suave, el $L_1$ sin limitador en las mismas cuatro mallas es $2.0251\times 10^{-3}$ , $5.0589\times 10^{-4}$ (2.001), $1.2630\times 10^{-4}$ (2.002), $3.1562\times 10^{-5}$ (2.001). Con minmod encendido la columna es $5.6574\times 10^{-3}$ , $2.0125\times 10^{-3}$ (1.491), $5.7467\times 10^{-4}$ (1.808), $1.5736\times 10^{-4}$ (1.869). En N=100 el cociente de esos dos errores es el 2.8 ya citado. En N=800 siguen lado a lado el valor sin limitador $3.1562\times 10^{-5}$ y el valor con limitador $1.5736\times 10^{-4}$ . El orden 1.869 es la subida de vuelta hacia 2, no un cambio del orden de diseño.

El mismo código, el mismo esquema, una solución suave#

Las ecuaciones de Euler compresibles también admiten una solución exacta suave fácil de escribir: presión uniforme, velocidad uniforme y una perturbación solo en densidad. Es una onda de contacto pura, así que la forma viaja sin cambios.

\rho(x,t) = \rho_0(x - u_0 t), \qquad u \equiv u_0, \qquad p \equiv p_0

Aquí $\rho_0$ es el perfil inicial de densidad, $u_0$ la velocidad uniforme y $p_0$ la presión uniforme. Basta con dar una vuelta completa en un dominio periódico y comparar con el campo inicial. Es el mismo papel que cumple el vórtice isentrópico en dos dimensiones.

Convecté una campana gaussiana una vuelta con $u_0 = 1$ . Con el limitador apagado y pendiente central simple, los órdenes salieron 2.001, 2.002, 2.001. El orden de diseño aparece intacto. Al empinarse la campana hacia un pulso rectangular, la línea del esquema no se toca y el orden medido se asienta hacia 1. El control de eso queda justo detrás de la tabla impresa.

Dónde minmod recorta 2.001 hasta 1.491#

Los códigos de producción corren con el limitador encendido. Minmod anula la pendiente cuando las dos diferencias tienen signos distintos.

\sigma_i = \operatorname{minmod}(\rho_i - \rho_{i-1},\; \rho_{i+1} - \rho_i)

Aquí $\sigma_i$ es la pendiente de reconstrucción en la celda $i$ . En la cima de la campana la diferencia izquierda es positiva y la derecha negativa. La condición se dispara y $\sigma_i = 0$ . Esa única celda cae a reconstrucción de primer orden.

Los órdenes medidos fueron 1.491 → 1.808 → 1.869. Van subiendo de vuelta hacia 2. El orden no cambió; el error de primer orden de un extremo simplemente todavía no se diluyó. En la misma malla, apagar el limitador deja el error 2.8 veces más pequeño.

La consecuencia práctica es una sola regla. La verificación de orden se hace con el limitador apagado o con un limitador que preserve extremos, y la producción corre con el limitador encendido. Si se mezclan las dos configuraciones en una misma tabla, ninguna fila es confiable.

Un solo solver, cuatro números#

Pasé Sod y la onda suave por el mismo solver. Las referencias son la solución exacta de Riemann y una gaussiana trasladada.

import numpy as np
G = 1.4
 
def riemann_star_pressure(rl, ul, pl, rr, ur, pr):
    al, ar = np.sqrt(G*pl/rl), np.sqrt(G*pr/rr)
    def fk(p, r, u, a, pk):
        if p > pk:
            A, B = 2.0/((G+1)*r), (G-1)/(G+1)*pk
            f = (p-pk)*np.sqrt(A/(p+B))
            df = np.sqrt(A/(B+p))*(1.0-0.5*(p-pk)/(B+p))
        else:
            f = 2*a/(G-1)*((p/pk)**((G-1)/(2*G))-1.0)
            df = 1.0/(r*a)*(p/pk)**(-(G+1)/(2*G))
        return f, df
    p = 0.5*(pl+pr)
    for _ in range(80):
        f1, d1 = fk(p, rl, ul, al, pl)
        f2, d2 = fk(p, rr, ur, ar, pr)
        dp = (f1+f2+(ur-ul))/(d1+d2)
        p = max(p-dp, 1e-10)
        if abs(dp) < 1e-13: break
    f1, _ = fk(p, rl, ul, al, pl)
    f2, _ = fk(p, rr, ur, ar, pr)
    return p, 0.5*(ul+ur)+0.5*(f2-f1)
 
def sample_exact_riemann(x, t, rl, ul, pl, rr, ur, pr):
    ps, us = riemann_star_pressure(rl, ul, pl, rr, ur, pr)
    al, ar = np.sqrt(G*pl/rl), np.sqrt(G*pr/rr)
    s = x/t
    rho = np.empty_like(s)
    for i, si in enumerate(s):
        if si <= us:                                  # a la izquierda del contacto
            if ps > pl:                               # choque a la izquierda
                rsl = rl*((ps/pl+(G-1)/(G+1))/((G-1)/(G+1)*ps/pl+1))
                ss = ul-al*np.sqrt((G+1)/(2*G)*ps/pl+(G-1)/(2*G))
                rho[i] = rl if si < ss else rsl
            else:                                     # abanico a la izquierda
                asl = al*(ps/pl)**((G-1)/(2*G))
                if si < ul-al: rho[i] = rl
                elif si > us-asl: rho[i] = rl*(ps/pl)**(1/G)
                else:
                    c = 2/(G+1)+(G-1)/((G+1)*al)*(ul-si)
                    rho[i] = rl*c**(2/(G-1))
        else:
            if ps > pr:
                rsr = rr*((ps/pr+(G-1)/(G+1))/((G-1)/(G+1)*ps/pr+1))
                ss = ur+ar*np.sqrt((G+1)/(2*G)*ps/pr+(G-1)/(2*G))
                rho[i] = rr if si > ss else rsr
            else:
                asr = ar*(ps/pr)**((G-1)/(2*G))
                if si > ur+ar: rho[i] = rr
                elif si < us+asr: rho[i] = rr*(ps/pr)**(1/G)
                else:
                    c = 2/(G+1)-(G-1)/((G+1)*ar)*(ur-si)
                    rho[i] = rr*c**(2/(G-1))
    return rho
 
def to_prim(U):
    r = U[0]; u = U[1]/r; p = (G-1)*(U[2]-0.5*r*u*u)
    return r, u, p
 
def hllc_flux(rL, uL, pL, rR, uR, pR):
    aL, aR = np.sqrt(G*pL/rL), np.sqrt(G*pR/rR)
    EL = pL/(G-1)+0.5*rL*uL**2; ER = pR/(G-1)+0.5*rR*uR**2
    sL = np.minimum(uL-aL, uR-aR); sR = np.maximum(uL+aL, uR+aR)
    sM = (pR-pL+rL*uL*(sL-uL)-rR*uR*(sR-uR))/(rL*(sL-uL)-rR*(sR-uR))
    FL = np.array([rL*uL, rL*uL**2+pL, uL*(EL+pL)])
    FR = np.array([rR*uR, rR*uR**2+pR, uR*(ER+pR)])
    UL = np.array([rL, rL*uL, EL]); UR = np.array([rR, rR*uR, ER])
    UsL = rL*(sL-uL)/(sL-sM)*np.array([np.ones_like(rL), sM, EL/rL+(sM-uL)*(sM-pL/(rL*(sL-uL)))])
    UsR = rR*(sR-uR)/(sR-sM)*np.array([np.ones_like(rR), sM, ER/rR+(sM-uR)*(sM-pR/(rR*(sR-uR)))])
    F = np.where(sL >= 0, FL, np.where(sM >= 0, FL+sL*(UsL-UL),
        np.where(sR >= 0, FR+sR*(UsR-UR), FR)))
    return F
 
def minmod_slope(q, limited):
    dL = q[1:-1]-q[:-2]; dR = q[2:]-q[1:-1]
    if not limited:
        return 0.5*(dL+dR)
    return np.where(dL*dR <= 0, 0.0, np.where(np.abs(dL) < np.abs(dR), dL, dR))
 
def residual(U, h, periodic, limited):
    V = np.vstack(to_prim(U))                      # 3 x (N+4)
    sl = np.zeros_like(V)
    for k in range(3):
        sl[k, 1:-1] = minmod_slope(V[k], limited)
    VL = V[:, 1:-1]+0.5*sl[:, 1:-1]                # cara derecha de la celda i
    VR = V[:, 2:]-0.5*sl[:, 2:]                    # cara izquierda de la celda i+1
    F = hllc_flux(VL[0], VL[1], VL[2], VR[0], VR[1], VR[2])
    R = np.zeros_like(U)
    R[:, 2:-2] = -(F[:, 1:-1]-F[:, :-2])/h
    return R
 
def fill_ghost(U, periodic):
    if periodic:
        U[:, :2] = U[:, -4:-2]; U[:, -2:] = U[:, 2:4]
    else:
        U[:, 0] = U[:, 2]; U[:, 1] = U[:, 2]
        U[:, -1] = U[:, -3]; U[:, -2] = U[:, -3]
 
def evolve(U, h, tend, periodic, limited, cfl=0.4):
    t = 0.0
    while t < tend-1e-14:
        fill_ghost(U, periodic)
        r, u, p = to_prim(U[:, 2:-2])
        dt = cfl*h/np.max(np.abs(u)+np.sqrt(G*p/r))
        dt = min(dt, tend-t)
        U1 = U+dt*residual(U, h, periodic, limited); fill_ghost(U1, periodic)
        U = 0.5*(U+U1+dt*residual(U1, h, periodic, limited))
        t += dt
    fill_ghost(U, periodic)
    return U
 
def order_table(ns, errs):
    return [None]+[np.log(errs[i-1]/errs[i])/np.log(ns[i]/ns[i-1]) for i in range(1, len(ns))]
 
def run_sod(ns, tend=0.2):
    e1, ei = [], []
    for N in ns:
        h = 1.0/N
        x = (np.arange(N)+0.5)*h
        r = np.where(x < 0.5, 1.0, 0.125); u = np.zeros(N)
        p = np.where(x < 0.5, 1.0, 0.1)
        U = np.zeros((3, N+4))
        U[:, 2:-2] = np.vstack([r, r*u, p/(G-1)+0.5*r*u*u])
        U = evolve(U, h, tend, False, True)
        re = sample_exact_riemann(x-0.5, tend, 1.0, 0.0, 1.0, 0.125, 0.0, 0.1)
        d = np.abs(U[0, 2:-2]-re)
        e1.append(np.sum(d)*h); ei.append(np.max(d))
    return e1, ei
 
def run_wave(ns, limited, tend=1.0):
    e1 = []
    for N in ns:
        h = 1.0/N
        x = (np.arange(N)+0.5)*h
        prof = lambda z: 1.0+0.5*np.exp(-60.0*(np.mod(z-0.5, 1.0)-0.5)**2)
        r = prof(x); u = np.ones(N); p = np.ones(N)
        U = np.zeros((3, N+4))
        U[:, 2:-2] = np.vstack([r, r*u, p/(G-1)+0.5*r*u*u])
        U = evolve(U, h, tend, True, limited)
        e1.append(np.sum(np.abs(U[0, 2:-2]-prof(x-tend)))*h)
    return e1
 
 
def rankine_hugoniot(m1):
    """cocientes de salto del choque estacionario (caso choque-vortice)"""
    rr = ((G+1)*m1**2)/((G-1)*m1**2+2)
    pr = (2*G*m1**2-(G-1))/(G+1)
    m2 = np.sqrt((1+0.5*(G-1)*m1**2)/(G*m1**2-0.5*(G-1)))
    return rr, pr, m2
 
NS = [100, 200, 400, 800]
s1, si = run_sod(NS)
wl, wc = run_wave(NS, True), run_wave(NS, False)
o_s1, o_si, o_wl, o_wc = (order_table(NS, e) for e in (s1, si, wl, wc))
fmt = lambda o: "    -  " if o is None else f"{o:7.3f}"
print("Sod shock tube (two discontinuities), t = 0.2")
print("    N      L1 error   order     Linf error  order")
for i, N in enumerate(NS):
    print(f" {N:4d}   {s1[i]:.4e} {fmt(o_s1[i])}   {si[i]:.4e} {fmt(o_si[i])}")
print()
print("smooth density wave (exact = translated Gaussian), t = 1.0")
print("    N      L1 minmod  order     L1 unlim.   order")
for i, N in enumerate(NS):
    print(f" {N:4d}   {wl[i]:.4e} {fmt(o_wl[i])}   {wc[i]:.4e} {fmt(o_wc[i])}")
ps, us = riemann_star_pressure(1.0, 0.0, 1.0, 0.125, 0.0, 0.1)
print(f"\nSod star state   p* = {ps:.6f}  u* = {us:.6f}")
rr, pr, m2 = rankine_hugoniot(1.2)
print(f"RH at M1 = 1.2   rho2/rho1 = {rr:.4f}  p2/p1 = {pr:.4f}  M2 = {m2:.4f}")
amp = (G-1)*25.0/(8*G*np.pi**2)
for r in (2.5, 5.0):
    print(f"vortex tail r = {r:.1f}  exp(1-r^2) = {np.exp(1-r*r):.3e}  dT(beta=5) = {amp*np.exp(1-r*r):.3e}")

Salida:

Sod shock tube (two discontinuities), t = 0.2
    N      L1 error   order     Linf error  order
  100   6.5089e-03     -     6.7837e-02     -  
  200   3.7579e-03   0.792   8.3037e-02  -0.292
  400   2.1461e-03   0.808   8.3428e-02  -0.007
  800   1.2420e-03   0.789   8.5667e-02  -0.038
 
smooth density wave (exact = translated Gaussian), t = 1.0
    N      L1 minmod  order     L1 unlim.   order
  100   5.6574e-03     -     2.0251e-03     -  
  200   2.0125e-03   1.491   5.0589e-04   2.001
  400   5.7467e-04   1.808   1.2630e-04   2.002
  800   1.5736e-04   1.869   3.1562e-05   2.001
 
Sod star state   p* = 0.303130  u* = 0.927453
RH at M1 = 1.2   rho2/rho1 = 1.3416  p2/p1 = 1.5133  M2 = 0.8422
vortex tail r = 2.5  exp(1-r^2) = 5.248e-03  dT(beta=5) = 4.747e-04
vortex tail r = 5.0  exp(1-r^2) = 3.775e-11  dT(beta=5) = 3.415e-12

Al subir sharpness $m$ , la campana pasa de gaussiana a casi rectangular y el orden medido baja de cerca de 2 hacia 1.

sharpness m = 2 N = 200 smear = 3 cells 
m = 2 gives a measured L1 order of 2.00 — the scheme never changed.

La línea sin limitador, la que dio 2.001, se desliza hacia 1 junto con $m$ . Con show Linf el mismo cálculo recorre de 2 a 0. La línea del esquema no cambia.

Los valores $p^{*} = 0.303130$ y $u^{*} = 0.927453$ son el estado estrella de Sod de los libros de texto, lo que sitúa la velocidad del choque en 1.7522. Si el código reproduce esos números es lo que realmente está mirando el tubo de choque de Sod.

El papel que cumple el vórtice isentrópico#

El caso estándar para medir orden en dos dimensiones es la convección del vórtice isentrópico. Se tiende una corriente uniforme sobre un dominio periódico $[0,10]^2$ y se suma encima una perturbación de vórtice.

\delta u = -\frac{\beta}{2\pi}(y - y_c)\,e^{(1-r^2)/2}, \qquad \delta v = \frac{\beta}{2\pi}(x - x_c)\,e^{(1-r^2)/2}

Aquí $\beta$ es la intensidad del vórtice, $r$ la distancia al centro y $(x_c, y_c)$ las coordenadas del centro. La temperatura sale de la relación isentrópica.

T = 1 - \frac{(\gamma-1)\beta^2}{8\gamma\pi^2}e^{1-r^2}, \qquad \rho = T^{1/(\gamma-1)}, \qquad p = \rho^{\gamma}

Este campo es una solución exacta que se traslada a la velocidad de la corriente libre. Es suave y no tiene discontinuidades, así que medir orden sobre él es legítimo.

La trampa está en las colas. La perturbación decae como $e^{1-r^2}$ , pero nunca llega a cero. Si se recorta en la frontera del dominio, la solución "exacta" deja de ser exacta.

En $[0,10]^2$ la distancia a la frontera es 5. Allí $e^{1-r^2}$ vale $3.8 \times 10^{-11}$ y, con $\beta = 5$ , la perturbación de temperatura es $3.4 \times 10^{-12}$ . Eso es holgadamente seguro para un esquema de segundo orden. Al reducir el dominio a $[0,5]^2$ , en $r = 2.5$ esas mismas cantidades pasan a $5.2 \times 10^{-3}$ y $4.7 \times 10^{-4}$ . Un esquema de alto orden choca con ese piso mientras todavía está refinando, y ahí el orden se derrumba.

Choque–vórtice no mide orden#

El tercer caso es de otra naturaleza. Un choque estacionario de $M = 1.2$ se coloca a la izquierda, el estado derecho sale de las relaciones de Rankine–Hugoniot y un vórtice isentrópico se empuja a través del choque.

Este problema no tiene solución cerrada, así que no se puede extraer ningún orden.

Lo que sí se puede cotejar es el salto. En $M_1 = 1.2$ y $\gamma = 1.4$ los cocientes son $\rho_2/\rho_1 = 1.3416$ , $p_2/p_1 = 1.5133$ , $M_2 = 0.8422$ . Si el código no reproduce esto, toda imagen aguas abajo del choque pierde sentido. Las integrales de masa y momento en la dirección periódica deberían mantenerse fijas a precisión de máquina. Un funcional, como la amplitud acústica o una integral de vorticidad, se mide en tres mallas y la más fina hace de referencia. Eso es convergencia de malla, no un orden de precisión.

Lo que se elige al escribir un número en un informe#

La frase "se confirmó precisión de segundo orden" esconde tres decisiones: qué solución, qué norma y si el limitador estaba encendido.

Caso	Solución exacta	Suavidad	Qué puede medir
Tubo de Sod	sí (Riemann)	discontinua	posición y velocidad del choque, monotonía, positividad
Vórtice isentrópico	sí (traslación)	suave	orden de convergencia de diseño
Choque–vórtice	no	mixta	conservación, consistencia RH, convergencia de malla

Si Sod reporta 0.79, lo que hay que corregir es la frase del informe, no el esquema. Si el vórtice reporta 2.00, la configuración del limitador va al lado. Medir el orden sobre una solución exacta suave es el mismo trabajo que contrastar un tubo con el perfil de Poiseuille. Elegir dónde medir es la mitad del trabajo.