Cuando el método explícito revienta, resolver implícito sin invertir la matriz — LU-SGS

Estabilidad del avance temporal implícito y la factorización LU-SGS, verificadas con código

La simulación divergió a las 3 de la madrugada. La última línea del registro siempre decía lo mismo: pressure residual = NaN. Bajé el CFL de 0.9 a 0.5, y luego a 0.3. Cada paso avanzaba tan poco tiempo que llegar al estado estacionario tomaría días. Ese muro le resulta familiar a cualquiera que haga cálculo estacionario compresible. Este artículo muestra el camino para superarlo —por qué el avance temporal implícito permite un paso de tiempo grande— y luego implementa LU-SGS, un método que resuelve ese sistema sin invertir nunca la matriz completa. Al final tendrás en la mano un solver de Burgers que sobrevive a un CFL de 5.

La región de estabilidad se abre al infinito#

La estabilidad del avance temporal se decide con una sola ecuación modelo.

\frac{du}{dt} = \lambda u

Aquí $\lambda$ es el autovalor que produce la discretización espacial (viene de las velocidades de onda y de la viscosidad) y $u$ es la amplitud de un modo. Tras un paso, la amplitud queda multiplicada por el factor de amplificación $G$ . Hace falta $|G|\le 1$ para evitar la explosión.

El avance explícito (Euler hacia adelante) da

G = 1 + \lambda\Delta t

de modo que $|1+\lambda\Delta t|\le 1$ : el punto $\lambda\Delta t$ debe quedar dentro de un disco pequeño centrado en $-1$ y de radio $1$ . Una malla viscosa o una celda delgada agranda $\lambda$ , y este sale rápido de ese disco.

El avance implícito (Euler hacia atrás) es lo contrario.

G = \frac{1}{1 - \lambda\Delta t}

Para todo modo físico amortiguado con $\mathrm{Re}(\lambda)<0$ , $|G|<1$ . Permanece estable por grande que sea $\Delta t$ . La restricción de estabilidad se disuelve y el paso de tiempo solo necesita ser tan pequeño como exija la precisión.

Juega con los parámetros en la simulación de abajo.

Re(λΔt) = -1.00Im(λΔt) = 1.40

Cyan = stable region (|G| ≤ 1). Forward Euler and RK4 are bounded blobs — push λΔt past their edge and the mode blows up. Backward Euler and Trapezoidal cover the entire left half plane, so any decaying physical mode stays stable no matter how large Δt grows.

Elige Euler hacia adelante y arrastra $\mathrm{Re}(\lambda\Delta t)$ hasta $-2.5$ : el punto sale de la región cian y cambia a UNSTABLE. Cambia a Euler hacia atrás y todo el semiplano izquierdo se llena de cian, así que el mismo punto vuelve a ser estable.

Invertir la matriz completa es un costo impagable#

Nada es gratis. El método implícito pide el residuo en el estado futuro.

\frac{U_i^{n+1} - U_i^n}{\Delta t} = -R_i(U^{n+1})

$R_i$ es la divergencia de flujo (residuo) de la celda $i$ y $U$ es la variable conservada. Al linealizar $R(U^{n+1})$ alrededor del estado actual,

R_i(U^{n+1}) \approx R_i^n + \frac{\partial R_i}{\partial U}\,\Delta U

queda un sistema lineal para $\Delta U = U^{n+1}-U^n$ .

\left( \frac{I}{\Delta t} + \frac{\partial R}{\partial U} \right) \Delta U = -R^n

El problema es la matriz $A$ del lado izquierdo. Con cincuenta millones de celdas, $A$ tiene cincuenta millones de dimensiones. La inversa directa es imposible en memoria y en tiempo. Un solver de Krylov como GMRES también sirve, pero hay que almacenar la matriz de forma explícita y necesita un precondicionador. En los años noventa, los analistas de flujo compresible encontraron un camino más barato: aproximar la matriz sin almacenamiento, sin multiplicaciones y con dos barridos.

LU-SGS — separar en D·L·U y barrer dos veces#

Descompón $A$ en un triángulo inferior por bloques $L$ , una diagonal por bloques $D$ y un triángulo superior por bloques $U$ .

A = D + L + U

LU-SGS (Gauss-Seidel simétrico inferior-superior) aproxima esta $A$ con el producto

(D + L)\,D^{-1}\,(D + U)\,\Delta U = -R^n

que descarta el término $L D^{-1} U$ creado por el producto $LU$ . Cuando $\Delta U$ es pequeño, es despreciable. La resolución se reduce a dos sustituciones triangulares en secuencia.

Barrido hacia adelante — triángulo inferior, de adelante hacia atrás:

\Delta U_i^{*} = D_i^{-1}\!\left( -R_i^n - L_i\,\Delta U_{i-1}^{*} \right)

Barrido hacia atrás — triángulo superior, de atrás hacia adelante:

\Delta U_i = \Delta U_i^{*} - D_i^{-1}\,U_i\,\Delta U_{i+1}

La clave es con qué se llenan $D$ , $L$ y $U$ . Yoon y Jameson (1988) aproximaron el jacobiano no lineal del flujo con su radio espectral $\rho$ . Al separar el upwind de primer orden como $A^{\pm} = \tfrac{1}{2}(a \pm |a|)$ , la diagonal queda $D_i = I/\Delta t + |a_i|/\Delta x$ y las no diagonales solo cargan los jacobianos de cara separados. Nunca construyes la matriz entera: bastan unos pocos escalares (o un bloque pequeño) por celda.

Python — avanzar Burgers implícito con una sola factorización#

