[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

Los autovalores valen 1 solo con el complemento de Schur exacto. Una resolución floja de S por AMG vuelve a abrir el conteo externo.

Se corrió un melt pool láser en forma implícita. En el primer paso de tiempo, GMRES dio 1000 vueltas y el residuo se quedó en 1e-2. Al bajar el Mach a 1e-3, la velocidad del sonido pasa a ser 1000 veces la del flujo. El jacobiano enfermó en esa razón. Weston, Nourgaliev, Delplanque y Barker (2019) cortaron el estancamiento usando el complemento de Schur velocidad–presión como precondicionador. Estas notas pasan ese precondicionador a un toy de 2 campos y anotan dónde se pega el número de iteraciones al cambiar el Mach.

Al bajar el control de Mach, solo la onda acústica de la pista superior se acelera, a $c=1/M$ . El trazador convectivo de abajo se queda en $u=1$ .

Mach number M = 3.2e-2

explicit acoustic CFL forces Δt ∝ M·Δx — the smaller M, the stiffer the coupled system.

En $M=10^{-2}$ la onda acústica da 100 vueltas a la pista mientras el trazador avanza una muesca. En $M=10^{-3}$ la razón es 1000:1. La misma brecha de la noche en que el residuo se pegó en 1e-2.

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. Lo que entrega el artículo es el complemento de Schur vP-vT. La forma en bloques del acoplamiento físico se copia en la matriz del precondicionador.

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 autovalores muy esparcidos.

Al colocar $M$ en los cuatro puntos del código, $\beta=1/M^2$ sale directo. $M=1$ da $\beta=1$ . $M=10^{-1}$ da $\beta=10^{2}$ . $M=10^{-2}$ da $\beta=10^{4}$ . $M=10^{-3}$ da $\beta=10^{6}$ . Una entrada del bloque presión–velocidad crece un millón de veces. La malla es $n=64$ y el espaciado $dx=1/64$ en los cuatro casos. El tamaño de la matriz sigue en $2n=128$ . Solo cambia el espectro escalado por $\beta$ .

Un CFL acústico aplasta el paso según $\Delta t \propto M\Delta x$ . Bajar $M$ de 1 a $10^{-3}$ reduce el paso admisible a la milésima parte. Lo totalmente implícito suelta ese límite y entrega el número de condición a $\beta$ . La misma rigidez, amplificada al revés por un precondicionador, está en el precondicionador por bloques con CFL acústico 10.

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.

J\boldsymbol{\kappa} \approx \frac{F(x + h\boldsymbol{\kappa}) - F(x)}{h}

$F$ es el residuo no lineal, $\boldsymbol{\kappa}$ un vector de Krylov y $h$ un número finito pequeño. No guardar la matriz reduce memoria. Sin precondicionador, esa diferencia camina los autovalores esparcidos en bajo Mach.

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.

\frac{\partial F(x)}{\partial U}\frac{\partial U}{\partial W}\,\delta W = -F(x)

$\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$ .

M = \begin{bmatrix} M_{vv} & M_{vP} & M_{vT} \\ M_{Pv} & M_{PP} & M_{PT} \\ M_{Tv} & M_{TP} & M_{TT} \end{bmatrix}

Cada bloque corresponde a un campo primitivo (velocidad, presión, temperatura). Ese orden es el trampolín de la factorización.

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 lo que queda es el complemento de Schur velocidad–presión.

S_{vP} = M_{PP} - M_{Pv}\,M_{vv}^{-1}\,M_{vP}

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

El código para en rtol=1e-8, con tope maxiter=500 y reinicio 2n=128. Autovalores en un punto cierran el subespacio de Krylov en dos pasos. Si no, la corrida toca 500 antes.

Python — contar iteraciones de GMRES al bajar el Mach#

La temperatura queda fuera. El corazón se reproduce 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}")

Sin precondicionador, el número de iteraciones trepa a decenas y cientos a medida que $M$ encoge, y en $M=10^{-3}$ ni siquiera entra en el tope de 500. El lado de Schur factoriza $S=I-\beta D_v G$ una vez por LU, y GMRES baja el residuo a 1e-11 en dos o tres pasos. Conteos de Schur iguales en $M=1$ ( $\beta=1$ ) y en $M=10^{-3}$ ( $\beta=10^{6}$ ) son la cara numérica de autovalores apilados en 1.

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

La brecha 100:1 en $M=10^{-2}$ es $c/u=1/M$ . El jacobiano recibe el cuadrado, $\beta=10^{4}$ . En $M=10^{-3}$ , $\beta=10^{6}$ . En forma explícita, $\Delta t \propto M\Delta x$ .

Interactivo — dónde el precondicionador derrumba el residuo#

Ahora el residuo real de GMRES. Al mover el control de Mach, las dos curvas — con y sin precondicionamiento — se redibujan.

Mach number M = 1.0e-2 (β = 1/M² = 1.0e+4)

unpreconditioned: 0 iters to 1.0e+0 · Schur: 0 iters to 1.0e+0

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 es máxima en $M=10^{-3}$ . El bloque escalado por $\beta=10^{6}$ queda absorbido entero por un precondicionador copiado del acoplamiento físico.

Tres cosas que generaron dudas al programarlo#

Un $S$ exacto solo es gratis en un toy. Aquí $S=I-\beta D_v G$ se factorizó por LU y la convergencia llegó 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 de un conteo plano en Mach se esconde aquí.

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 frenar la convergencia. Si el precondicionador se simplifica, lo decide la física.

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 el primer tropiezo al pasar del enfoque segregado (SIMPLE) de OpenFOAM a uno totalmente acoplado e implícito.

Cierre — un precondicionador no es gratis#

La noche en que GMRES se detuvo en 1000 iteraciones y un residuo de 1e-2 no se cerró apretando más el jacobiano. Hacía falta un precondicionador que copiara el propio acoplamiento. El complemento de Schur es la ecuación de presión con la velocidad eliminada, y resolverlo de forma exacta agrupa los autovalores en 1. En un $S_{vP}$ real, $M_{vv}^{-1}$ es denso, así que el artículo lo aproxima con un V-cycle de AMG. Un ciclo áspero vuelve a abrir el conteo externo desde dos o tres pasos.

No guardar el jacobiano y quedarse solo con el producto matriz–vector es la misma bifurcación que la diferenciación automática por sobrecarga de operadores. JFNK sustituye ese producto por la diferencia de Fréchet $h$ . Un $h$ desviado detiene a Newton aunque el precondicionador esté bien. La medición siguiente, en el mismo toy, es cuánto se aleja el conteo externo de dos o tres pasos cuando $S$ se resuelve de forma floja con AMG.