Cuando Lax–Wendroff tembló en un escalón — Un esquema de advección TVD de 2° orden con flux limiters

Dentro de la región de Sweby el escalón no crea picos nuevos. En una onda suave el orden sigue partiéndose según el limitador, de cerca de 1.3 a 2.0.

En un caso de prueba aparecieron jorobas extrañas delante y detrás del escalón. Estaba resolviendo la ecuación de advección lineal (linear advection, pura traslación con un campo de velocidades constante) con Lax–Wendroff, y la condición inicial era un simple pulso rectangular. Al ser un esquema de 2° orden debía resultar más nítido que el upwind de 1° orden, pero el resultado mostraba undershoot delante del escalón y overshoot detrás. El registro decía CFL = 0.5, dentro del límite de estabilidad.

El problema no era inestabilidad. El algoritmo no era monótono.

Los botones del limitador arrancan en lax-wendroff. Con el deslizador de CFL en 0.5, la curva naranja forma jorobas alrededor del escalón. Al pasar a superbee, esas jorobas se retiran.

CFL0.50steps: 0

La naranja es la solución numérica. La línea punteada cian es la exacta. minmod queda romo a ambos lados, sin joroba. superbee puede verse más afilado que el escalón original. Cerca de CFL 0.9 hasta los limitadores quietos se extienden más. La estabilidad no cancela la difusión numérica.

El origen de las oscilaciones — la pared de Godunov#

La ecuación de advección lineal 1D es

\partial_t q + u\,\partial_x q = 0

Con $u$ constante, la solución es simplemente el perfil inicial trasladado hacia la derecha. Numéricamente se define un flujo $f_{i-1/2}$ en cada cara de celda y se escribe

q_i^{n+1} = q_i^n - \frac{\Delta t}{\Delta x}\bigl(f_{i+1/2} - f_{i-1/2}\bigr)

Esta es la flux-conserving form, que preserva magnitudes conservadas (masa, momento, energía) hasta la precisión de máquina. Lax–Wendroff es un esquema lineal que añade una corrección Taylor de 2° orden a esta forma.

Pero un teorema probado por Godunov en 1959 es implacable: no existe ningún esquema lineal para la advección lineal de coeficientes constantes que sea simultáneamente monótono y de 2° orden o superior. Hay que renunciar a uno de los dos. Lax–Wendroff eligió el 2° orden y sacrificó la monotonía — de ahí las oscilaciones en el escalón.

La condición TVD: desigualdad de Harten#

Harten (1983) ofreció una salida: un esquema cuya "variación total disminuye".

\mathrm{TV}(q) = \sum_i |q_{i+1} - q_i|, \qquad \mathrm{TV}(q^{n+1}) \le \mathrm{TV}(q^n)

El lado izquierdo mide el total de ondulaciones de la solución numérica. Si este valor nunca crece, no pueden aparecer nuevos extremos (overshoot, undershoot). Los esquemas que satisfacen esto se llaman TVD.

Para rodear la pared de Godunov, el esquema debe volverse no lineal — su forma debe cambiar según la suavidad local de la solución. 2° orden en regiones suaves, upwind de 1° orden cerca de gradientes pronunciados. Ese conmutador es el flux limiter.

Flux Limiter en una sola línea#

El flujo de Lax–Wendroff equivale al upwind de 1° orden más una corrección:

f_{i-1/2} = u\,q_{i-1} + \tfrac{1}{2}\,|u|\!\left(1 - \left|\tfrac{u\Delta t}{\Delta x}\right|\right) \varphi(r)\,(q_i - q_{i-1})

donde $u$ es la velocidad de advección, $\Delta x$ el ancho de celda, $\varphi(r)$ el limitador y $r$ la razón entre pendientes vecinas:

r_{i-1/2} = \frac{q_{i-1} - q_{i-2}}{q_i - q_{i-1}}

Con una solución suave $r \approx 1$ ; cerca de los extremos $r < 0$ o $r \gg 1$ . El limitador observa ese $r$ y ajusta la intensidad de la corrección.

Casos particulares:

$\varphi(r) = 0$ : donor-cell (upwind de 1° orden)
$\varphi(r) = 1$ : Lax–Wendroff
$\varphi(r) = r$ : Beam–Warming

En el diagrama de Sweby, cualquier limitador que satisfaga $0 \le \varphi(r) \le \min(2r,\,2)$ es TVD. Cuatro limitadores representativos dentro de ese trapecio. Apagar con $\varphi(r) = 0$ deja el upwind de primer orden. Lo que ese upwind embota, a cambio de matar la oscilación, está en el intercambio cuando FTCS empieza a temblar.

