Snes with constrains updated

include/bout/field2d.hxx

@@ -601,6 +601,9 @@ inline Field2D DC(const Field2D& f) { return f; }
 /// Wrapper around member function f.timeDeriv()
 inline Field2D& ddt(Field2D& f) { return *(f.timeDeriv()); }
 
+/// Alias for ddt that is more readable for constraints
+inline Field2D& residual(Field2D& f) { return *(f.timeDeriv()); }
+
 /// toString template specialisation
 /// Defined in utils.hxx
 template <>

include/bout/field3d.hxx

@@ -995,6 +995,9 @@ inline void invalidateGuards([[maybe_unused]] Field3D& var) {}
 /// Wrapper around member function f.timeDeriv()
 inline Field3D& ddt(Field3D& f) { return *(f.timeDeriv()); }
 
+/// Alias for ddt that is more readable for constraints
+inline Field3D& residual(Field3D& f) { return *(f.timeDeriv()); }
+
 /// toString template specialisation
 /// Defined in utils.hxx
 template <>

include/bout/solver.hxx

@@ -258,8 +258,8 @@ public:
   virtual void add(Vector3D& v, const std::string& name,
                    const std::string& description = "");
 
-  /// Returns true if constraints available
-  virtual bool constraints() { return has_constraints; }
+  /// Returns true if this solver supports constraint variables
+  virtual bool constraints() { return supports_constraints; }
 
   /// Add constraint functions (optional). These link a variable v to
   /// a control parameter C_v such that v is adjusted to keep C_v = 0.
@@ -549,7 +549,7 @@ protected:
   };
 
   /// Can this solver handle constraints? Set to true if so.
-  bool has_constraints{false};
+  bool supports_constraints{false};
   /// Has init been called yet?
   bool initialised{false};
   /// If calling user RHS for the first time

manual/sphinx/user_docs/time_integration.rst

@@ -524,8 +524,9 @@ Timestepping Modes
 
 The solver supports several timestepping strategies controlled by ``equation_form``:
 
-**Backward Euler (default)**
-   Standard implicit backward Euler method. Good for general timestepping.
+**Rearranged Backward Euler (default)**
+   Standard implicit backward Euler method written in a rearranged form that is
+   robust when driving a system to steady state.
 
    .. code-block:: ini
 
@@ -555,6 +556,65 @@ The solver supports several timestepping strategies controlled by ``equation_for
    This uses the same form as rearranged_backward_euler, but the time step
    can be different for each cell.
 
+Constraints (DAEs)
+~~~~~~~~~~~~~~~~~~
+
+BOUT++ can define algebraic constraints in a physics model using the
+``Solver::constraint(...)`` API. With the SNES solver these are treated as a
+differential-algebraic equation (DAE) system:
+
+- Differential variables include the usual timestepping terms for the selected
+  ``equation_form``.
+- Algebraic variables keep only the constraint residual ``G(x) = 0``.
+
+This is supported for all SNES equation forms. For ``direct_newton`` the full
+system is solved directly as a steady-state nonlinear problem, while for the
+timestepping forms only the differential variables receive the timestep terms.
+
+When constraints are enabled, the SNES solver can optionally split the
+preconditioner into differential and algebraic blocks using PETSc
+``fieldsplit``. The split names are:
+
+- ``diff``: differential variables
+- ``alg``: algebraic (constraint) variables
+
+Constraint splitting requires ``matrix_free = false``. ``matrix_free_operator``
+may still be used because the preconditioner matrix is still assembled.
+
+Example:
+
+.. code-block:: ini
+
+   [solver]
+   type = snes
+   equation_form = backward_euler
+   pc_type = fieldsplit
+
+   [petsc]
+   pc_fieldsplit_type = additive
+   fieldsplit_diff_ksp_type = preonly
+   fieldsplit_diff_pc_type = jacobi
+   fieldsplit_alg_ksp_type = preonly
+   fieldsplit_alg_pc_type = jacobi
+
+A more tailored setup can use different preconditioners for the two blocks, for
+example ILU on the differential variables and Hypre BoomerAMG on the algebraic
+variables:
+
+.. code-block:: ini
+
+   [solver]
+   type = snes
+   equation_form = backward_euler
+   matrix_free = false
+   pc_type = fieldsplit
+
+   [petsc]
+   pc_fieldsplit_type = additive
+   fieldsplit_diff_pc_type = ilu
+   fieldsplit_alg_pc_type = hypre
+   fieldsplit_alg_pc_hypre_type = boomeramg
+
 Adaptive Timestepping
 ~~~~~~~~~~~~~~~~~~~~~
 
@@ -1422,19 +1482,20 @@ similar way to time integrated variables. For example
 
     Field3D phi;
     ...
-    solver->constraint(phi, ddt(phi), "phi");
+    solver->constraint(phi, residual(phi), "phi");
 
 The first argument is the variable to be solved for (constrained). The
 second argument is the field to contain the residual (error). In this
-example the time derivative field ``ddt(phi)`` is used, but it could
-be another `Field3D` variable. The solver will attempt to
+example ``residual(phi)`` is used, which is an alias for the same
+storage as ``ddt(phi)`` but makes the algebraic role clearer. It could
+also be another ``Field3D`` variable. The solver will attempt to
 find a solution to the first argument (``phi`` here) such that the
