CMPS Theory & User Reference Manual › Sparse Linear Solver

9 Sparse Linear Solver

The global coupled correction has the block-sparse form \(\mathbf A\,\delta\mathbf q=\mathbf B\). Each matrix block represents the sensitivity of one cell residual to its own or a neighboring cell's primitive-variable correction. The GUI linear-solver options change how this system is solved; they do not redefine the finite-volume governing equations.

9.1 Block-sparse structure

For cell \(P\), \[\mathbf{A}_{PP}\delta\mathbf{q}_P +\sum_{N\in\mathcal{N}(P)}\mathbf{A}_{PN}\delta\mathbf{q}_N =\mathbf{B}_P. \](9.1) The dense block structure is an advantage: coupling among momentum, continuity, energy, turbulence, particles, and species is retained locally, while sparsity is governed by mesh connectivity rather than by the square of the total cell count.

The active CMPS linear-algebra formulation provides a CPU path based on Intel MKL and a GPU path based on NVIDIA CUDA. Both operate on the same assembled coupled block system. Linear-solver options are runtime-configurable; this guide does not assume that one Krylov/preconditioner combination is universally optimal for every equation set.

For a current linear iterate \(\delta\mathbf q^{(k)}\), the algebraic residual is

\[\mathbf r^{(k)}=\mathbf B-\mathbf A\,\delta\mathbf q^{(k)}.\]

A normalized linear convergence measure can therefore be written as

\[\eta_k=\frac{\lVert\mathbf r^{(k)}\rVert_2}{\max(\lVert\mathbf B\rVert_2,\epsilon)}.\]

9.2 CMPS scalable solver layer

The coupled matrix is solved through the CMPS scalable linear-solver interface. On CPUs, sparse/block algebra and solver kernels are provided through the Intel MKL formulation; on NVIDIA GPUs, the corresponding accelerator path uses CUDA. The interface exposes several solver families without changing the finite-volume assembly: conventional Krylov methods with local preconditioners, algebraic multigrid (AMG) as a solver, and AMG-preconditioned Krylov methods.

This separation is important. The nonlinear CFD layer constructs the same block system \[\mathbf{A}\,\delta\mathbf{q}=\mathbf{B},\] while the scalable-solver layer determines how that system is inverted. Changing from ILU-preconditioned GMRES to AMG-preconditioned FGMRES, for example, does not change the residual or its AD Jacobian.

9.3 Algebraic multigrid (AMG)

AMG constructs a hierarchy directly from the algebraic operator rather than from the CFD mesh hierarchy. Denote the fine-grid matrix by \(\mathbf{A}_0\). AMG generates successively smaller operators \[\mathbf{A}_{\ell+1}=\mathbf{R}_{\ell}\mathbf{A}_{\ell}\mathbf{P}_{\ell}, \](9.2) where \(\mathbf{P}_{\ell}\) is the algebraic prolongation and \(\mathbf{R}_{\ell}\) is restriction. A multilevel application combines smoothing of high-frequency algebraic error on each level with coarse-grid correction of error components that are slow for the smoother.

CMPS exposes an AMG-only path, denoted by the formulation option AGGAMG. In this mode the multilevel hierarchy itself provides the linear correction, without requiring an outer Krylov iteration. The formulation provides configurable AMG depth; the present CMPS option set includes a default maximum of five AMG levels unless changed by the case settings. The same CMPS AMG abstraction can be used by the MKL CPU solver path or by the CUDA accelerator path according to the selected build and execution target.

AMG is distinct from the agglomeration/geometric multigrid described in Chapter 18. The latter coarsens the physical control-volume problem and can rediscretize the governing equations. AMG coarsens the assembled linear algebra of one implicit solve. CMPS can therefore use geometric/agglomeration multigrid at the nonlinear CFD level while still using AMG or an AMG-preconditioned Krylov method inside a level’s linear solve.

One multigrid correction step starts from a smoothed approximation \(\mathbf x_\ell\), forms the level residual and restricts it,

\[\mathbf r_\ell=\mathbf b_\ell-\mathbf A_\ell\mathbf x_\ell,\qquad \mathbf r_{\ell+1}=\mathbf R_\ell\mathbf r_\ell,\]

then solves or approximately solves the coarse error equation

\[\mathbf A_{\ell+1}\mathbf e_{\ell+1}=\mathbf r_{\ell+1},\]

and prolongs the correction back to the finer level,

\[\mathbf x_\ell\leftarrow\mathbf x_\ell+\mathbf P_\ell\mathbf e_{\ell+1}.\]

