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

La corrida divergió a las 3 de la madrugada. La última línea del registro: pressure residual = NaN. El CFL bajó de 0.9 a 0.5 y luego a 0.3. Cada paso era tan corto que el estado estacionario pedía días. Muro habitual del cálculo estacionario compresible. El avance temporal implícito rompe ese muro por el lado del paso. El cálculo que no invierte la matriz completa es LU-SGS. El choque de Burgers sigue ahí con CFL 5.

El interruptor de esquema pasa de Euler hacia adelante a Euler hacia atrás, y $\mathrm{Re}(\lambda\Delta t)$ se arrastra hasta $-2.5$. La región estable cambia de un disco pequeño al semiplano izquierdo.

Con Euler hacia adelante el punto sale del disco cian y pasa a UNSTABLE. Con Euler hacia atrás todo el semiplano izquierdo queda cian, y la misma coordenada vuelve a ser estable.

La región de estabilidad se abre al infinito#

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

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

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.

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.

El disco explícito termina en 0 y en $-2$ sobre el eje real. Centro $-1$, radio 1. $\mathrm{Re}(\lambda\Delta t)=-2.5$ queda a 1.5 del centro, fuera del disco, y por eso el marcador lee UNSTABLE. Euler hacia atrás usa $G=1/(1-\lambda\Delta t)$. Si $\mathrm{Re}(\lambda)<0$, el denominador supera 1 y $|G|$ se mantiene por debajo de 1 aunque $\Delta t$ crezca. El registro de madrugada recortó el CFL por 0.9, 0.5 y 0.3 para encerrar $\lambda\Delta t$ en ese disco. En una malla viscosa $\lambda$ crece primero. El disco se rompe antes de que la precisión pida un paso más corto.

Invertir la matriz completa es un costo impagable#

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

$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,

queda un sistema lineal para $\Delta U = U^{n+1}-U^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$.

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

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:

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

Con qué se llenan $D$, $L$ y $U$ decide la calidad de la aproximación. 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 se construye la matriz entera: bastan unos pocos escalares (o un bloque pequeño) por celda.

El producto aproximado es $(D+L)\,D^{-1}\,(D+U)$. Al desarrollarlo aparece $L D^{-1} U$ junto a $A$, y ese término se descarta. Como $D_i = I/\Delta t + |a_i|/\Delta x$, un $\Delta t$ pequeño infla $I/\Delta t$ y el término descartado se encoge. Subir el CFL a 5 agranda $\Delta t$ y encoge $I/\Delta t$: crece el error de la aproximación. El avance explícito sigue dando NaN en dos o tres pasos en ese tramo. La salida implícita lo cruza en 49 pasos con $u$ cerrado en $[0.012, 0.988]$, en $n=200$ celdas hasta $t=0.22$ y $\Delta x=1/200$. A diferencia de una comparación de convergencia Krylov que guarda la matriz, aquí una división diagonal y dos barridos ocupan ese lugar.

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]

El mismo problema con Euler explícito hacia adelante a CFL 5 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).

Con $d=\nu\Delta t/\Delta x^2=1$ la diagonal vale $1+2d=3$ y cada no diagonal vale $-d=-1$. La diagonal supera la suma de los dos vecinos, y unos pocos barridos bajan el residuo diez órdenes. Con $d=40$ la diagonal es 81 y las no diagonales son $-40$ dos veces, así que la suma de la fila es $81-80=1$ y casi se cancela. Por eso la curva se acuesta. La estabilidad comprada con un $\Delta t$ grande y el número de barridos que pide la relajación interna van en el mismo libro.

Al mover el control de $d$ desde cerca de 1 hasta 40, cambia la pendiente de la curva de residuo con el mismo número de barridos.

Con $d$ cerca de 1 el residuo cae diez órdenes en un puñado de barridos. Al llevar $d$ hasta 40 (un $\Delta t$ implícito enorme) la misma curva se aplana. La estabilidad se compró 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 juntas sin costo.

Si el radio viscoso sale de la diagonal, el barrido muere#

Sin $\nu/\Delta x^2$ en la diagonal $D$, $D$ se encoge en una celda viscosa. El radio espectral queda solo en el $|a|/\Delta x$ convectivo. El barrido diverge enseguida. Si el jacobiano viscoso es un laplaciano 1D, el dominio diagonal es el balance entre $1+2d$ y $-d$, y quitar el término viscoso vacía el lado izquierdo de ese balance. Cuando la no linealidad es fuerte, no se suma todo el $\Delta U$ de una vez. Se usa $U^{n+1}=U^n+\omega\,\Delta U$ con $\omega\in[0.5,0.8]$. Las oscilaciones cerca de los choques bajan.

La dirección también es una decisión. Repetir solo el barrido hacia adelante manda la información hacia índices crecientes. Alternar adelante y atrás deja simétricos a $L$ y $U$. En una cara de partición MPI tiene que entrar una comunicación entre los dos barridos, o el $\Delta U$ del bloque vecino no llega. Cómo se numeran las celdas cambia ese orden de barrido. El registro en que LU-SGS se vuelve lento tras renumerar con RCM es el otro lado de la misma decisión.