Limitador	Definición	Carácter
minmod	$\max(0,\min(1,r))$	El más cauto. Conserva escalones; en zonas suaves es algo torpe
superbee	$\max(0,\min(2r,1),\min(r,2))$	El más afilado. Puede convertir ondas suaves en escaleras
van Leer	$(r+	r
MC	$\max(0,\min(\tfrac{1+r}{2},2,2r))$	Opción por defecto muy sólida

Comparación directa en Python#

Con vectorización basada en roll, cabe en menos de 40 líneas. 200 celdas, $\mathrm{CFL} = 0.5$ , 500 pasos (condición de contorno periódica; la onda recorre el dominio unas 2.5 veces).

import numpy as np
 
def limiter_phi(name, r):
    r = np.where(np.isfinite(r), r, 0.0)
    if name == 'donor':    return np.zeros_like(r)
    if name == 'lw':       return np.ones_like(r)
    if name == 'minmod':   return np.maximum(0, np.minimum(1, r))
    if name == 'superbee':
        return np.maximum.reduce([np.zeros_like(r),
                                  np.minimum(2*r, 1),
                                  np.minimum(r, 2)])
    if name == 'vanleer':  return (r + np.abs(r)) / (1 + np.abs(r) + 1e-30)
    if name == 'mc':
        return np.maximum(0, np.minimum.reduce([0.5*(1+r),
                                                2*np.ones_like(r),
                                                2*r]))
 
def advect_limited_fvm(q, cfl, name):
    """Advección 1D en volúmenes finitos, periódica, u = +1."""
    dq  = q - np.roll(q, 1)                  # cara i-1/2: q_i - q_{i-1}
    r   = np.roll(dq, 1) / (dq + 1e-30)      # (q_{i-1}-q_{i-2}) / (q_i-q_{i-1})
    phi = limiter_phi(name, r)
    flx = np.roll(q, 1) + 0.5 * phi * (1 - cfl) * dq
    return q - cfl * (np.roll(flx, -1) - flx)
 
N, cfl, steps = 200, 0.5, 500
x  = np.arange(N) / N
q0 = ((x > 0.2) & (x < 0.4)).astype(float)
for name in ['lw', 'minmod', 'superbee', 'vanleer', 'mc']:
    q = q0.copy()
    for _ in range(steps):
        q = advect_limited_fvm(q, cfl, name)
    l1 = np.abs(q - q0).sum() / N
    print(f'{name:<9s}  L1={l1:.4f}  overshoot={q.max()-1:.3f}  undershoot={q.min():.3f}')

Resultados típicos:

Limitador	Error $L_1$	Overshoot	Undershoot
lax-wendroff	0.056	+0.181	−0.133
minmod	0.081	0.000	0.000
superbee	0.042	+0.002	0.000
van Leer	0.061	0.000	0.000
MC	0.051	0.000	0.000

Solo lax-wendroff muestra un overshoot cercano al 20%; los demás quedan en cero. El ranking de error también llama la atención: superbee tiene el menor $L_1$ . Eso significa que mantiene el escalón bien definido, a cambio de distorsionar ondas seno suaves en forma de escaleras.

Si se repite la prueba con una onda suave ( $q_0 = \sin(2\pi x)$ , 500 pasos), el orden de convergencia baja a unos 1.5 para minmod, 2.0 para van Leer/MC y 1.3 para superbee. La etiqueta "2° orden" se desliza cada vez que el limitador se activa cerca de un extremo y cae localmente al 1° orden.

En r=1 los cuatro limitadores devuelven la misma phi#

Al meter $r=1$ en la tabla, minmod, superbee, van Leer y MC dan 1. $\varphi=1$ es Lax–Wendroff. En el centro suave aplican la misma corrección. Se separan solo cuando $r$ se aleja de 1.

Con $r=2$ , minmod se detiene en 1. superbee se pega a 2, el borde superior de Sweby. van Leer vale $4/3$ . MC vale 1.5. Ese borde es la razón de que superbee deje el $L_1$ del escalón en 0.042, el menor de la tabla. La misma elección baja el orden de la onda seno a cerca de 1.3. van Leer y MC se quedan cerca de 2.0. minmod cae cerca de 1.5. En un extremo con $r \le 0$ , los cuatro bajan a 0 y pasan a donor-cell.

El coeficiente de la corrección ya está fijo en el código. Con velocidad 1 y $\mathrm{CFL}=0.5$ , el factor delante del salto es 0.25. Incluso con $\varphi=1$ solo se suma un cuarto del salto vecino. Quinientos pasos recorren $0.5\times 500=250$ celdas. En una malla de 200 eso es 1.25 vueltas del contorno periódico. El overshoot +0.181 es el 18 por ciento de un pulso de altura 1. El undershoot −0.133 es lo que el suelo baja bajo cero. El $L_1$ de minmod, 0.081, es unas 1.9 veces el 0.042 de superbee. Borrar la joroba cuesta ese emborronado extra.

Checklist práctico#

Para captura de choques, el punto de partida es minmod o van Leer. La prioridad es cero oscilación.
Para LES y turbulencia a gran escala, superbee es riesgoso. Convierte estructuras suaves en escaleras y distorsiona el espectro de energía. van Leer o MC son el compromiso.
Bajar a upwind de 1° orden en celdas de frontera por falta de stencil reduce la precisión global a 1° orden. Conviene conservar el stencil del limitador en la frontera, o bien usar ghost cells.
El limitedLinear 1 de OpenFOAM es un limitador de un parámetro de la familia van Leer. Coeficiente 0 = upwind de 1° orden; 1.0 = región TVD completa.
El orden de convergencia se verifica siempre con una solución suave. Medir pendientes $L_1$ sobre un caso con discontinuidad hace que la discontinuidad (no el limitador) marque el orden.

superbee en un escalón baja la onda seno al orden 1.3#

En la tabla el $L_1$ de superbee es 0.042, el menor. El overshoot es +0.002, casi invisible sobre el escalón. El mismo limitador sobre una onda seno baja el orden a cerca de 1.3. van Leer y MC se quedan cerca de 2.0. La elección que protege el salto gasta el orden de la onda suave.

Si el escalar es una densidad, el undershoot −0.133 de lax-wendroff es la puerta a una densidad negativa. El limitador que cierra esa puerta está en donde se escapa un valor positivo. Para un escalón pasivo bastan minmod o van Leer. Para una densidad, se abre la puerta siguiente.