For an algebraically generated hierarchy, the coarse operator has the standard triple-product form

\[\mathbf A_{\ell+1}=\mathbf R_\ell\mathbf A_\ell\mathbf P_\ell.\]

9.4 AMG-preconditioned Krylov solvers

The CMPS AMGPCKSS family uses AMG as the preconditioner of a Krylov-space method rather than as a stand-alone solve. With right preconditioning, \[\mathbf{A}\mathbf{M}_{\mathrm{AMG}}^{-1}\mathbf{y}=\mathbf{B}, \qquad \delta\mathbf{q}=\mathbf{M}_{\mathrm{AMG}}^{-1}\mathbf{y}. \](9.3) The purpose of \(\mathbf{M}_{\mathrm{AMG}}^{-1}\) is to reduce the long-range/low-frequency error components that make a purely local ILU or block-Jacobi preconditioner deteriorate as the mesh and MPI partition count grow.

The scalable-solver options include flexible Krylov choices, notably FGMRES and a flexible BiCGSTAB-type path, so that the effective AMG preconditioner is allowed to vary between iterations. This is appropriate when the multilevel application contains level-dependent iterative smoothing or when the preconditioner is not mathematically stationary.

The AMG configuration exposes, among other controls:

AMG and preconditioner-AMG level counts (currently configurable up to the selected case limit, with five levels used by the existing default configuration),
V-, W-, and F-cycle selection,
aggregation/coarsening thresholds,
numbers of pre- and post-smoothing sweeps,
ILU, block-Jacobi, symmetric Gauss–Seidel, or additive-Schwarz-type relaxation choices where supported by the selected execution mode,
Krylov smoothers including GMRES, FGMRES, BiCGSTAB, and TFQMR families, and
ILU fill level for local subdomain factorizations.

The distinction between AGGAMG and AMGPCKSS should therefore be kept explicit: the former uses aggregation AMG as the primary linear solve, whereas the latter embeds AMG inside an outer Krylov iteration.

With a possibly iteration-dependent AMG preconditioner, a flexible Krylov method applies

\[\mathbf z_j=\mathbf M_j^{-1}\mathbf v_j,\qquad \mathbf w_j=\mathbf A\mathbf z_j,\]

and stores the preconditioned basis vectors so that the correction has the form

\[\delta\mathbf q_m=\delta\mathbf q_0+\mathbf Z_m\mathbf y_m,\qquad \mathbf Z_m=[\mathbf z_1,\ldots,\mathbf z_m].\]

This differs from fixed-preconditioner GMRES because \(\mathbf M_j\) is allowed to change between Krylov iterations.

9.5 Local block preconditioning

For MPI calculations, a simpler scalable option is a block-Jacobi decomposition in which each MPI rank approximately solves its local diagonal subproblem. In CMPS this is represented by the PCKSS family, with local ILU available inside each block. If \[\mathbf{A}=\mathbf{D}+\mathbf{E},\] where \(\mathbf{D}=\mathrm{blockdiag}(\mathbf{A}_{11},\ldots,\mathbf{A}_{N_rN_r})\) contains rank-local diagonal submatrices, then a block-Jacobi preconditioner is \[\mathbf{M}_{\mathrm{BJ}}^{-1} =\mathrm{blockdiag}(\widetilde{\mathbf{A}}_{11}^{-1},\ldots, \widetilde{\mathbf{A}}_{N_rN_r}^{-1}),\] where each local inverse can be approximated by ILU. This has low communication cost but does not by itself provide the global coarse correction available from AMG.

For a domain decomposition with rank-local diagonal blocks \(\mathbf A_{rr}\), the block-Jacobi preconditioner is

\[\mathbf M_{BJ}=\operatorname{diag}(\mathbf A_{11},\mathbf A_{22},\ldots,\mathbf A_{N_rN_r}),\]

with the preconditioner application defined by independent local solves

\[\mathbf A_{rr}\mathbf z_r=\mathbf r_r,\qquad r=1,\ldots,N_r.\]

9.6 GMRES

For a nonsymmetric Jacobian, GMRES constructs an approximate correction in the Krylov space \[\mathcal{K}_m(\mathbf{M}^{-1}\mathbf{A},\mathbf{M}^{-1}\mathbf{r}_0) =\mathrm{span}\left\{ \mathbf{z}_0, (\mathbf{M}^{-1}\mathbf{A})\mathbf{z}_0, \ldots, (\mathbf{M}^{-1}\mathbf{A})^{m-1}\mathbf{z}_0 \right\},\] where \(\mathbf{M}\) is the preconditioner. The method minimizes the residual norm over the current Krylov space and is therefore suitable for the nonsymmetric block matrices produced by upwind finite-volume fluxes.

