General concepts » Solvers

Direct Solvers

Iterative Solvers

SparseLU

The SparseLU solver performs a sparse LU decomposition. It works for any non-singular matrix, making it the most general-purpose choice:

#include <Rodin/Solver.h>
#include <Rodin/Geometry.h>
#include <Rodin/Variational.h>

using namespace Rodin;
using namespace Rodin::Geometry;
using namespace Rodin::Variational;

int main()
{
  Mesh mesh;
  mesh = mesh.UniformGrid(Polytope::Type::Triangle, {16, 16});
  mesh.getConnectivity().compute(1, 2);

  P1 Vh(mesh);
  TrialFunction u(Vh);
  TestFunction  v(Vh);

  Problem poisson(u, v);
  poisson = Integral(Grad(u), Grad(v))
          - Integral(RealFunction(1.0), v)
          + DirichletBC(u, Zero());

  // Inline usage — construct solver with the problem, then solve
  Solver::SparseLU(poisson).solve();

  return 0;
}

SparseLU is a good default choice when:

• The system is not too large (up to roughly 100,000 DOFs)
• You need a robust solver that works for any matrix structure
• The system matrix is non-symmetric or indefinite

UMFPack

UMFPack is an alternative direct solver from the SuiteSparse library. It is especially effective for general unsymmetric sparse systems and often outperforms SparseLU on larger problems:

#include <Rodin/Solver/UMFPack.h>

Solver::UMFPack(problem).solve();

UMFPack requires the SuiteSparse library to be installed. It is the recommended direct solver for:

• Saddle-point systems (e.g. Stokes, Navier–Stokes)
• Coupled multi-physics problems
• Systems where SparseLU is too slow

Conjugate Gradient (CG)

The CG method is optimal for symmetric positive definite (SPD) systems, which arise from many common PDEs (Poisson, elasticity, heat equation):

// Inline usage — construct with the problem and solve immediately
Solver::CG(poisson).solve();

The CG solver supports several configuration options:

Solver::CG solver(poisson);
solver.setTolerance(1e-10);           // Relative convergence threshold
solver.setMaxIterations(1000);        // Maximum number of iterations
solver.solve();

if (solver.success())
  std::cout << "CG converged." << std::endl;

CG is the recommended solver when:

• The system matrix is symmetric positive definite
• The system is large (100,000+ DOFs)
• You want memory-efficient iterative solving

Hyperelastic Example

For hyperelastic problems, the Newton system is built from Solid::MaterialTangent (tangent stiffness) and Solid::InternalForce (nonlinear residual), both linearized about the current displacement:

#include <Rodin/Solid.h>
#include <Rodin/Solver/NewtonSolver.h>
#include <Rodin/Solver/SparseLU.h>

Solid::NeoHookean law(lambda, mu);
GridFunction u(Vh);   // accumulated displacement

Solid::MaterialTangent tangent(law, du, v);
tangent.setDisplacement(u);

Solid::InternalForce residual(law, v);
residual.setDisplacement(u);

Problem newton(du, v);
newton = tangent + residual
       - Integral(bodyForce, v)
       + DirichletBC(du, VectorFunction{Zero(), Zero()}).on(fixedBC);

SparseLU linearSolver(newton);
NewtonSolver solver(linearSolver);
solver.setMaxIterations(50)
      .setAbsoluteTolerance(1e-10);
solver.solve(u);

The Newton solver is used for:

• Hyperelastic solid mechanics (with Solid::InternalForce / MaterialTangent)
• Nonlinear PDEs (Navier–Stokes, nonlinear diffusion)
• Any problem where the bilinear form depends on the current solution

Inline Solving

For quick one-shot use, construct the solver with the problem and call solve() in a single expression. This is the most common pattern in Rodin:

Solver::CG(poisson).solve();
Solver::SparseLU(elasticity).solve();

Configured Solving

When you need to tune parameters, construct the solver object first, configure it, then call solve():

Solver::CG solver(poisson);
solver.setMaxIterations(500);
solver.setTolerance(1e-12);
solver.solve();

Solver Reuse via problem.solve()

You can also pass a pre-configured solver to problem.solve(). This is useful when the same solver is applied to multiple problems (e.g.\ in an optimization loop):

Solver::SparseLU solver;
for (size_t i = 0; i < nIterations; ++i)
{
  // ... update problem ...
  problem.solve(solver);
}

KSP (Krylov Subspace)

KSP is the base PETSc linear solver wrapper. It provides access to the full PETSc Krylov solver library and can be configured both programmatically and via PETSc command-line options.

#include <Rodin/PETSc.h>

// Use the KSP base solver directly
Solver::KSP(problem).solve();

// Configure via PETSc command-line options:
//   ./myprogram -ksp_type gmres -pc_type ilu -ksp_rtol 1e-8 -ksp_monitor

PETSc CG

CG<PETSc::Math::LinearSystem> is a convenience specialization that sets the KSP type to KSPCG. It is the recommended iterative solver for SPD systems when using PETSc:

PETSc::Variational::TrialFunction u(Vh);
PETSc::Variational::TestFunction  v(Vh);

Problem poisson(u, v);
poisson = Integral(Grad(u), Grad(v))
        - Integral(f, v)
        + DirichletBC(u, Zero());

Solver::CG(poisson).solve();

PETSc GMRES

GMRES<PETSc::Math::LinearSystem> sets the KSP type to KSPGMRES. Suitable for non-symmetric and indefinite systems:

Solver::GMRES(stokes).solve();

SNES (Nonlinear)

SNES wraps the PETSc Scalable Nonlinear Equations Solver for Newton-type methods. It uses an associated KSP for the inner linear solve at each Newton step:

// Create the KSP linear solver
Solver::KSP ksp(problem);
ksp.setTolerances(1e-10, 1e-14, 1e5, 200);

// Wrap it in SNES for Newton iteration
Solver::SNES snes(ksp);
snes.setTolerances(1e-8, 1e-10, 1e-6, 50, 10000);
snes.solve(x);

Multi-Field Problems

PETSc solvers support block systems with multiple trial/test function pairs. The resulting linear system can be solved with field-split preconditioning configured via PETSc command-line options:

PETSc::Variational::TrialFunction u(Vh_velocity);
PETSc::Variational::TrialFunction p(Vh_pressure);

PETSc::Variational::TestFunction v(Vh_velocity);
PETSc::Variational::TestFunction q(Vh_pressure);

Problem stokes(u, p, v, q);
stokes = Integral(Jacobian(u), Jacobian(v))
       - Integral(p, Div(v))
       + Integral(Div(u), q)
       - Integral(f, v)
       + DirichletBC(u, g);

Solver::KSP(stokes).solve();
// Run with: -ksp_type gmres -pc_type fieldsplit -ksp_monitor

PETSc Solver Summary