CFL = 1 tarda diez mil pasos, CFL = 100 explota — Switched Evolution Relaxation

El residual estacionario es quisquilloso con el CFL. En 1, diez mil iteraciones siguen arriba. En 100, el primer paso revienta la solución. SER (Switched Evolution Relaxation) deja que el residual elija el paso entre esos dos extremos. Si baja, el paso crece. Si salta, el paso se encoge.

Al mover el control de CFL fijo de 1 a 20, el residual ámbar salta a diverged, mientras SER solo expone el CFL inicial y el exponente.

Con $\mathrm{CFL}_0=1$ y $p=1$ la curva cian se quiebra hacia abajo al final. Subir $p$ a 1.5 la hace más empinada, pero en cuanto el residual inicial oscila, la curva cian oscila con él.

En un solver estacionario el número de CFL (el paso de tiempo adimensionalizado por el tamaño de malla y la velocidad de onda) ya no trata de precisión. Es una sola perilla que sostiene a la vez la velocidad de convergencia y la estabilidad.

El estado estacionario se persigue con tiempo falso#

La ecuación estacionaria es una línea: $R(U) = 0$. Aquí $R$ reúne convección, difusión y fuentes en un residual espacial. Resolverla directo implica romper un sistema algebraico no lineal. En su lugar, casi todos los códigos de CFD añaden un término de tiempo falso.

Aquí $\tau$ es seudotiempo (pseudo-time) sin significado físico. Cuando $\tau \to \infty$, $\partial U / \partial \tau \to 0$, así que $R(U)=0$: justo la solución estacionaria buscada. Se descarta la precisión temporal y solo se apunta al estado estacionario. Esto es la continuación seudotransitoria (pseudo-transient continuation).

Al avanzar el tiempo falso de forma implícita (backward Euler) y linealizar el residual, cada iteración queda así:

Aquí $\partial R/\partial U$ es el Jacobiano de flujo y $\Delta U = U^{n+1}-U^n$. Cada iteración resuelve una vez la matriz del lado izquierdo para obtener $\Delta U$.

Empezar pequeño, aterrizar en Newton#

Esa única línea contiene toda la intuición de SER. Llevemos $\Delta\tau$ a ambos extremos.

Cuando $\Delta\tau \to 0$, el término $I/\Delta\tau$ domina la matriz y $\Delta U \approx -\Delta\tau\, R$. Eso es avance explícito puro. Robustísimo pero lento como una tortuga.

Cuando $\Delta\tau \to \infty$, el término $I/\Delta\tau$ desaparece y solo queda $(\partial R/\partial U)\,\Delta U = -R$. Eso es exactamente el método de Newton, con convergencia cuadrática (el error se encoge por su cuadrado) cerca de la solución. Pero una mala estimación inicial lo hace divergir.

Así que la estrategia ideal es clara. Mantener $\Delta\tau$ pequeño al principio para sobrevivir a la no linealidad y luego agrandarlo cuando el residual cae y se está cerca de la solución, para absorber la convergencia cuadrática de Newton. El detalle es "cuándo y cuánto". SER cede ese juicio al residual.

SER — agrandar el CFL con el inverso del residual#

La regla que Mulder y van Leer propusieron en 1985 es asombrosamente simple. Agrandar el CFL actual por la razón entre el residual inicial y el actual.

$\mathrm{CFL}_0$ es el CFL inicial, $\lVert R_n \rVert$ es la norma L2 del $n$-ésimo residual, $p$ es un exponente normalmente cercano a 1 y $\mathrm{CFL}_{\max}$ pone tope al desbocamiento. Si el residual cae diez veces (con $p=1$), el CFL crece diez veces. El paso de seudotiempo se convierte a partir de este CFL y la velocidad de onda local.

El denominador suma la velocidad de onda convectiva $|u_i|$ y la difusiva $2\nu/\Delta x$. Puede diferir por celda, así que se combina de modo natural con el avance temporal local (local time-stepping).

Con N=80 el espaciado es 1/81. Con $\nu=0.02$, $2\nu/\Delta x$ vale 3.24, más que $|u_i|$ en la estimación lineal. Un CFL fijo de 3 deja $\Delta\tau$ por debajo de 0.01. El tope 1e5 lleva $\Delta\tau$ a los cientos, así que $I/\Delta\tau$ casi desaparece. El código reparte el máximo de $|u|$ entre todas las celdas, también en el choque, donde la velocidad ronda 0. Es más estricto que el paso local de la fórmula.

Python — resolver Burgers estacionario de dos maneras#

Resolvemos la ecuación de Burgers viscosa estacionaria $u u_x = \nu u_{xx}$. Las condiciones de borde $u(0)=1$, $u(1)=-1$ crean una capa de choque estacionaria en el centro. Un paso implícito de seudotiempo resuelve una matriz tridiagonal con el algoritmo de Thomas. La misma sustitución hacia atrás que vuelve incondicionalmente estable a la difusión 1D es donde este lado izquierdo suma el término temporal.

import numpy as np
 
def steady_residual(u, nu, dx):
    """R(u) = u u_x - nu u_xx, bordes u(0)=1, u(1)=-1."""
    um = np.concatenate(([1.0], u[:-1]))   # vecino izquierdo
    up = np.concatenate((u[1:], [-1.0]))   # vecino derecho
    conv = u * (up - um) / (2 * dx)
    diff = nu * (up - 2 * u + um) / dx**2
    return conv - diff
 
def thomas_solve(a, b, c, d):
    """Resuelve un sistema tridiagonal (sub/main/super diag a,b,c, rhs d)."""
    n = len(d)
    b, d = b.copy(), d.copy()
    for i in range(1, n):
        m = a[i] / b[i - 1]
        b[i] -= m * c[i - 1]
        d[i] -= m * d[i - 1]
    x = np.empty(n)
    x[-1] = d[-1] / b[-1]
    for i in range(n - 2, -1, -1):
        x[i] = (d[i] - c[i] * x[i + 1]) / b[i]
    return x
 
