Linear solvers
Every solve in the library goes through a small protocol, set_operator(A) then
solve(x, b), so the forward, adjoint, incremental and prior solves can each be
given a different solver.
for attr in ("solver", "solver_fwd_inc", "solver_adj_inc"):
setattr(pde, attr, hm.KrylovSolver(comm, "cg", "amg"))
Letting the size choose
auto_solver() picks a direct solve where one
is affordable and Krylov where it is not:
for attr in ("solver", "solver_fwd_inc", "solver_adj_inc"):
setattr(pde, attr, hm.auto_solver(Vu, comm))
The criterion is the global unknown count, not the rank count. A rule like
comm.size == 1 has it backwards: a one-rank run on a fine mesh is exactly where
a replicated factorization runs out of memory, while a four-rank run on a coarse one
is where an exact solve is both affordable and worth having. The threshold is the
replicated solver’s own refusal size, so the choice never lands on a solver that
would then decline; pass max_direct to change it. Where petsc4py is importable
the distributed factorization is preferred, since it is exact without replicating
the fill.
Krylov with a hypre preconditioner
When the Jacobian is symmetric positive definite, say so instead of setting the
solvers by hand: PDEVariationalProblem(..., spd_jacobian=True) (also on the
time-dependent problem) makes every default solver CG with BoomerAMG. Measured on
the Poisson problem with the same iteration counts, a warm solve is 5-8 % cheaper on
the host at 48^3-64^3 and 27-30 % cheaper with hypre on a device (0.178 against
0.244 s at 64^3, 0.604 against 0.859 s at 96^3 on four L40S), because CG holds no
Krylov basis to orthogonalize against. It is opt-in: the symmetry probe behind
symmetric_jacobian="auto" cannot tell definite from indefinite, and CG on an
indefinite Jacobian fails where GMRES converges.
KrylovSolver wraps MFEM’s CG, GMRES,
FGMRES, BiCGSTAB and MINRES with BoomerAMG, ILU, Jacobi, Gauss-Seidel, ParaSails or
Euclid. For vector-valued problems pass systems_dim (and elasticity=True
with the space) so BoomerAMG uses its systems-of-PDEs mode. The elasticity options exist
with hypre on the host only: on a GPU MFEM has none, and the solver keeps the systems
options with a warning. BoomerAMG also has to know
which row belongs to which component, and it is told: component after component
(systems_bynodes, the default, which is the ordering of this library’s spaces unless
another is asked for; read from the space where one is given) or interleaved
(systems_bynodes=False, for mfem.Ordering.byVDIM). Told the wrong one, it pairs
rows of different components and a solve takes two to five times the iterations.
parameters["amg_relax_type"]
selects hypre’s relaxation (-1, the default, keeps MFEM’s choice, l1-Jacobi on a
device); 16, Chebyshev, measured 12 % faster per CG + AMG solve of a 2.1 M-dof
SPD operator on L40S cards, on one card and on four alike, with half the CG
iterations. It is not the default because it needs an SPD operator (the geothermal
model’s Jacobian is not) and the gain is a constant, not a change in scaling.
parameters["amg_options"] takes further BoomerAMG settings by hypre’s names
(HIPPYMFEM_AMG for every new solver). HIPPYMFEM_AMG_CG (hm.config.amg_cg)
adds settings for new CG solvers alone, after those: the operator of a CG solve is
symmetric positive definite, where l1-Jacobi with any weight below two converges, so
HIPPYMFEM_AMG_CG=relax_wt=1.8 is safe there. It saved 21 to 31 % of the iterations
of twelve SPD operators on one L40S with hypre on the card: the model problem, a
coefficient contrast of 1e6, anisotropy of 1e-3 aligned with the mesh and rotated, a
mesh stretched twenty times, quadratic tetrahedra, linear elasticity with systems AMG up
to a Poisson ratio of 0.45, and the prior’s operator. The savings were the same on two
cards with twice the elements. The same weight in HIPPYMFEM_AMG, which reaches
GMRES too, stalled GMRES on advection-diffusion at a mesh Peclet number of 5 to 6
(1000 iterations with the residual at 0.3 to 0.6, where MFEM’s settings converged in 150
to 192).
CG with a hypre preconditioner runs in hypre’s own PCG (MFEM’s HyprePCG), not in
MFEM’s CGSolver (HIPPYMFEM_HYPRE_PCG=0 or hm.config.hypre_pcg = False for the
latter). The two make the same iteration: the same iterates, counts and stopping test
from a zero initial guess. (Where hypre’s PCG leaves an iteration that is none, with a
preconditioner or a matrix that is not positive definite or a right-hand side that is
not a number, the solver takes hypre’s own verdict and a residual that is a number,
and reports the failure that MFEM’s CG reports.) The difference is in the V-cycle.
hypre’s PCG marks the
vector it hands to BoomerAMG as zero, and the first relaxation on the finest level then
skips its matrix-vector product; MFEM zeroes the vector without saying so. With
l1-Jacobi relaxation on a GPU that is one product of four per iteration: a solve of 2.1
million dofs to 1e-12 in 24 iterations took 0.103 s with MFEM’s CG and 0.089 s with
hypre’s on an H100, and 0.274 and 0.227 s on an L40S. The solvers of one operator (the
forward and the two incremental ones) share one PCG object, which sets the
preconditioner up once. A solve from a nonzero initial guess, a preconditioner kept
from an older matrix (pc_reuse) and a transposed operator keep MFEM’s CG.
With a single-precision build of hypre next to the double-precision one
(HIPPYMFEM_HYPRE_SINGLE), the Jacobian of a PDE problem, its BoomerAMG hierarchy and
the CG solves with it live in that library; Single precision in the GPU guide
has the measurements, and A single-precision hypre the build.
One behaviour differs deliberately from MFEM’s default. MFEM’s hypre wrappers use
ABORT_HYPRE_ERRORS, so a failed AMG setup ends the whole MPI job. Inside an
optimizer that is the wrong response: a line search routinely proposes a parameter
for which the operator is hopeless (exp(m) spanning hundreds of decades is
finite but defeats AMG’s coarsening), and the right answer is to reject that step.
KrylovSolver sets WARN_HYPRE_ERRORS and treats a “converged” result that
contains non-finite values as non-convergence, so the step becomes a catchable
RuntimeError the line search handles.
Exact solves, including in parallel
LUSolver is the stand-in for hIPPYlib’s
PETScLUSolver. On one rank it factorizes in place with scipy’s SuperLU; in
parallel it gathers the matrix onto rank 0, factorizes it there and scatters each
solution (replicate=True factorizes on every rank instead, so that no rank waits
on rank 0). That is exact and rank-count independent, verified to 1e-15 by
comparing a functional of the solution across 1, 2 and 4 rank runs, at the price of
a serial factorization: one factorization of the whole matrix in time,
O(nnz(L+U)) on rank 0 in memory. It refuses matrices above
max_global_size (400 000 unknowns by default) rather than silently exhausting
memory.
Why bother: an exact solve makes the reduced Hessian exact, so a spectrum or a tolerance study measures the discretization instead of a Krylov tolerance.
Genuinely distributed direct solvers
PETScLUSolver uses the best
factorization the local PETSc was built with (MUMPS, SuperLU_dist, STRUMPACK or
PaStiX), and otherwise PETSc’s redundant with an inner lu, which is exact as
well but serial; package_used says which one a solver got. The petsc4py on
PyPI configures its PETSc without any of the four, and a conda-forge PETSc brings
an MPI of its own that cannot share a process with the one PyMFEM links, so
tools/install_petsc_mumps.sh builds PETSc from source against the system MPI
with MUMPS, ScaLAPACK, METIS and ParMETIS, and petsc4py against it (a venv made
with --system-site-packages on top of the PyMFEM interpreter keeps the PyPI
petsc4py untouched). What MUMPS costs on a 3D problem, measured on the P2
Laplacian-plus-mass of an \(n^3\) hexahedral mesh with that build (system OpenMPI,
one node): at \(32^3\) (274 625 unknowns) on four ranks the factorization takes
105 s and a solve 0.14 s, against 0.84 s for CG with BoomerAMG to 1e-12; at
\(40^3\) (531 441 unknowns) on eight ranks 189 s and 0.18 s against 0.84 s. The
factorization pays for itself after a few hundred solves, which a Laplace
approximation with a hundred Hessian actions is close to, and the solve is exact
(residual 1e-14) where the Krylov solve has a tolerance.
PETScKrylovSolver exposes PETSc’s Krylov
methods and preconditioners on the same matrices.
Warning
petsc4py must be imported before PyMFEM. In the reference environment the
reverse order fails with libmpi_mpifh.so.40: undefined symbol:
mpi_conversion_fn_null_, because the system MPI Fortran library is found first.
Set HIPPYMFEM_PETSC=1 before importing hippymfem, or import petsc4py at the
top of the script. hm.petsc_available() returns (True, "") or
(False, reason) so a script can report the reason rather than fail
mysteriously.
If MFEM itself was built with MUMPS, SuperLU_dist or STRUMPACK, those are reachable
through assemble_native_matrix and MFEM’s own solver classes;
hm.mfem_config() says which are present.
A solve that only preconditions
A Krylov solve stopped at a tolerance is a different map for every right-hand side.
Where a solve is the preconditioner of another iteration, that iteration either needs
it far more accurate than its own tolerance or needs one fixed linear map of any
accuracy. ChebyshevSolver is the second: it
takes a KrylovSolver whose operator is set
and applies its preconditioner (a multigrid cycle) degree times, each from the
iterate before it, with the coefficients of the Chebyshev iteration on the interval of
the cycle’s eigenvalues, which it estimates once. The result is symmetric and positive
definite where the cycle is, and within error_bound of the solve in the energy norm
(1.6 % with four cycles and 0.14 % with six for a cycle with which CG takes 22
iterations to 1e-13). The prior’s precision solver made of two of them is the
preconditioner that ReducedSpaceNewtonCG takes with cg_preconditioner_cycles
(Finding the MAP point).
ks = hm.KrylovSolver(comm, "cg", "amg")
ks.set_operator(A)
fast = hm.ChebyshevSolver(ks, degree=5)
fast.solve(x, b) # five cycles
print(fast.interval, fast.error_bound)
Choosing
situation |
solver |
|---|---|
large problem, well-conditioned operator |
|
the same, as a preconditioner only |
|
exact Hessian or spectrum needed |
|
exact and larger |
|
mass-matrix-like operator, cheap |
|