[Reseña de artículo] GMRES corrió 1000 veces y no se detenía — precondicionamiento por bloques físico para un melt pool a toda velocidad

Corrí una simulación implícita de un melt pool (charco fundido) por láser. En el primer paso de tiempo, GMRES dio 1000 vueltas y el residuo se quedó clavado en 1e-2. Al bajar el número de Mach a 1e-3, el síntoma empeoró — la velocidad del sonido era 1000 veces la del flujo, y el jacobiano estaba enfermo. Weston, Nourgaliev, Delplanque y Barker (2019) resolvieron justo esto. El truco es usar el complemento de Schur velocidad–presión como precondicionador. Hoy reproduzco su núcleo en un toy de 2 campos y observo cómo un solo precondicionador devuelve la convergencia independiente del Mach.

Dónde se sitúa este artículo#

• Autores: B. Weston, R. Nourgaliev, J.-P. Delplanque, A. T. Barker (LLNL / UC Davis)
• Revista: Journal of Computational Physics, 397, 108847 (2019)
• Título: Preconditioning a Newton-Krylov solver for all-speed melt pool flow physics
• En una línea: resolver las Navier–Stokes compresibles a toda velocidad (all-speed) con cambio de fase mediante JFNK totalmente implícito, y domar la rigidez de bajo Mach con un precondicionador de complemento de Schur por bloques en variables primitivas.

Los precondicionadores habituales — SOR por bloque de elemento, AMG monolítico, Gauss–Seidel por bloques — disparan todos el número de iteraciones cuando el paso de tiempo o el Mach llegan a un extremo. El aporte del artículo es el complemento de Schur vP-vT, un precondicionador que copia directamente la estructura del acoplamiento físico.

Por qué el jacobiano enferma cuando M tiende a cero#

"A toda velocidad" significa que un mismo código resuelve desde lo incompresible ($M\to 0$) hasta lo supersónico. El problema es el límite de bajo Mach. La razón entre la velocidad acústica $c$ y la del flujo $u$ se abre como $c/u = 1/M$. En $M=10^{-3}$, el sonido va 1000 veces más rápido que el fluido.

Al ir totalmente implícito, esas dos escalas se mezclan dentro de un único jacobiano. El bloque de acoplamiento presión–velocidad crece como $1/M^2$ y el número de condición de la matriz se dispara. Un solver de Krylov se arrastra sobre una matriz cuyos autovalores están muy esparcidos. Por eso GMRES no se detenía.

JFNK (Jacobian-Free Newton–Krylov) nunca forma el jacobiano de forma explícita; solo aproxima el producto matriz–vector con una diferencia de Fréchet.

$F$ es el residuo no lineal, $\boldsymbol{\kappa}$ un vector de Krylov y $h$ un número finito pequeño. No almacenar la matriz ahorra memoria, pero sin precondicionador la convergencia muere a bajo Mach. Precondicionar no es opcional aquí: es todo el juego.

Soltar las variables conservativas por las primitivas $(P,v,T)$#

El primer movimiento es un cambio de variables. En vez del conjunto conservativo $U=(\rho, \rho v, E)$, se resuelve por las primitivas $W=(P, v, T)$, porque están mucho mejor condicionadas a bajo Mach.

$\partial U/\partial W$ es el jacobiano del cambio de variables. El residuo $F$ se sigue escribiendo para satisfacer las leyes de conservación, así que cambiar de variables no rompe la conservación de masa, momento ni energía. Ahora se ordenan los grados de libertad por campo en una matriz por bloques $3\times3$.

Cada bloque corresponde a un campo primitivo (velocidad, presión, temperatura). Esta estructura es el trampolín del siguiente paso.

LU por bloques y el complemento de Schur velocidad–presión#

El acoplamiento presión–temperatura ($M_{PT}, M_{TP}$) es débil. Al descartar ambos, la factorización LU por bloques sale limpia y el núcleo se condensa en el complemento de Schur velocidad–presión.

$S_{vP}$ es la ecuación efectiva para la presión tras eliminar la velocidad — el mismo esqueleto que la ecuación de corrección de presión de la familia SIMPLE. Toma el precondicionador $P$ con esta forma triangular superior por bloques y los autovalores de $P^{-1}A$ se agrupan cerca de 1. Con el $S$ exacto, el factor $U$ de $A=LU$ es literalmente $P$, de modo que $P^{-1}A$ es semejante a una triangular inferior unipotente: todos los autovalores valen exactamente 1. GMRES termina en dos o tres pasos.

Python — contar iteraciones de GMRES al bajar el Mach#

Deja fuera la temperatura y reproduce el corazón en un toy de 2 campos $(v,P)$. El sistema por bloques es $A=\begin{bmatrix} I & \beta G \\ D_v & I\end{bmatrix}$ con $\beta=1/M^2$. El precondicionador es el complemento de Schur exacto $S=I-\beta D_v G$.

import numpy as np
from scipy.sparse import diags, eye, bmat
from scipy.sparse.linalg import gmres, LinearOperator, splu
 
n = 64
dx = 1.0 / n
 
def fwd_diff():
    # diferencia hacia adelante periodica (operador de gradiente de presion)
    D = diags([-1.0, 1.0], [0, 1], shape=(n, n)).tolil()
    D[-1, 0] = 1.0
    return (D / dx).tocsr()
 
