General concepts » PETSc: Parallel Solvers and Distributed Linear Algebra

Introduction

Rodin's PETSc module integrates the Portable, Extensible Toolkit for Scientific Computation (PETSc) into the finite element pipeline. It provides:

• PETSc-backed trial/test functions, grid functions, bilinear forms, and linear forms that store data in PETSc Vec and Mat objects.
• KSP (Krylov subspace) and SNES (nonlinear) solver wrappers.
• Sequential, OpenMP, and MPI assembly strategies targeting PETSc objects.
• Seamless support for both Context::Local (single-process) and Context::MPI (distributed) mesh contexts.

Include the module via: #include <Rodin/PETSc.h>

PETSc Variational Components

The PETSc variational wrappers live in the Rodin::PETSc::Variational namespace. They are thin CTAD wrappers around the core Rodin classes that trigger PETSc-specific template deduction.

P1 Vh(mesh);

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

// Automatically selected via CTAD
BilinearForm a(u, v);  // assembles into ::Mat
a = Integral(Grad(u), Grad(v));

LinearForm L(v);       // assembles into ::Vec
L = Integral(f, v);

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

// Assembly + solve
Solver::CG(poisson).solve();

// The solution is automatically scattered to u.getSolution()

PETSc Solvers

// Direct usage
Solver::KSP(problem).solve();

// Configured usage
Solver::KSP ksp(problem);
ksp.setType(KSPGMRES);
ksp.setTolerances(
  1e-10,    // rtol:   relative decrease in residual norm
  1e-14,    // abstol: absolute residual norm threshold
  1e5,      // dtol:   divergence tolerance
  1000      // maxIt:  maximum iterations
);
ksp.setPreconditioner(precondMat);
ksp.solve();

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

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

Solver::KSP ksp(problem);
Solver::SNES snes(ksp);
snes.setType(SNESNEWTONLS);
snes.setTolerances(
  1e-8,     // abstol
  1e-10,    // rtol
  1e-6,     // stol (step norm)
  50,       // maxIt (nonlinear iterations)
  10000     // maxF (function evaluations)
);
snes.solve(x);

Linear System

PETSc::Math::LinearSystem couples the PETSc matrix, RHS vector, and solution vector:

auto& axb = problem.getLinearSystem();
auto& A = axb.getOperator();   // ::Mat
auto& b = axb.getVector();     // ::Vec (RHS)
auto& x = axb.getSolution();   // ::Vec (solution)

// The Problem class automatically sets up field splits for multi-field
// problems. Configure via command line:
//   -pc_type fieldsplit
//   -pc_fieldsplit_type schur
//   -fieldsplit_velocity_ksp_type cg
//   -fieldsplit_pressure_ksp_type preonly

Assembly Strategies

The PETSc module provides three assembly strategies:

The assembly strategy is selected automatically via the Assembly::Default<TrialContext, TestContext> trait based on the mesh context. Users do not normally need to select it explicitly.

Multi-Field Problems

PETSc problems support an arbitrary number of coupled trial/test function pairs. This is essential for saddle-point systems like Stokes, which involve velocity, pressure, and (optionally) a Lagrange multiplier for the mean-pressure constraint.

// FE spaces: H1(2) for velocity, H1(1) for pressure, P0g for LM
// The first argument to H1 is a compile-time order using std::integral_constant
H1  uh(std::integral_constant<size_t, 2>{}, mesh, mesh.getSpaceDimension());
H1  ph(std::integral_constant<size_t, 1>{}, mesh);
P0g lh(mesh);

PETSc::Variational::TrialFunction u(uh); u.setName("u");
PETSc::Variational::TrialFunction p(ph); p.setName("p");
PETSc::Variational::TrialFunction l(lh); l.setName("lambda");

PETSc::Variational::TestFunction v(uh);
PETSc::Variational::TestFunction q(ph);
PETSc::Variational::TestFunction m(lh);

Problem stokes(u, p, l, v, q, m);
stokes = Integral(Jacobian(u), Jacobian(v))
       - Integral(p, Div(v))
       + Integral(Div(u), q)
       + Integral(l, q)
       + Integral(p, m)
       + 1e-12 * Integral(p, q)
       + 1e-12 * Integral(l, m)
       - Integral(f, v)
       + DirichletBC(u, u_exact);

stokes.assemble();
Solver::KSP(stokes).solve();

Sequential PETSc Example

This example solves a 3D Poisson equation on hexahedra using PETSc CG:

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

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

int main(int argc, char** argv)
{
  PetscInitialize(&argc, &argv, PETSC_NULLPTR, PETSC_NULLPTR);

  Mesh mesh;
  mesh = mesh.UniformGrid(Polytope::Type::Hexahedron, {16, 16, 16});
  mesh.scale(1.0 / 15);
  mesh.getConnectivity().compute(2, 3);
  mesh.getConnectivity().compute(3, 2);

  RealFunction f = 1;
  P1 Vh(mesh);

  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();

  u.getSolution().save("Poisson.gf", IO::FileFormat::MFEM);
  mesh.save("Poisson.mesh", IO::FileFormat::MFEM);

  PetscFinalize();
}

MPI + PETSc Example

See the MPI guide for the full distributed workflow. The key difference from sequential usage is:

1. Create a Context::MPI and initialize PETSc.
2. Partition and distribute the mesh using Sharder.
3. Reconcile shared entities.
4. Use PETSc::Variational::TrialFunction / TestFunction as usual.
5. The assembly and solve are automatically distributed.

Command-Line Configuration

All PETSc solvers support command-line configuration via PETSc's options database. This is especially powerful for experimenting with preconditioners and solver parameters without recompiling:

See Also

• PETSc namespace reference
• KSP solver
• SNES nonlinear solver
• MPI module
• Solvers guide