Restarted GMRES limits the basis length to control memory. A smaller restart reduces storage but can slow convergence when important low-frequency/error components are discarded at restart.

The Arnoldi process constructs an orthonormal basis \(\mathbf V_m\) and upper-Hessenberg matrix \(\overline{\mathbf H}_m\) such that

\[\mathbf A\mathbf V_m=\mathbf V_{m+1}\overline{\mathbf H}_m.\]

GMRES chooses the Krylov coefficients from the least-squares problem

\[\mathbf y_m=\arg\min_{\mathbf y}\left\|\beta\mathbf e_1-\overline{\mathbf H}_m\mathbf y\right\|_2,\qquad \beta=\lVert\mathbf r_0\rVert_2,\]

and updates

\[\delta\mathbf q_m=\delta\mathbf q_0+\mathbf V_m\mathbf y_m.\]

9.7 Incomplete-LU preconditioning

An ILU factorization approximates \[\mathbf{A}\approx\mathbf{L}\mathbf{U}\] without allowing all fill generated by exact Gaussian elimination. ILU(0) retains the original sparsity pattern; ILU(\(k\)) permits progressively larger levels of fill. More fill generally improves the preconditioner but increases setup cost, memory, and communication/parallel complexity.

For a cell-coupled CFD matrix it is important to preserve the dense equation block when applying a local factorization. Treating each scalar equation independently can destroy the very pressure-velocity/energy coupling introduced by the fully coupled discretization.

The preconditioner application consists of triangular substitutions

\[\mathbf L\mathbf y=\mathbf r,\qquad \mathbf U\mathbf z=\mathbf y,\qquad \mathbf z\approx\mathbf A^{-1}\mathbf r.\]

For block ILU, every retained scalar entry in these expressions represents a dense coupled variable block; this preserves local pressure–velocity–energy/species coupling.

9.8 Linear versus nonlinear accuracy

The linear residual is not the physical CFD residual. It measures how accurately the current linearized correction equation is solved. If the nonlinear state is still far from the solution or the pseudo-time term is strongly stabilizing the system, an excessively tight Krylov tolerance may spend work solving an approximation that will immediately be rebuilt. Conversely, a very loose solve can prevent the outer iteration from following a useful Newton direction.

A practical strategy therefore coordinates:

pseudo-time CFL/diagonal stabilization,
nonlinear residual reduction,
Krylov relative/absolute tolerance,
maximum Krylov iterations, and
preconditioner strength.

No asymptotic complexity claim such as “quadratic becomes linear” should be inferred solely from the storage format; actual cost depends on mesh connectivity, block size, fill, Krylov iteration count, and parallel communication.

9.9 Multigrid relationship

Geometric/agglomeration multigrid does not replace the fine-grid linear algebra with an unrelated solver. It provides a hierarchy that attacks error components on different spatial scales. Each active coarse level can reuse the same coupled residual/Jacobian machinery and its configured linear solver, with smoothing/coarse-solve effort chosen for the multigrid role. Chapter 18 describes this separation.

9.10 CMPS scalable-solver option map

The Advanced Linear Solver panel separates the outer family, local preconditioner, algebraic-multigrid relaxation/smoothing, cycle type and the outer flexible Krylov method used when AAMG acts as a preconditioner. The customer-facing GUI names are PCKSS, AAMG and AAMGKSS; their detailed controls and defaults are listed in the Linear Solver GUI chapter.

Role	Current choices
Scalable solver family	PCKSS; aggregation AMG; AMG-preconditioned Krylov
Local PCKSS preconditioner	Block Jacobi; block ILU
PCKSS Krylov method	GMRES; BiCGSTAB; TFQMR
AMG level relaxation	ILU; symmetric Gauss–Seidel; additive Schwarz
AMG level Krylov smoother	None; GMRES; BiCGSTAB; TFQMR
AMG-preconditioned outer Krylov	FGMRES; flexible BiCGSTAB
AMG cycle	V; F; W; FW

Stand-alone AAMG and the AAMG preconditioner have separate controls for level count, relaxation, smoothing, fill level, threshold, smoothing count and cycle. The current defaults use five levels, two smoothing applications, ILU level relaxation, no extra Krylov smoother and an F-cycle. AAMGKSS uses FGMRES by default. Treat these as starting values rather than universal optima; total time-to-solution is the relevant tuning metric.