G = fwd_diff()
Dv = (-G.T).tocsr()          # divergencia = -G^T  ->  Dv*G ~ laplaciano
 
def all_speed_system(M):
    beta = 1.0 / M**2        # factor de rigidez de bajo Mach
    A = bmat([[eye(n), beta * G], [Dv, eye(n)]]).tocsr()
    return A, beta
 
def schur_preconditioner(beta):
    # complemento de Schur velocidad-presion S = M_PP - M_Pv M_vv^{-1} M_vP = I - beta Dv G
    S = (eye(n) - beta * (Dv @ G)).tocsc()
    luS = splu(S)
    def apply(r):
        ru, rp = r[:n], r[n:]
        zp = luS.solve(rp)               # S z_p = r_p
        zu = ru - beta * (G @ zp)        # z_u = r_u - beta G z_p
        return np.concatenate([zu, zp])
    return LinearOperator((2 * n, 2 * n), matvec=apply)
 
def count_iters(M, use_prec):
    A, beta = all_speed_system(M)
    b = np.random.default_rng(0).standard_normal(2 * n)
    P = schur_preconditioner(beta) if use_prec else None
    it = {"k": 0}
    x, info = gmres(A, b, M=P, rtol=1e-8, maxiter=500, restart=2 * n,
                    callback=lambda xk: it.__setitem__("k", it["k"] + 1))
    return it["k"]
 
print(f"{'Mach':>8}  {'no-prec':>8}  {'Schur':>6}")
for M in [1.0, 1e-1, 1e-2, 1e-3]:
    print(f"{M:>8.0e}  {count_iters(M, False):>8d}  {count_iters(M, True):>6d}")

El resultado es contundente. Sin precondicionador, el número de iteraciones trepa a decenas y cientos a medida que $M$ encoge, y en $M=10^{-3}$ ni siquiera converge dentro del tope. El precondicionador de Schur termina en dos o tres pasos sin importar $M$. La "convergencia independiente del Mach" que reporta el artículo nace de este precondicionador de una sola línea.

Interactivo — acústica y convección, dos velocidades que se separan#

Primero, mira de dónde sale la rigidez. Baja el número de Mach en la simulación de abajo. La onda acústica de la pista superior corre a $c=1/M$, y el trazador convectivo de abajo a $u=1$.

Bájalo a $M=10^{-2}$ y la onda acústica da 100 vueltas a la pista mientras el trazador convectivo avanza una sola muesca. Esa brecha de velocidad 100:1 se traslada directamente al número de condición del jacobiano. De forma explícita, el CFL acústico aplastaría el paso de tiempo hasta $\Delta t \propto M\Delta x$.

Interactivo — dónde el precondicionador derrumba el residuo#

Ahora observa el residuo real de GMRES. Mueve el deslizador de Mach abajo y las dos curvas — con y sin precondicionamiento — se redibujan.

La curva naranja (sin precondicionador) se aplana y se estanca al bajar $M$. La curva cian (precondicionada con Schur) se desploma a 1e-11 en dos o tres pasos, esté donde esté $M$. La brecha entre ambas es más dramática en $M=10^{-3}$ — es el precondicionador, construido a partir del acoplamiento físico, absorbiendo la rigidez por completo.

Tres cosas que me hicieron dudar al programarlo#

Primero, un $S$ exacto solo es gratis en un toy. Aquí resolví $S=I-\beta D_v G$ de forma exacta por LU y convergí en dos pasos. En la práctica, el $M_{vv}^{-1}$ dentro de $S_{vP}=M_{PP}-M_{Pv}M_{vv}^{-1}M_{vP}$ es una inversa densa que no se puede formar. El artículo resuelve $S$ de forma aproximada (un V-cycle de AMG, por ejemplo), y la calidad de esa aproximación gobierna el número de iteraciones externas. El costo real de la "independencia del Mach" se esconde justo aquí.

Segundo, hasta dónde descartar el acoplamiento. El artículo dice que descartar $M_{PT}, M_{TP}$ no cuesta robustez. Es un supuesto ajustado a la física del melt pool. Con un acoplamiento térmico compresible fuerte — deflagración, combustión — el mismo truncamiento puede matar la convergencia. Simplificar el precondicionador depende de la física.

Tercero, la diferencia de Fréchet $h$ de JFNK. Un $h$ demasiado pequeño deja entrar el redondeo; demasiado grande deja que el error de truncamiento contamine el producto matriz–vector. Aun con un buen precondicionador, un $h$ mal ajustado estanca a Newton. Es la primera trampa al pasar del enfoque segregado (SIMPLE) de OpenFOAM a uno totalmente acoplado e implícito.

Cierre — un precondicionador no es gratis#

La lección de esa noche en que GMRES no se detenía fue simple. La rigidez de bajo Mach no se resuelve apretando más el jacobiano. Hace falta un precondicionador que copie la propia estructura del acoplamiento. El complemento de Schur es la ecuación efectiva para la presión, y resolverlo de forma exacta agrupa los autovalores en 1. Pero el $S$ del mundo real es una aproximación, y cuán barata y exactamente lo resuelves es la próxima batalla. La próxima vez agregaré una capa más a este mismo toy y observaré qué pasa con el número de iteraciones externas cuando $S$ se aproxima con AMG.