[Paper Review] GMRES Ran 1000 Times and Still Wouldn't Stop — Physics-Based Block Preconditioning for an All-Speed Melt Pool

I ran a laser melt pool implicitly. On the first time step GMRES spun 1000 times and the residual sat at 1e-2. Drop Mach to 1e-3 and the sound speed is 1000 times the flow. The Jacobian got sick on that ratio. Weston, Nourgaliev, Delplanque, and Barker (2019) broke the stall with a velocity–pressure Schur complement as the preconditioner. These notes move that preconditioner into a 2-field toy and record where the iteration count sticks as Mach changes.

Drop the Mach slider and only the acoustic wave on the top track speeds up, at $c=1/M$. The convective tracer below stays at $u=1$.

At $M=10^{-2}$ the acoustic wave laps the track 100 times while the tracer moves one notch. At $M=10^{-3}$ the ratio is 1000:1. Same gap as the night the residual stuck at 1e-2.

Where This Paper Stands#

• Authors: B. Weston, R. Nourgaliev, J.-P. Delplanque, A. T. Barker (LLNL / UC Davis)
• Journal: Journal of Computational Physics, 397, 108847 (2019)
• Title: Preconditioning a Newton-Krylov solver for all-speed melt pool flow physics
• One line: Solve all-speed compressible Navier–Stokes with phase change by fully implicit JFNK, and tame the low-Mach stiffness with a primitive-variable block Schur complement preconditioner.

The usual preconditioners — element-block SOR, monolithic AMG, block Gauss–Seidel — all blow up in iteration count when the time step or Mach number goes to an extreme. What the paper actually ships is the vP-vT Schur complement. The block shape of the physical coupling is copied into the preconditioner.

Why the Jacobian Gets Sick as M Goes to Zero#

"All-speed" means one code solves everything from incompressible ($M\to 0$) to supersonic. The trouble is the low-Mach limit. The ratio of the acoustic speed $c$ to the flow speed $u$ widens as $c/u = 1/M$. At $M=10^{-3}$, sound moves 1000 times faster than the fluid.

Go fully implicit and these two scales mix inside one Jacobian. The pressure–velocity coupling block grows like $1/M^2$ and the matrix's condition number explodes. A Krylov solver crawls on a matrix whose eigenvalues are smeared wide. That's why GMRES wouldn't stop.

Stamp $M$ at the four values in the code and $\beta=1/M^2$ falls out. $M=1$ gives $\beta=1$. $M=10^{-1}$ gives $\beta=10^{2}$. $M=10^{-2}$ gives $\beta=10^{4}$. $M=10^{-3}$ gives $\beta=10^{6}$. One entry in the pressure–velocity block grows by a million. The grid is $n=64$ and the spacing is $dx=1/64$ in every case. The matrix size stays $2n=128$. Only the spectrum scaled by $\beta$ changes.

An acoustic CFL crushes the step as $\Delta t \propto M\Delta x$. Drop $M$ from 1 to $10^{-3}$ and the allowed step shrinks by 1000. Fully implicit lifts that limit and hands the condition number to $\beta$. The same stiffness, amplified the other way by a preconditioner, is written up in block preconditioning at acoustic CFL 10.

JFNK (Jacobian-Free Newton–Krylov) never forms the Jacobian explicitly; it approximates only the matrix–vector product with a Fréchet difference.

$F$ is the nonlinear residual, $\boldsymbol{\kappa}$ a Krylov vector, $h$ a small finite number. Skipping the stored matrix cuts memory. Without a preconditioner, that difference just walks the smeared eigenvalues at low Mach.

Dropping Conservative Variables for the Primitives $(P,v,T)$#

The first move is a change of variables. Instead of the conservative set $U=(\rho, \rho v, E)$, solve for the primitives $W=(P, v, T)$, because the primitives are far better conditioned at low Mach.

$\partial U/\partial W$ is the change-of-variables Jacobian. The residual $F$ is still written to satisfy the conservation laws, so swapping variables does not break mass, momentum, or energy conservation. Now order the degrees of freedom by field into a $3\times3$ block matrix.

Each block corresponds to one primitive field (velocity, pressure, temperature). That ordering is the springboard for the factorization.

Block LU and the Velocity–Pressure Schur Complement#

The pressure–temperature coupling ($M_{PT}, M_{TP}$) is weak. Drop both and the block LU factorization comes out clean. What remains is the velocity–pressure Schur complement.

