Solver Toolkit

PETSc

Portable, Extensible Toolkit for Scientific Computation. Provides scalable linear/nonlinear solvers, time integrators, and optimisation algorithms for PDEs on distributed-memory systems. The backbone solver library used by FEniCS, deal.II, MOOSE, and many others.

C / C++ / Fortran / Python

MPI Distributed

GPU (CUDA/HIP/SYCL)

Open Source (BSD)

Solver Hierarchy

PETSc uses a composable, hierarchical solver architecture. Each level can be configured at runtime via command-line options:

ApplicationPhysics-specific code

TS / TAOTime integration / Optimisation

SNESNonlinear solvers

KSPKrylov linear solvers

PCPreconditioners

Mat / VecMatrices & Vectors

Runtime Configurability

One of PETSc's most powerful features — change solver algorithms without recompiling:

# Run with CG solver + algebraic multigrid preconditioner
mpirun -np 16 ./myapp -ksp_type cg -pc_type gamg

# Switch to GMRES + ILU(3) at runtime
mpirun -np 16 ./myapp -ksp_type gmres -pc_type ilu -pc_factor_levels 3

# Enable convergence monitoring
mpirun -np 16 ./myapp -ksp_monitor -ksp_converged_reason

Convergence Criterion

KSP solvers converge when the relative residual drops below a tolerance:

Default rtol = 10⁻⁵. Configurable via -ksp_rtol.

Core Components

Linear Solvers (KSP)

Krylov subspace methods: GMRES, CG, BiCGStab, MINRES, TFQMR, and more. Pluggable preconditioner interface supports Jacobi, SOR, ILU, ICC, ASM, multigrid (PCMG), and algebraic multigrid (GAMG, BoomerAMG via Hypre).

Nonlinear Solvers (SNES)

Newton methods (line search, trust region), nonlinear CG, Anderson mixing, NGMRES, and quasi-Newton (L-BFGS). Jacobian-free Newton-Krylov (JFNK) for matrix-free applications.

Time Steppers (TS)

Explicit (Euler, RK), implicit (backward Euler, Crank-Nicolson, BDF), IMEX (ARK), and adaptive time-stepping with error control. Adjoint and tangent linear model support for sensitivity analysis.

Distributed Arrays (DM)

Structured (DMDA) and unstructured (DMPlex) mesh management with built-in partitioning. Automatic parallel data layout, ghost value exchange, and interpolation between mesh levels for multigrid.

Example: Linear Solve

solve.c

// PETSc: Solve Ax = b with CG + Jacobi preconditioner
#include <petsc.h>

int main(int argc, char **argv) {
  PetscInitialize(&argc, &argv, NULL, NULL);
  
  Mat A; Vec x, b;
  KSP ksp; PC pc;
  PetscInt n = 1000;

  // Create matrix and vectors
  MatCreate(PETSC_COMM_WORLD, &A);
  MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
  MatSetFromOptions(A);
  MatSetUp(A);
  
  // ... assemble tridiagonal system ...

  VecCreate(PETSC_COMM_WORLD, &b);
  VecSetSizes(b, PETSC_DECIDE, n);
  VecSetFromOptions(b);
  VecDuplicate(b, &x);

  // Create solver
  KSPCreate(PETSC_COMM_WORLD, &ksp);
  KSPSetOperators(ksp, A, A);
  KSPSetType(ksp, KSPCG);
  
  // Set preconditioner
  KSPGetPC(ksp, &pc);
  PCSetType(pc, PCJACOBI);
  
  // Solve
  KSPSolve(ksp, b, x);

  // Clean up
  KSPDestroy(&ksp);
  MatDestroy(&A); VecDestroy(&x); VecDestroy(&b);
  PetscFinalize();
  return 0;
}

Resources: Visit petsc.org for documentation, tutorials, and the extensive manual. VELSTROM integration tutorials for PETSc-based earth science solvers are under development.