-second argument (``ddt(phi)``) is zero to within tolerances.
+second argument (``residual(phi)`` here) is zero to within tolerances.
 
 In the RHS function the residual should be calculated. In this example
 (``examples/constraints/drift-wave-constraint``) we have::
 
-    ddt(phi) = Delp2(phi) - Vort;
+    residual(phi) = Delp2(phi) - Vort;
 
 so the time integration solver includes the algebraic constraint
 ``Delp2(phi) = Vort`` i.e. (:math:`\nabla_\perp^2\phi = \omega`).

src/solver/impls/arkode/arkode.cxx

@@ -162,7 +162,7 @@ ArkodeSolver::ArkodeSolver(Options* opts)
           (*options)["optimize"].doc("Use ARKode optimal parameters").withDefault(false)),
 #endif
       suncontext(createSUNContext(BoutComm::get())) {
-  has_constraints = false; // This solver doesn't have constraints
+  supports_constraints = false; // This solver doesn't have constraints
 
   // Add diagnostics to output
   add_int_diagnostic(nsteps, "arkode_nsteps", "Cumulative number of internal steps");

src/solver/impls/cvode/cvode.cxx

@@ -181,7 +181,7 @@ CvodeSolver::CvodeSolver(Options* opts)
               .doc("When to save CVODE Jacobians: output or linear_setup")
               .withDefault(bout::CvodeJacobianExportTrigger::linear_setup)),
       suncontext(createSUNContext(BoutComm::get())) {
-  has_constraints = false; // This solver doesn't have constraints
+  supports_constraints = false; // This solver doesn't have constraints
   canReset = true;
 
   if ((*options)["use_precon"].isSet()) {

src/solver/impls/ida/ida.cxx

@@ -85,7 +85,7 @@ IdaSolver::IdaSolver(Options* opts)
                        .doc("N_Vector backend to use: sundials or manyvector")
                        .withDefault(NVectorType::Sundials)),
       suncontext(createSUNContext(BoutComm::get())) {
-  has_constraints = true; // This solver has constraints
+  supports_constraints = true; // This solver has constraints
 }
 
 IdaSolver::~IdaSolver() {

src/solver/impls/imex-bdf2/imex-bdf2.cxx

@@ -84,7 +84,7 @@ IMEXBDF2::IMEXBDF2(Options* opt)
       verbose((*options)["verbose"]
                   .doc(" More outputs at each timestep")
                   .withDefault(false)) {
-  has_constraints = true; ///< This solver can handle constraints
+  supports_constraints = true; ///< This solver can handle constraints
 }
 
 IMEXBDF2::~IMEXBDF2() {
@@ -154,22 +154,22 @@ int IMEXBDF2::init() {
                n2Dvars(), neq, nlocal);
 
   // Check if there are any constraints
-  have_constraints = false;
+  has_constraint_variables = false;
 
   for (int i = 0; i < n2Dvars(); i++) {
     if (f2d[i].constraint) {
-      have_constraints = true;
+      has_constraint_variables = true;
       break;
     }
   }
   for (int i = 0; i < n3Dvars(); i++) {
     if (f3d[i].constraint) {
-      have_constraints = true;
+      has_constraint_variables = true;
       break;
     }
   }
 
-  if (have_constraints) {
+  if (has_constraint_variables) {
     is_dae.reallocate(nlocal);
     // Call the Solver function, which sets the array
     // to zero when not a constraint, one for constraint
@@ -1273,7 +1273,7 @@ PetscErrorCode IMEXBDF2::snes_function(Vec x, Vec f, bool linear) {
 
   // G(x) now in fdata
 
-  if (!have_constraints) {
+  if (!has_constraint_variables) {
     // No constraints, so simple loop over all variables
 
     for (int i = 0; i < nlocal; i++) {

src/solver/impls/imex-bdf2/imex-bdf2.hxx

@@ -200,8 +200,8 @@ private:
   int linear_fails;    ///< Number of linear (KSP) convergence failures
   int nonlinear_fails; ///< Numbef of nonlinear (SNES) convergence failures
 
-  bool have_constraints;  ///< Are there any constraint variables?
-  Array<BoutReal> is_dae; ///< If using constraints, 1 -> DAE, 0 -> AE
+  bool has_constraint_variables; ///< Are there any constraint variables?
+  Array<BoutReal> is_dae;        ///< If using constraints, 1 -> DAE, 0 -> AE
 
   MatFDColoring fdcoloring; ///< Matrix coloring context, used for finite difference
                             ///< Jacobian evaluation

src/solver/impls/pvode/pvode.cxx

@@ -101,7 +101,7 @@ PvodeSolver::PvodeSolver(Options* opts)
           (*options)["mxstep"].doc("Maximum number of internal steps").withDefault(500)),
       abstol((*options)["atol"].doc("Absolute tolerance").withDefault(1.0e-12)),
       reltol((*options)["rtol"].doc("Relative tolerance").withDefault(1.0e-5)) {
-  has_constraints = false; ///< This solver doesn't have constraints
+  supports_constraints = false; ///< This solver doesn't have constraints
 }
 
 PvodeSolver::~PvodeSolver() {

src/solver/impls/slepc/slepc.cxx

@@ -148,7 +148,7 @@ std::string formatEig(BoutReal reEig, BoutReal imEig) {
 }
 
 SlepcSolver::SlepcSolver(Options* options) {
-  has_constraints = false;
+  supports_constraints = false;
   initialised = false;
   stIsShell = PETSC_FALSE;