$S_{vP}$ is the effective equation for pressure after eliminating velocity — the same skeleton as the pressure-correction equation in the SIMPLE family. Take the preconditioner $P$ as this block-upper-triangular form and the eigenvalues of $P^{-1}A$ cluster near 1. With the exact $S$, the $U$ factor of $A=LU$ literally equals $P$, so $P^{-1}A$ is similar to a unipotent lower triangular matrix — every eigenvalue is exactly 1. GMRES finishes in two or three steps.

The code stops at rtol=1e-8, with maxiter=500 and restart 2n=128. Eigenvalues on one point close the Krylov space in two steps. Otherwise the run hits 500 first.

Python — Counting GMRES Iterations as Mach Drops#

Leave out temperature and reproduce the heart in a 2-field $(v,P)$ toy. The block system is $A=\begin{bmatrix} I & \beta G \\ D_v & I\end{bmatrix}$ with $\beta=1/M^2$. The preconditioner is the exact Schur complement $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():
    # periodic forward difference (the pressure-gradient operator)
    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()          # divergence = -G^T  ->  Dv*G ~ Laplacian
 
def all_speed_system(M):
    beta = 1.0 / M**2        # low-Mach stiffness factor
    A = bmat([[eye(n), beta * G], [Dv, eye(n)]]).tocsr()
    return A, beta
 
def schur_preconditioner(beta):
    # velocity-pressure Schur complement 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}")

Without a preconditioner the iteration count climbs into the tens and hundreds as $M$ shrinks, and at $M=10^{-3}$ it never gets inside the cap of 500. The Schur side factors $S=I-\beta D_v G$ once by LU, and GMRES drops the residual to 1e-11 in two or three steps. Equal Schur counts at $M=1$ ($\beta=1$) and $M=10^{-3}$ ($\beta=10^{6}$) are the numerical face of eigenvalues piled at 1.

Interactive — Acoustic and Convective, Two Speeds Pulling Apart#

The 100:1 gap at $M=10^{-2}$ is $c/u=1/M$. The Jacobian receives the square, $\beta=10^{4}$. At $M=10^{-3}$, $\beta=10^{6}$. Explicitly, $\Delta t \propto M\Delta x$.

Interactive — Where the Preconditioner Collapses the Residual#

Now the actual GMRES residual. Move the Mach slider and the two curves — with and without preconditioning — redraw.

The orange curve (no preconditioner) flattens and stalls as $M$ drops. The cyan curve (Schur preconditioned) plunges to 1e-11 in two or three steps no matter where $M$ is. The gap is widest at $M=10^{-3}$. The block scaled by $\beta=10^{6}$ is absorbed whole by a preconditioner copied from the physical coupling.

Three Things I Doubted While Coding It#

An exact $S$ is only free in a toy. Here $S=I-\beta D_v G$ was factored by LU and the solve converged in two steps. In practice, the $M_{vv}^{-1}$ inside $S_{vP}=M_{PP}-M_{Pv}M_{vv}^{-1}M_{vP}$ is a dense inverse you cannot form. The paper solves $S$ approximately (one AMG V-cycle, say), and the quality of that approximation drives the outer iteration count. The cost of a Mach-flat iteration count hides here.

How far to drop the coupling. The paper says dropping $M_{PT}, M_{TP}$ costs no robustness. That assumption is tuned to melt pool physics. With strong compressible thermal coupling — deflagration, combustion — the same truncation can stall convergence. Whether the preconditioner gets simplified is a physics call.

The JFNK Fréchet difference $h$. Too small an $h$ lets round-off in, too large lets truncation error pollute the matrix–vector product. Even with a strong preconditioner, a mistuned $h$ stalls Newton. It is the first snag when you leave OpenFOAM's segregated (SIMPLE) approach for a fully coupled implicit one.

Closing — A Preconditioner Is Not Free#

The night GMRES stopped at 1000 iterations and a residual of 1e-2 did not end by pushing the Jacobian harder. The preconditioner had to copy the coupling itself. The Schur complement is the pressure equation with velocity eliminated, and an exact solve piles the eigenvalues at 1. In a real $S_{vP}$, $M_{vv}^{-1}$ is dense, so the paper approximates it with one AMG V-cycle. A rough cycle opens the outer count back up from two or three steps.

Skipping the stored Jacobian and keeping only the matrix–vector product is the same fork as operator-overloading automatic differentiation. JFNK replaces that product with the Fréchet difference $h$. A bad $h$ stops Newton even when the preconditioner is right. Next measurement on this toy: how far the outer count moves off two or three steps when $S$ is solved loosely by AMG.