El problema de juguete es la formación de un choque en la ecuación de Burgers 1D. Una onda senoidal suave se empina hasta un choque con el tiempo. El método explícito fijaría el CFL por debajo de 1. El paso implícito LU-SGS sobrevive a CFL 5.

import numpy as np
 
def upwind_flux_residual(u, dx):
    """Flujo upwind de primer orden para el residuo R_i de d(u^2/2)/dx (periodico)."""
    f = 0.5 * u * u
    a = u                       # velocidad de onda local df/du
    f_left  = np.where(a >= 0, np.roll(f, 1), f)
    a_right = np.roll(a, -1)
    f_right = np.where(a_right >= 0, f, np.roll(f, -1))
    return (f_right - f_left) / dx
 
def lu_sgs_burgers_step(u, dt, dx):
    """Un paso de tiempo: avanzar el sistema implicito con una pasada LU-SGS."""
    n = u.size
    a  = u
    ap = 0.5 * (a + np.abs(a))         # A+ (contribucion del triangulo inferior)
    am = 0.5 * (a - np.abs(a))         # A- (contribucion del triangulo superior)
    D  = 1.0 / dt + np.abs(a) / dx     # bloque diagonal
    R  = upwind_flux_residual(u, dx)
 
    dus = np.zeros(n)                  # barrido adelante: (D+L) dU* = -R
    for i in range(n):
        lower = ap[i - 1] / dx * dus[i - 1]
        dus[i] = (-R[i] + lower) / D[i]
 
    du = dus.copy()                   # barrido atras: (D+U) dU = D dU*
    for i in range(n - 1, -1, -1):
        upper = am[(i + 1) % n] / dx * du[(i + 1) % n]
        du[i] = dus[i] - upper / D[i]
    return u + du
 
def drive_burgers_shock(n=200, cfl=5.0, tmax=0.22):
    dx = 1.0 / n
    x  = (np.arange(n) + 0.5) * dx
    u  = 0.5 + 0.5 * np.sin(2 * np.pi * x)   # condicion inicial suave
    t, steps = 0.0, 0
    while t < tmax:
        dt = cfl * dx / np.max(np.abs(u))    # CFL>1 permitido (implicito)
        dt = min(dt, tmax - t)
        u = lu_sgs_burgers_step(u, dt, dx)
        t += dt; steps += 1
    return x, u, steps
 
x, u, steps = drive_burgers_shock(cfl=5.0)
print(f"CFL=5 -> {steps} steps,  u in [{u.min():.3f}, {u.max():.3f}]")
# CFL=5 -> 49 steps,  u in [0.012, 0.988]

Corre el mismo problema con Euler explícito hacia adelante a CFL 5 y produce NaN en dos o tres pasos. La versión implícita alcanza el choque en 49 pasos estables. Sin un solver de Riemann pesado por cara: una división diagonal y dos barridos hicieron todo el trabajo.

Barrido tras barrido, el residuo se desploma#

Una sola aplicación de LU-SGS da una solución aproximada. Para resolver con más exactitud, repites el barrido sobre el mismo sistema. La velocidad de convergencia depende entonces de la dominancia diagonal de la matriz. En el caso dominado por la viscosidad —donde el jacobiano viscoso se vuelve un laplaciano 1D— la diagonal es $1+2d$ y las no diagonales son $-d$ . Aquí $d=\nu\Delta t/\Delta x^2$ es el número de difusión (el paso de tiempo adimensional de la difusión viscosa).

Mira abajo cómo cae el residuo barrido a barrido.

diffusion # d = 6.0N = 48

Each sweep is one forward + one backward Gauss-Seidel pass. At d ≈ 1 the residual drops ten orders in a handful of sweeps. Crank d up to 40 (a huge implicit Δt) and the same curve flattens — the matrix is stiffer, so LU-SGS needs more sweeps or a stronger preconditioner to keep up.

Mantén $d$ cerca de 1 y el residuo cae diez órdenes en un puñado de barridos. Sube $d$ hasta 40 (un $\Delta t$ implícito enorme) y la misma curva se aplana. Compraste estabilidad agrandando el paso, pero la matriz se volvió más rígida, así que la relajación interna necesita más barridos. Estabilidad y convergencia interna no llegan gratis a la vez.

Trampas que evitar al escribir el código#

Tres cosas atrapan a la gente en la práctica.

No olvides el radio espectral viscoso en la diagonal. Cuando el término viscoso es grande, debes sumar una contribución $\nu/\Delta x^2$ a la diagonal $D$ . Si la omites, $D$ se encoge y el barrido diverge.

Ten la subrrelajación a mano. Cuando la no linealidad es fuerte, no sumes todo el $\Delta U$ de golpe. Usa $U^{n+1}=U^n+\omega\,\Delta U$ con $\omega\in[0.5,0.8]$ . Apaga las oscilaciones cerca de los choques.

Alterna la dirección del barrido. Repetir solo el barrido hacia adelante deja que la información fluya en un solo sentido. Alternar adelante y atrás mantiene viva la simetría y acelera la convergencia. Cruzando una frontera de partición MPI, inserta una comunicación entre barridos para que la actualización del bloque vecino se propague.

Lo que quiero dejarte: el método implícito abre la región de estabilidad a todo el semiplano izquierdo y concede un $\Delta t$ grande. La matriz gigante que cuesta la aproxima LU-SGS con una separación $D\!\cdot\!L\!\cdot\!U$ y dos barridos, sin almacenamiento ni multiplicaciones. La estabilidad es casi gratis, pero la convergencia interna se paga en proporción a la rigidez de la matriz.