def ptc_update(u, nu, dx, dt):
    """Un paso implícito de seudotiempo: (I/dt + dR/du) du = -R."""
    n = len(u)
    um = np.concatenate(([1.0], u[:-1]))
    up = np.concatenate((u[1:], [-1.0]))
    a = -u / (2 * dx) - nu / dx**2                       # subdiagonal
    b = (up - um) / (2 * dx) + 2 * nu / dx**2 + 1.0 / dt  # diagonal principal
    c = u / (2 * dx) - nu / dx**2                        # superdiagonal
    du = thomas_solve(a, b, c, -steady_residual(u, nu, dx))
    return u + du
 
def ser_schedule(r0, rn, cfl0, p, cfl_max):
    return min(cfl0 * (r0 / rn) ** p, cfl_max)
 
def converge(nu=0.02, ser=True, cfl0=1.0, p=1.0, cfl_max=1e5,
             N=80, tol=1e-9, max_iter=600):
    dx = 1.0 / (N + 1)
    x = np.linspace(dx, 1 - dx, N)
    u = 1 - 2 * x                          # estimacion inicial lineal
    r0 = np.linalg.norm(steady_residual(u, nu, dx)) / np.sqrt(N)
    rn = r0
    for it in range(1, max_iter + 1):
        cfl = ser_schedule(r0, rn, cfl0, p, cfl_max) if ser else cfl0
        dt = cfl * dx / (np.abs(u).max() + 2 * nu / dx)
        u = ptc_update(u, nu, dx, dt)
        rn = np.linalg.norm(steady_residual(u, nu, dx)) / np.sqrt(N)
        if rn / r0 < tol:
            return u, it
    return u, max_iter
 
for tag, kw in [("fixed CFL=3", dict(ser=False, cfl0=3.0)),
                ("SER  p=1.0", dict(ser=True, cfl0=1.0, p=1.0))]:
    _, iters = converge(**kw)
    print(f"{tag:14s} -> {iters:4d} iters")

La salida se ve así.

fixed CFL=3    ->  588 iters
SER  p=1.0     ->   34 iters

Al bajar el control del tope, el cronograma muestra dónde las 34 iteraciones impresas vuelven a alargarse. Cian es el CFL, rosa es el residual.

Mientras el residual está plano al inicio, el CFL se queda cerca de $\mathrm{CFL}_0$. En el instante en que el residual empieza a colapsar, el CFL se dispara casi en vertical. Esa es la esencia de SER: cuanto mejor se vuelve la solución, mayor es el paso que permite. Con el tope en 20 el CFL choca con el techo, pierde la aceleración tipo Newton del final y el conteo de iteraciones vuelve a subir.

Para alcanzar una solución estacionaria de la misma precisión, el CFL fijo necesita cientos de iteraciones y SER apenas decenas. La brecha se abre en la fase final, porque SER lleva el CFL a los miles y de hecho se transforma en el método de Newton.

CFL fijo y SER, lado a lado#

588 partido por 34 da cerca de 17. La tolerancia es una razón de residuales bajo 1e-9. Con $p=1$, CFL inicial 1 y tope 1e5, el techo llega cuando la razón alcanza 1e5, es decir 1e-5. Quedan cuatro décadas hasta 1e-9, todas con CFL 1e5, ya del lado de Newton. El CFL fijo de 3 gasta las 588 iteraciones sin ese tramo.

El cronograma cabe en una línea. $\mathrm{CFL}_n = \min(\mathrm{CFL}_0 (\lVert R_0\rVert/\lVert R_n\rVert)^p,\ \mathrm{CFL}_{\max})$. Si el residual baja, el paso crece, y se detiene en el tope.

Trampas al encender SER en producción#

Conviene fijar siempre $\mathrm{CFL}_{\max}$. Si se deja en infinito, un solo pico del residual desboca el CFL y mata la solución.

Cuando el residual oscila de forma no monótona (non-monotone), $R_0/R_n$ puede caer por debajo de 1 y encoger el CFL. No es un error: es una válvula de seguridad. Si el residual vuelve a empeorar, acortar el paso es lo correcto.

$p$ no es aceleración gratis. Un $p>1$ agranda el CFL enormemente ante una caída pequeña del residual, elevando el riesgo. Un arranque robusto es $p=1$, $\mathrm{CFL}_0 \approx 1$.

Si el Jacobiano es inexacto (Jacobiano aproximado o sin matriz), la convergencia cuadrática de Newton cae a primer orden con CFL grande. Entonces conviene bajar $\mathrm{CFL}_{\max}$ y comprar estabilidad.

Un tope vacío mata antes al Newton tardío#

Dejar $\mathrm{CFL}_{\max}$ vacío hace que un solo pico del residual empuje el CFL tras la razón de residuales sin techo. La iteración implícita de seudotiempo se vuelve el método de Newton cuando $\Delta\tau\to\infty$. Ese orden cuadrático aparece solo cerca de la solución. Saltar ahí al principio hace que la matriz empuje la solución al revés. El tope 1e5 del código es una puerta: el paso entra en esa zona solo después de que el residual haya caído cien mil veces.

La misma puerta hace falta cuando el Jacobiano no se puede formar con exactitud en cada iteración. Si el orden cuadrático del CFL grande ya bajó a primer orden, la opción que queda es un tope más bajo. Aproximar esa matriz con barridos, sin ensamblarla, es el tema de LU-SGS después de que un paso explícito revienta. El paso implícito se queda, y la matriz no tiene que almacenarse.