Generic MATRIX-FREE Krylov solver loops – richardson_solve, cg_solve, bicgstab_solve, gmres_solve – operating on pops::MultiFab with the operator supplied as a CALLBACK (ApplyFn). More...

#include <pops/core/foundation/types.hpp>
#include <pops/mesh/storage/mf_arith.hpp>
#include <pops/mesh/storage/multifab.hpp>
#include <pops/numerics/elliptic/linear/krylov_result.hpp>
#include <cmath>
#include <functional>
#include <stdexcept>
#include <vector>

Include dependency graph for generic_krylov.hpp:

Go to the source code of this file.

Namespaces
namespace	pops

namespace	pops::detail

Typedefs
using	pops::ApplyFn = std::function< void(MultiFab &out, const MultiFab &in)>
	Matrix-free operator callback: `out <- A(in)`.

Functions
Real	pops::detail::krylov_dot (const MultiFab &x, const MultiFab &y)
	Krylov inner product x.y, COLLECTIVE (all_reduce_sum).

Real	pops::detail::krylov_l2_norm (const MultiFab &x)
	GLOBAL L2 norm sqrt(sum x.x), collective (all_reduce_sum).

void	pops::detail::require_max_iter (int max_iters)
	Guards a dynamic solver loop: a non-positive iteration budget is a configuration error.

KrylovResult	pops::richardson_solve (const ApplyFn &A, MultiFab &phi, const MultiFab &rhs, Real omega, Real rel_tol, int max_iters)
	Richardson iteration x <- x + omega (b - A x), solving A x = b.

KrylovResult	pops::cg_solve (const ApplyFn &A, MultiFab &phi, const MultiFab &rhs, Real rel_tol, int max_iters)
	Conjugate Gradient, solving A x = b for an SPD operator A.

KrylovResult	pops::bicgstab_solve (const ApplyFn &A, const ApplyFn &precond, MultiFab &phi, const MultiFab &rhs, Real rel_tol, int max_iters)
	Preconditioned BiCGStab, solving A x = b for a general (possibly non-symmetric) operator A.

KrylovResult	pops::gmres_solve (const ApplyFn &A, const ApplyFn &precond, MultiFab &phi, const MultiFab &rhs, Real rel_tol, int max_iters, int restart=30)
	Left-preconditioned restarted GMRES(m), solving A x = b for a GENERAL (possibly NON-symmetric) operator A.

Variables
constexpr Real	pops::detail::kKrylovTiny = Real(1e-300)
	Tiny breakdown guard for the BiCGStab scalar recurrences (division by ~0).

Detailed Description

Generic MATRIX-FREE Krylov solver loops – richardson_solve, cg_solve, bicgstab_solve, gmres_solve – operating on pops::MultiFab with the operator supplied as a CALLBACK (ApplyFn).

Layer: include/pops/numerics/elliptic/linear. Role: the reusable Krylov core. Where TensorKrylovSolver (krylov_solver.hpp) hardwires the operator and preconditioner to a GeometricMG (the 5-point Laplacian matvec + an MG V-cycle preconditioner on the symmetric part), this header LIFTS the operator/preconditioner to generic std::function callbacks: the caller supplies ANY matrix-free apply out <- A(in). The BiCGStab math and its fixed scratch structure are COPIED from TensorKrylovSolver; only the matvec / preconditioner indirection differs. Richardson and CG are added alongside. These loops are the dynamic solver primitives that a later slice wires into the compiled time-program (ProgramContext / codegen); this header is pure C++ numerics and has no Python / program dependency.

Convention: each loop solves A x = b (A = the A callback). phi is the unknown IN/OUT: the incoming value is the initial guess (warm start), the outgoing value is the solution. rhs is the right-hand side, never modified. Convergence is the RELATIVE L2 residual ||r|| <= rel_tol * ||b|| (||.|| = sqrt of the Krylov inner product, a GLOBAL L2 norm), or max_iters reached (returns converged=false, best effort). When ||b|| == 0 the base is taken as 1.

MULTI-COMPONENT (vector / state) FIELDS: every inner product goes through detail::krylov_dot, which reduces over ALL components (pops::dot_all) when the field has ncomp>1 and over component 0 only (pops::dot, BIT-IDENTICAL) when it is scalar. A vector solve thus drives the residual / search- direction norms and the CG / BiCGStab scalar recurrences over every component – a component-0-only dot would converge on component 0 alone and leave the others unsolved. The pointwise updates (saxpy / lincomb) already span all components, so no other change is needed; the scalar (ncomp==1) path is unchanged.

COLLECTIVE / ALL-RANKS CONTRACT: every reduction goes through krylov_dot -> pops::dot / dot_all, which perform a COLLECTIVE all_reduce_sum and MUST run on EVERY rank (including a rank with no box), otherwise MPI deadlocks. Because the residual driving the stopping test is collective, all ranks observe the SAME residual and break at the SAME iteration: the loops NEVER short-circuit a dot() and the trip count is identical on every rank. CG and BiCGStab use the same all_reduce'd dot, so they are also rank-divergence free.

ALLOC-ONCE / REUSE DISCIPLINE: every scratch MultiFab is allocated ONCE at the START of each call, co-distributed with phi (its BoxArray / DistributionMapping / ncomp / nghost), and reused across the iteration loop – there is NO heap allocation inside the loop. Richardson needs r (1 scratch); CG needs r, p, Ap (3); BiCGStab needs r, rhat, p, v, s, t, plus phat/shat when a preconditioner is supplied (mirrors TensorKrylovSolver's fixed footprint).

Constraints: max_iters <= 0 throws std::invalid_argument (a dynamic loop with no iteration budget is a configuration error; spec error 13). std::function lives only at the HOST loop level (never inside a device kernel), so it does not violate the device-clean kernel rules.

include/pops/numerics/elliptic/linear/generic_krylov.hpp File Reference¶

Namespaces

Typedefs

Functions

Variables

Detailed Description