diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx
index 13a26faa38..b522c593ab 100644
--- a/include/bout/field2d.hxx
+++ b/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 <>
diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 62486dded9..bf0a666771 100644
--- a/include/bout/field3d.hxx
+++ b/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 <>
diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx
index 70ab05a56a..f0f277ab18 100644
--- a/include/bout/solver.hxx
+++ b/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
diff --git a/manual/sphinx/user_docs/time_integration.rst b/manual/sphinx/user_docs/time_integration.rst
index cf0d08112c..ff3d5942ae 100644
--- a/manual/sphinx/user_docs/time_integration.rst
+++ b/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`).
diff --git a/src/solver/impls/arkode/arkode.cxx b/src/solver/impls/arkode/arkode.cxx
index 2035824054..c5dab45dd9 100644
--- a/src/solver/impls/arkode/arkode.cxx
+++ b/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");
diff --git a/src/solver/impls/cvode/cvode.cxx b/src/solver/impls/cvode/cvode.cxx
index 2d572087ca..38a27bc998 100644
--- a/src/solver/impls/cvode/cvode.cxx
+++ b/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()) {
diff --git a/src/solver/impls/ida/ida.cxx b/src/solver/impls/ida/ida.cxx
index 00e91f17c1..8a28a72a59 100644
--- a/src/solver/impls/ida/ida.cxx
+++ b/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() {
diff --git a/src/solver/impls/imex-bdf2/imex-bdf2.cxx b/src/solver/impls/imex-bdf2/imex-bdf2.cxx
index 24b4d028f9..26721da95b 100644
--- a/src/solver/impls/imex-bdf2/imex-bdf2.cxx
+++ b/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++) {
diff --git a/src/solver/impls/imex-bdf2/imex-bdf2.hxx b/src/solver/impls/imex-bdf2/imex-bdf2.hxx
index 9435695a2d..1d6c3f2824 100644
--- a/src/solver/impls/imex-bdf2/imex-bdf2.hxx
+++ b/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
diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx
index ba4dce5849..045d83ed49 100644
--- a/src/solver/impls/pvode/pvode.cxx
+++ b/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() {
diff --git a/src/solver/impls/slepc/slepc.cxx b/src/solver/impls/slepc/slepc.cxx
index 9887f184eb..63e3272bb6 100644
--- a/src/solver/impls/slepc/slepc.cxx
+++ b/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;
 
diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index 543bacb841..eb786822a3 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -20,10 +20,12 @@
 #include <bout/utils.hxx>
 
 #include <algorithm>
+#include <array>
 #include <cmath>
 #include <cstddef>
 #include <fmt/format.h>
 #include <string>
+#include <tuple>
 #include <vector>
 
 #include "petscerror.h"
@@ -73,6 +75,45 @@ PetscErrorCode FormScaledFunctionForColoring(void* UNUSED(snes), Vec x, Vec f,
   return static_cast<SNESSolver*>(ctx)->scaled_rhs_function(x, f, true);
 }
 
+/// Apply a PETSc vector operation either to whole vectors or to matching
+/// constrained subvectors selected by ``indices``.
+///
+/// This keeps the constrained and unconstrained paths in ``snes_function()``
+/// identical while ensuring any acquired subvectors are always restored.
+template <typename Func, typename... Args>
+PetscErrorCode withOptionalSubvectors(Func operation, IS indices, Args... args) {
+  if (indices == nullptr) {
+    return operation(args...);
+  }
+
+  constexpr std::size_t N = sizeof...(Args);
+  std::array<Vec, N> vectors{args...};
+  std::array<Vec, N> subvectors{};
+
+  PetscErrorCode ierr = PETSC_SUCCESS;
+  std::size_t acquired = 0;
+  for (; acquired < N; ++acquired) {
+    ierr = VecGetSubVector(vectors[acquired], indices, &subvectors[acquired]);
+    if (ierr != PETSC_SUCCESS) {
+      break;
+    }
+  }
+
+  if (ierr == PETSC_SUCCESS) {
+    ierr = std::apply(operation, subvectors);
+  }
+
+  for (std::size_t i = acquired; i > 0; --i) {
+    PetscErrorCode restore_ierr =
+        VecRestoreSubVector(vectors[i - 1], indices, &subvectors[i - 1]);
+    if (ierr == PETSC_SUCCESS) {
+      ierr = restore_ierr;
+    }
+  }
+
+  return ierr;
+}
+
 PetscErrorCode snesPCapply(PC pc, Vec x, Vec y) {
   // Get the context
   SNESSolver* s;
@@ -384,7 +425,18 @@ SNESSolver::SNESSolver(Options* opts)
                         .withDefault<bool>(false)),
       jacobian_export_kind((*options)["jacobian_export_kind"]
                                .doc("Which Jacobian to save: system, scaled, or rhs")
-                               .withDefault(bout::JacobianExportKind::system)) {}
+                               .withDefault(bout::JacobianExportKind::system)) {
+  supports_constraints = true; // This solver can handle constraints
+}
+
+SNESSolver::~SNESSolver() {
+  if (is_diff != nullptr) {
+    ISDestroy(&is_diff);
+  }
+  if (is_alg != nullptr) {
+    ISDestroy(&is_alg);
+  }
+}
 
 void SNESSolver::exportMatrixAndMetadata(bout::JacobianExportKind kind, Mat jacobian) {
   Solver::writeOnceJacobianMetadata("snes");
@@ -457,6 +509,28 @@ int SNESSolver::init() {
   output_info.write("\t3d fields = {:d}, 2d fields = {:d} neq={:d}, local_N={:d}\n",
                     n3Dvars(), n2Dvars(), neq, nlocal);
 
+  // Check if there are any constraints
+
+  for (int i = 0; i < n2Dvars(); i++) {
+    if (f2d[i].constraint) {
+      has_constraint_variables = true;
+      break;
+    }
+  }
+  for (int i = 0; i < n3Dvars(); i++) {
+    if (f3d[i].constraint) {
+      has_constraint_variables = true;
+      break;
+    }
+  }
+
+  if (has_constraint_variables) {
+    is_dae.reallocate(nlocal);
+    // Call the Solver function, which sets the array
+    // to one when not a constraint, zero for constraint
+    set_id(std::begin(is_dae));
+  }
+
   // Initialise fields for storing residual of nonlinear solves
   if (diagnose) {
     for (const auto& f : f2d) {
@@ -521,6 +595,32 @@ int SNESSolver::init() {
   local_residual_2d = 0.0;
   global_residual = 0.0;
 
+  if (has_constraint_variables) {
+    // CreatePETSc-native index sets representing the two parts of your DAE.
+    PetscInt istart, iend;
+    PetscCall(VecGetOwnershipRange(snes_x, &istart, &iend));
+    ASSERT2(iend - istart == nlocal);
+
+    std::vector<PetscInt> diff_idx;
+    std::vector<PetscInt> alg_idx;
+    diff_idx.reserve(nlocal);
+    alg_idx.reserve(nlocal);
+
+    for (PetscInt i = 0; i < nlocal; ++i) {
+      const PetscInt gi = istart + i;
+      if (is_dae[i] > 0.5) { // differential
+        diff_idx.push_back(gi);
+      } else { // algebraic constraint (i.e. phi)
+        alg_idx.push_back(gi);
+      }
+    }
+
+    PetscCall(ISCreateGeneral(BoutComm::get(), diff_idx.size(), diff_idx.data(),
+                              PETSC_COPY_VALUES, &is_diff));
+    PetscCall(ISCreateGeneral(BoutComm::get(), alg_idx.size(), alg_idx.data(),
+                              PETSC_COPY_VALUES, &is_alg));
+  }
+
   // Nonlinear solver interface (SNES)
   output_info.write("Create SNES\n");
   SNESCreate(BoutComm::get(), &snes);
@@ -600,6 +700,9 @@ int SNESSolver::init() {
   SNESSetForceIteration(snes, PETSC_TRUE);
 #endif
 
+  // Enable checking for domain errors in Jacobian evaluation
+  SNESSetCheckJacobianDomainError(snes, PETSC_TRUE);
+
   // Get KSP context from SNES
   KSP ksp;
   SNESGetKSP(snes, &ksp);
@@ -652,6 +755,17 @@ int SNESSolver::init() {
     }
   }
 
+  if (has_constraint_variables && !matrix_free && pc_type == "fieldsplit") {
+    output_info.write("Using PCFieldSplit preconditioner for DAE system\n");
+
+    // Use PETSc fieldsplit
+    PetscCall(PCSetType(pc, PCFIELDSPLIT));
+
+    // Give PETSc the index sets
+    PetscCall(PCFieldSplitSetIS(pc, "diff", is_diff));
+    PetscCall(PCFieldSplitSetIS(pc, "alg", is_alg));
+  }
+
   // Get runtime options
   lib.setOptionsFromInputFile(snes);
 
@@ -1577,34 +1691,69 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) {
   // Call the RHS function
   if (scaled_rhs_function(x, f, linear) != PETSC_SUCCESS) {
     // Tell SNES that the input was out of domain
-    SNESSetFunctionDomainError(snes);
+    if (linear) {
+      // During Jacobian evaluation
+      SNESSetJacobianDomainError(snes);
+    } else {
+      // During function evaluation
+      SNESSetFunctionDomainError(snes);
+    }
     // Note: Returning non-zero error here leaves vectors in locked state
     return 0;
   }
 
+  ASSERT2(!has_constraint_variables || is_diff != nullptr);
+
   switch (equation_form) {
   case BoutSnesEquationForm::rearranged_backward_euler: {
     // Rearranged Backward Euler
-    // f = (x0 - x)/Δt + f
-    // First calculate x - x0 to minimise floating point issues
-    VecWAXPY(delta_x, -1.0, x0, x); // delta_x = x - x0
-    VecAXPY(f, -1. / dt, delta_x);  // f <- f - delta_x / dt
+    // F = (x0 - x)/Δt + f
+    // Algebraic:     F = G(x)  (already stored in f by rhs_function)
+
+    // Apply the transient term to all variables, or only the differential part
+    // when constraints are present.
+    PetscCall(withOptionalSubvectors(
+        [this](Vec x_part, Vec x0_part, Vec delta_x_part, Vec f_part) -> PetscErrorCode {
+          PetscCall(VecWAXPY(delta_x_part, -1.0, x0_part, x_part));
+          PetscCall(VecAXPY(f_part, -1.0 / dt, delta_x_part));
+          return PETSC_SUCCESS;
+        },
+        has_constraint_variables ? is_diff : nullptr, x, x0, delta_x, f));
     break;
   }
   case BoutSnesEquationForm::pseudo_transient: {
     // Pseudo-transient timestepping. Same as Rearranged Backward Euler
     // except that Δt is a vector
-    // f = (x0 - x)/Δt + f
-    VecWAXPY(delta_x, -1.0, x0, x);
-    VecPointwiseDivide(delta_x, delta_x, dt_vec); // delta_x /= dt
-    VecAXPY(f, -1., delta_x);                     // f <- f - delta_x
+    // F = (x0 - x)/Δt + f
+    // Algebraic:     F = G(x)  (already stored in f by rhs_function)
+
+    // As above, only the differential variables get the pseudo-time update in
+    // constrained runs.
+    PetscCall(withOptionalSubvectors(
+        [](Vec x_part, Vec x0_part, Vec delta_x_part, Vec f_part,
+           Vec dt_vec_part) -> PetscErrorCode {
+          PetscCall(VecWAXPY(delta_x_part, -1.0, x0_part, x_part));
+          PetscCall(VecPointwiseDivide(delta_x_part, delta_x_part, dt_vec_part));
+          PetscCall(VecAXPY(f_part, -1.0, delta_x_part));
+          return PETSC_SUCCESS;
+        },
+        has_constraint_variables ? is_diff : nullptr, x, x0, delta_x, f, dt_vec));
     break;
   }
   case BoutSnesEquationForm::backward_euler: {
-    // Backward Euler
-    // Set f = x - x0 - Δt*f
-    VecAYPX(f, -dt, x);   // f <- x - Δt*f
-    VecAXPY(f, -1.0, x0); // f <- f - x0
+    // Backward Euler:
+    // Differential:  F = x - x0 - dt*f
+    // Algebraic:     F = G(x)  (already stored in f by rhs_function)
+
+    // Only differential variables receive the backward-Euler update when DAE
+    // constraints are active.
+    PetscCall(withOptionalSubvectors(
+        [this](Vec x_part, Vec x0_part, Vec f_part) -> PetscErrorCode {
+          PetscCall(VecAYPX(f_part, -dt, x_part));
+          PetscCall(VecAXPY(f_part, -1.0, x0_part));
+          return PETSC_SUCCESS;
+        },
+        has_constraint_variables ? is_diff : nullptr, x, x0, f));
     break;
   }
   case BoutSnesEquationForm::direct_newton: {
diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx
index ddb25b9975..57a4c4e842 100644
--- a/src/solver/impls/snes/snes.hxx
+++ b/src/solver/impls/snes/snes.hxx
@@ -40,6 +40,7 @@ class SNESSolver;
 
 #include "mpi.h"
 
+#include <bout/array.hxx>
 #include <bout/bout_enum_class.hxx>
 #include <bout/bout_types.hxx>
 #include <bout/field2d.hxx>
@@ -78,7 +79,7 @@ BOUT_ENUM_CLASS(BoutSnesOutput,
 class SNESSolver : public Solver {
 public:
   explicit SNESSolver(Options* opts = nullptr);
-  ~SNESSolver() override = default;
+  ~SNESSolver() override;
 
   int init() override;
   int run() override;
@@ -247,6 +248,12 @@ private:
   int nlocal; ///< Number of variables on local processor
   int neq;    ///< Number of variables in total
 
+  bool has_constraint_variables{false}; ///< Are there any constraint variables?
+  Array<BoutReal> is_dae;               ///< If using constraints, 1 -> DAE, 0 -> AE
+
+  IS is_diff = nullptr; // is_dae == 1
+  IS is_alg = nullptr;  // is_dae == 0 (phi constraint and any other algebraics)
+
   PetscLib lib; ///< Handles initialising, finalising PETSc
   Vec snes_f;   ///< Used by SNES to store function
   Vec deriv; ///< Time derivative; only used if diagnose = true, otherwise will store in snes_f
diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx
index 8388759b74..29b50cd915 100644
--- a/src/solver/solver.cxx
+++ b/src/solver/solver.cxx
@@ -420,7 +420,7 @@ void Solver::constraint(Field2D& v, Field2D& C_v, std::string name) {
   }
 #endif
 
-  if (!has_constraints) {
+  if (!supports_constraints) {
     throw BoutException("ERROR: This solver doesn't support constraints\n");
   }
 
@@ -451,7 +451,7 @@ void Solver::constraint(Field3D& v, Field3D& C_v, std::string name) {
   }
 #endif
 
-  if (!has_constraints) {
+  if (!supports_constraints) {
     throw BoutException("ERROR: This solver doesn't support constraints\n");
   }
 
@@ -483,7 +483,7 @@ void Solver::constraint(Vector2D& v, Vector2D& C_v, std::string name) {
   }
 #endif
 
-  if (!has_constraints) {
+  if (!supports_constraints) {
     throw BoutException("ERROR: This solver doesn't support constraints\n");
   }
 
@@ -526,7 +526,7 @@ void Solver::constraint(Vector3D& v, Vector3D& C_v, std::string name) {
   }
 #endif
 
-  if (!has_constraints) {
+  if (!supports_constraints) {
     throw BoutException("ERROR: This solver doesn't support constraints\n");
   }
 
diff --git a/tests/integrated/CMakeLists.txt b/tests/integrated/CMakeLists.txt
index 7a8ac9ab89..59a3091971 100644
--- a/tests/integrated/CMakeLists.txt
+++ b/tests/integrated/CMakeLists.txt
@@ -41,6 +41,7 @@ add_subdirectory(test-restarting)
 add_subdirectory(test-slepc-solver)
 add_subdirectory(test-smooth)
 add_subdirectory(test-snb)
+add_subdirectory(test-snes-constraints)
 add_subdirectory(test-snes-save-jacobian)
 add_subdirectory(test-solver)
 add_subdirectory(test-squash)
diff --git a/tests/integrated/test-snes-constraints/CMakeLists.txt b/tests/integrated/test-snes-constraints/CMakeLists.txt
new file mode 100644
index 0000000000..fcffcc6898
--- /dev/null
+++ b/tests/integrated/test-snes-constraints/CMakeLists.txt
@@ -0,0 +1,6 @@
+bout_add_integrated_test(
+  test-snes-constraints
+  SOURCES test_snes_constraints.cxx
+  REQUIRES BOUT_HAS_PETSC
+  USE_DATA_BOUT_INP
+)
diff --git a/tests/integrated/test-snes-constraints/data/BOUT.inp b/tests/integrated/test-snes-constraints/data/BOUT.inp
new file mode 100644
index 0000000000..1042c514f3
--- /dev/null
+++ b/tests/integrated/test-snes-constraints/data/BOUT.inp
@@ -0,0 +1,20 @@
+MXG = 0
+
+[mesh]
+nx = 1
+ny = 8
+nz = 1
+
+[solver]
+type = snes
+snes_type = newtonls
+equation_form = backward_euler
+timestep = 0.2
+output_step = 0.2
+nout = 1
+timestep_control = fixed
+atol = 1e-14
+rtol = 1e-14
+stol = 1e-14
+max_nonlinear_iterations = 20
+ksp_type = gmres
diff --git a/tests/integrated/test-snes-constraints/test_snes_constraints.cxx b/tests/integrated/test-snes-constraints/test_snes_constraints.cxx
new file mode 100644
index 0000000000..4a77dff3b9
--- /dev/null
+++ b/tests/integrated/test-snes-constraints/test_snes_constraints.cxx
@@ -0,0 +1,45 @@
+#include <bout/bout_types.hxx>
+#include <bout/field3d.hxx>
+#include <bout/physicsmodel.hxx>
+#include <bout/unused.hxx>
+
+/// Solves a Differential-Algebraic Equation (DAE) system
+///
+/// du/dt = -u + phi
+/// phi - constraint_factor * u = 0
+///
+/// Eliminating phi gives the continuous reduced equation
+///   du/dt = (-1 + constraint_factor) * u
+///
+/// with exact solution
+///   u = u0 * exp( (-1 + constraint_factor) * t)
+///
+/// This integrated test advances a single implicit timestep with the SNES
+/// solver, so the checked value is the one-step backward-Euler result rather
+/// than the continuous exact solution.
+class TestSnesConstraints : public PhysicsModel {
+public:
+  Field3D u;
+  Field3D phi;
+
+  static constexpr BoutReal constraint_factor = 0.5;
+
+  int init(bool UNUSED(restarting)) override {
+    solver->add(u, "u");
+    solver->constraint(phi, residual(phi), "phi");
+
+    u = 1.0;
+    phi = constraint_factor;
+
+    return 0;
+  }
+
+  int rhs(BoutReal UNUSED(time)) override {
+    ddt(u) = -u + phi;
+    residual(phi) = phi - constraint_factor * u;
+
+    return 0;
+  }
+};
+
+BOUTMAIN(TestSnesConstraints);
diff --git a/tests/integrated/test-snes-constraints/test_snes_constraints.py b/tests/integrated/test-snes-constraints/test_snes_constraints.py
new file mode 100644
index 0000000000..1260fb372e
--- /dev/null
+++ b/tests/integrated/test-snes-constraints/test_snes_constraints.py
@@ -0,0 +1,64 @@
+from pathlib import Path
+
+import numpy as np
+import pytest
+from boutdata.collect import collect
+from boututils.run_wrapper import launch
+
+TEST_DIR = Path(".")
+TIMESTEP = 0.2
+CONSTRAINT_FACTOR = 0.5
+# For one backward-Euler step:
+#   (u1 - u0) / dt = -u1 + phi1,  phi1 = c * u1
+# so u1 = u0 / (1 + (1 - c) * dt)
+EXPECTED_U = 1.0 / (1.0 + (1.0 - CONSTRAINT_FACTOR) * TIMESTEP)
+EXPECTED_PHI = CONSTRAINT_FACTOR * EXPECTED_U
+
+
+def assert_constraint_solution():
+    u = np.asarray(collect("u", path="data", info=False))[-1]
+    phi = np.asarray(collect("phi", path="data", info=False))[-1]
+
+    np.testing.assert_allclose(u, EXPECTED_U, atol=1e-12, rtol=0.0)
+    np.testing.assert_allclose(phi, EXPECTED_PHI, atol=1e-12, rtol=0.0)
+    np.testing.assert_allclose(phi - CONSTRAINT_FACTOR * u, 0.0, atol=1e-12, rtol=0.0)
+
+
+@pytest.mark.parametrize(
+    "equation_form",
+    ["backward_euler", "rearranged_backward_euler", "pseudo_transient"],
+)
+def test_constraint_equation_forms(equation_form):
+    status, output = launch(
+        f"./test_snes_constraints solver:equation_form={equation_form}",
+        nproc=1,
+        pipe=True,
+    )
+
+    if status:
+        print(f"WARNING: status = {status}")
+        print(output)
+
+    assert_constraint_solution()
+
+
+def test_constraint_fieldsplit():
+    status, output = launch(
+        "./test_snes_constraints "
+        "solver:equation_form=backward_euler "
+        "solver:pc_type=fieldsplit "
+        "petsc:pc_fieldsplit_type=additive "
+        "petsc:fieldsplit_diff_ksp_type=preonly "
+        "petsc:fieldsplit_diff_pc_type=jacobi "
+        "petsc:fieldsplit_alg_ksp_type=preonly "
+        "petsc:fieldsplit_alg_pc_type=jacobi",
+        nproc=1,
+        pipe=True,
+    )
+
+    if status:
+        print(f"WARNING: status = {status}")
+        print(output)
+
+    assert "Using PCFieldSplit preconditioner for DAE system" in output
+    assert_constraint_solution()
diff --git a/tests/unit/solver/test_fakesolver.hxx b/tests/unit/solver/test_fakesolver.hxx
index c8600e4649..3567774420 100644
--- a/tests/unit/solver/test_fakesolver.hxx
+++ b/tests/unit/solver/test_fakesolver.hxx
@@ -3,15 +3,20 @@
 
 #include "gtest/gtest.h"
 
+#include "bout/bout_types.hxx"
+#include "bout/boutexception.hxx"
+#include "bout/field2d.hxx"
+#include "bout/options.hxx"
 #include "bout/solver.hxx"
 
 #include <algorithm>
+#include <iterator>
 #include <string>
 #include <vector>
 
 class FakeSolver : public Solver {
 public:
-  FakeSolver(Options* options) : Solver(options) { has_constraints = true; }
+  FakeSolver(Options* options) : Solver(options) { supports_constraints = true; }
   ~FakeSolver() = default;
 
   int run() override {
@@ -37,7 +42,7 @@ public:
   }
   bool init_called{false};
 
-  void changeHasConstraints(bool new_value) { has_constraints = new_value; }
+  void changeHasConstraints(bool new_value) { supports_constraints = new_value; }
 
   auto listField2DNames() -> std::vector<std::string> {
     std::vector<std::string> result{};