Refactor PETSc preconditioner for CVODE

CMakeLists.txt

@@ -355,6 +355,7 @@ set(BOUT_SOURCES
     ./src/solver/impls/snes/snes.hxx
     ./src/solver/impls/split-rk/split-rk.cxx
     ./src/solver/impls/split-rk/split-rk.hxx
+    ./src/solver/petsc_preconditioner.cxx
     ./src/solver/solver.cxx
     ./src/sys/adios_object.cxx
     ./src/sys/bout_types.cxx

include/bout/petsc_preconditioner.hxx

@@ -0,0 +1,94 @@
+/**************************************************************************
+ * Reusable PETSc coloring-based Jacobian preconditioner support
+ *
+ * This class factors out the duplicated code in PETSc-based solvers for:
+ * - building a Jacobian sparsity pattern (via stencil assumptions)
+ * - creating a MatFDColoring context for efficient finite-difference Jacobians
+ *
+ **************************************************************************/
+
+#ifndef BOUT_PETSC_PRECONDITIONER_H
+#define BOUT_PETSC_PRECONDITIONER_H
+
+#include "bout/build_defines.hxx"
+
+#if BOUT_HAS_PETSC
+
+#include "bout/petsc_interface.hxx"
+#include "bout/petsclib.hxx"
+
+#include <mpi.h>
+
+#include <petscmat.h>
+#include <petscvec.h>
+
+class Options;
+class Field3D;
+
+class PetscPreconditioner {
+public:
+  PetscPreconditioner() = default;
+  ~PetscPreconditioner() { reset(); }
+
+  PetscPreconditioner(const PetscPreconditioner&) = delete;
+  PetscPreconditioner& operator=(const PetscPreconditioner&) = delete;
+  PetscPreconditioner(PetscPreconditioner&&) = delete;
+  PetscPreconditioner& operator=(PetscPreconditioner&&) = delete;
+
+  /// Create and assemble a finite-difference Jacobian matrix `Jfd` with a sparsity
+  /// pattern derived from the solver's variable layout and a user-specified stencil.
+  ///
+  /// `index` should be generated by the calling solver (e.g. `Solver::globalIndex(0)`).
+  /// Safe to call multiple times; any existing PETSc objects are destroyed first.
+  PetscErrorCode createJacobianPattern(Field3D& index, Options& options, PetscInt nlocal,
+                                       int n2d, int n3d, MPI_Comm comm);
+
+  /// Recalculate the PETSc coloring and finite-difference coloring context
+  /// for the currently stored `Jfd` matrix.
+  template <class ColoringFunction>
+  PetscErrorCode updateColoring(ColoringFunction function, void* ctx) {
+    // Re-calculate the coloring
+    MatColoring coloring{nullptr};
+    PetscCall(MatColoringCreate(Jfd, &coloring));
+    PetscCall(MatColoringSetType(coloring, MATCOLORINGGREEDY));
+    PetscCall(MatColoringSetFromOptions(coloring));
+
+    // Calculate new index sets
+    ISColoring iscoloring{nullptr};
+    PetscCall(MatColoringApply(coloring, &iscoloring));
+    PetscCall(MatColoringDestroy(&coloring));
+
+    // Replace the old coloring with the new one
+    PetscCall(MatFDColoringDestroy(&fdcoloring));
+    PetscCall(MatFDColoringCreate(Jfd, iscoloring, &fdcoloring));
+    PetscCall(
+        MatFDColoringSetFunction(fdcoloring, bout::cast_MatFDColoringFn(function), ctx));
+    PetscCall(MatFDColoringSetFromOptions(fdcoloring));
+    PetscCall(MatFDColoringSetUp(Jfd, iscoloring, fdcoloring));
+    PetscCall(ISColoringDestroy(&iscoloring));
+
+    PetscFunctionReturn(PETSC_SUCCESS);
+  }
+
+  Mat jacobian() const { return Jfd; }
+  MatFDColoring coloring() const { return fdcoloring; }
+
+  void reset();
+
+private:
+  Mat Jfd{nullptr};
+  MatFDColoring fdcoloring{nullptr};
+};
+
+#else
+
+// PETSc is not available, but keep the type defined so headers can be included
+// unconditionally in PETSc-enabled compilation units.
+class PetscPreconditioner {
+public:
+  void reset() {}
+};
+
+#endif // BOUT_HAS_PETSC
+
+#endif // BOUT_PETSC_PRECONDITIONER_H

manual/sphinx/user_docs/time_integration.rst

@@ -33,7 +33,7 @@ needed to make the solver available.
 
 .. _tab-solvers:
 .. table:: Available time integration solvers
-	   
+
    +---------------+-----------------------------------------+------------------------+
    | Name          | Description                             | Compile options        |
    +===============+=========================================+========================+
@@ -68,7 +68,7 @@ given in table :numref:`tab-solveropts`.
 
 .. _tab-solveropts:
 .. table:: Time integration solver options
-	   
+
    +--------------------------+--------------------------------------------+-------------------------------------+
    | Option                   | Description                                | Solvers used                        |
    +==========================+============================================+=====================================+
@@ -144,6 +144,31 @@ iterations becomes large, this may be an indication that the system is
 poorly conditioned, and a preconditioner might help improve performance.
 See :ref:`sec-preconditioning`.
 
+CVODE preconditioning is controlled using ``solver:cvode_precon_method``:
+
+- ``none`` (default): Disable preconditioning.
+- ``auto``: Prefer a user-supplied preconditioner if provided, then PETSc
+  coloring if PETSc is available, otherwise use BBD.
+- ``user``: Require a user-supplied preconditioner.
+- ``petsc``: Require PETSc and use PETSc coloring.
+- ``bbd``: Force the built-in BBD preconditioner.
+
+For ``cvode_precon_method = petsc``, PETSc options for the internal KSP/PC can be
+set with the prefix ``cvode_petscpre_`` (either on the command line, or by putting
+prefixed keys into the ``[petsc]`` section). For example::
+
+    [petsc]
+    cvode_petscpre_ksp_type = preonly
+    cvode_petscpre_pc_type = hypre
+
+Two CVODE heuristics that control when the linear solver setup routine is called,
+and when the Jacobian/preconditioner are recomputed, can be adjusted with:
+
+- ``cvode_lsetup_frequency`` (default ``0``): Passed to ``CVodeSetLSetupFrequency``.
+  ``0`` uses the SUNDIALS default.
+- ``cvode_jac_eval_frequency`` (default ``0``): Passed to ``CVodeSetJacEvalFrequency``.
+  ``0`` uses the SUNDIALS default.
+
 CVODE can set constraints to keep some quantities positive, non-negative,
 negative or non-positive. These constraints can be activated by setting the
 option ``solver:apply_positivity_constraints=true``, and then in the section