From fdb23f76c921d59838da324bfb905fa2e31fc247 Mon Sep 17 00:00:00 2001 From: malamast Date: Tue, 13 Jan 2026 16:15:03 -0800 Subject: [PATCH 1/6] snes: changes to the solver that it can handle constraints for solution of algebraic equations. For now it is only supported by backward_euler option. We will modify the rest of the options later. --- src/solver/impls/snes/snes.cxx | 70 ++++++++++++++++++++++++++++++++-- src/solver/impls/snes/snes.hxx | 3 ++ 2 files changed, 69 insertions(+), 4 deletions(-) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 9a53aabfcb..af65e9d992 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -356,7 +356,9 @@ SNESSolver::SNESSolver(Options* opts) .withDefault(false)), asinh_vars((*options)["asinh_vars"] .doc("Apply asinh() to all variables?") - .withDefault(false)) {} + .withDefault(false)) { + has_constraints = true; ///< This solver can handle constraints +} int SNESSolver::init() { Solver::init(); @@ -377,6 +379,29 @@ 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 + have_constraints = false; + + for (int i = 0; i < n2Dvars(); i++) { + if (f2d[i].constraint) { + have_constraints = true; + break; + } + } + for (int i = 0; i < n3Dvars(); i++) { + if (f3d[i].constraint) { + have_constraints = true; + break; + } + } + + if (have_constraints) { + 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 PETSc components // Vectors @@ -503,6 +528,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); @@ -1403,7 +1431,13 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) { // Call the RHS function if (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; } @@ -1429,8 +1463,36 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) { 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 + if (!have_constraints) { + + VecAYPX(f, -dt, x); // f <- x - Δt*f + VecAXPY(f, -1.0, x0); // f <- f - x0 + + } else { + + const BoutReal* xdata = nullptr; + PetscCall(VecGetArrayRead(x, &xdata)); + + const BoutReal* x0data = nullptr; + PetscCall(VecGetArrayRead(x0, &x0data)); + + // Copy derivatives + BoutReal* fdata = nullptr; + PetscCall(VecGetArray(f, &fdata)); + + // Some constraints + for (int i = 0; i < nlocal; i++) { + if (is_dae[i] > 0.5) { // 1 -> differential, 0 -> algebraic + fdata[i] = xdata[i] - x0data[i] - dt * fdata[i]; + } + // Otherwise the constraint is that fdata[i] = rhs (residual of algebraic equation) + } + + // PetscCall(VecRestoreArrayRead(x, &xdata)); // NOTE: Do we need this? + // PetscCall(VecRestoreArrayRead(x0, &x0data)); // NOTE: Do we need this? + PetscCall(VecRestoreArray(f, &fdata)); + } + break; } case BoutSnesEquationForm::direct_newton: { diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index ddb447e865..20de6bd75f 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -217,6 +217,9 @@ private: int nlocal; ///< Number of variables on local processor int neq; ///< Number of variables in total + bool have_constraints; ///< Are there any constraint variables? + Array is_dae; ///< If using constraints, 1 -> DAE, 0 -> AE + PetscLib lib; ///< Handles initialising, finalising PETSc Vec snes_f; ///< Used by SNES to store function Vec snes_x; ///< Result of SNES From bdc5ce15ce7b23651c15839249d16b57b95629d1 Mon Sep 17 00:00:00 2001 From: malamast Date: Wed, 14 Jan 2026 14:00:32 -0800 Subject: [PATCH 2/6] snes: I added petsc index sets IS from is_dae vector so that we can work with petsc vectors instead of bout++ variables. I added an option to split the preconditioner so that the user can use a different PC_type for the algebraic equation. --- src/solver/impls/snes/snes.cxx | 82 +++++++++++++++++++++++++--------- src/solver/impls/snes/snes.hxx | 4 ++ 2 files changed, 66 insertions(+), 20 deletions(-) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index af65e9d992..9aa7447db0 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -452,6 +452,34 @@ int SNESSolver::init() { local_residual_2d = 0.0; global_residual = 0.0; + if (have_constraints) { + // 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 diff_idx; + std::vector 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)); + + have_is_maps = true; + } + // Nonlinear solver interface (SNES) output_info.write("Create SNES\n"); SNESCreate(BoutComm::get(), &snes); @@ -583,6 +611,25 @@ int SNESSolver::init() { } } + if (have_constraints && have_is_maps && !matrix_free && pc_type == "fieldsplit") { + output_info.write("Using PCFieldSplit preconditioner for DAE system\n"); + + // Use PETSc fieldsplit + PetscCall(PCSetType(pc, PCFIELDSPLIT)); + exit(0); + + // Give PETSc the index sets + PetscCall(PCFieldSplitSetIS(pc, "diff", is_diff)); + PetscCall(PCFieldSplitSetIS(pc, "alg", is_alg)); + + // Let the user configure from options (recommended) + // Example options you can set in input file: + // -pc_fieldsplit_type additive + // -fieldsplit_alg_pc_type hypre -fieldsplit_alg_pc_hypre_type boomeramg + // -fieldsplit_diff_pc_type ilu + // + } + // Get runtime options lib.setOptionsFromInputFile(snes); @@ -1461,8 +1508,10 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) { break; } case BoutSnesEquationForm::backward_euler: { - // Backward Euler - // Set f = x - x0 - Δt*f + // Backward Euler: + // Differential: F = x - x0 - dt*f + // Algebraic: F = G(x) (already stored in f by rhs_function) + if (!have_constraints) { VecAYPX(f, -dt, x); // f <- x - Δt*f @@ -1470,29 +1519,22 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) { } else { - const BoutReal* xdata = nullptr; - PetscCall(VecGetArrayRead(x, &xdata)); + ASSERT2(have_is_maps); + // Some constraints - const BoutReal* x0data = nullptr; - PetscCall(VecGetArrayRead(x0, &x0data)); + Vec x_diff, x0_diff, f_diff; - // Copy derivatives - BoutReal* fdata = nullptr; - PetscCall(VecGetArray(f, &fdata)); + PetscCall(VecGetSubVector(x, is_diff, &x_diff)); + PetscCall(VecGetSubVector(x0, is_diff, &x0_diff)); + PetscCall(VecGetSubVector(f, is_diff, &f_diff)); - // Some constraints - for (int i = 0; i < nlocal; i++) { - if (is_dae[i] > 0.5) { // 1 -> differential, 0 -> algebraic - fdata[i] = xdata[i] - x0data[i] - dt * fdata[i]; - } - // Otherwise the constraint is that fdata[i] = rhs (residual of algebraic equation) - } + PetscCall(VecAYPX(f_diff, -dt, x_diff)); // f_diff <- x_diff - dt*f_diff + PetscCall(VecAXPY(f_diff, -1.0, x0_diff)); // f_diff <- f_diff - x0_diff - // PetscCall(VecRestoreArrayRead(x, &xdata)); // NOTE: Do we need this? - // PetscCall(VecRestoreArrayRead(x0, &x0data)); // NOTE: Do we need this? - PetscCall(VecRestoreArray(f, &fdata)); + PetscCall(VecRestoreSubVector(x, is_diff, &x_diff)); + PetscCall(VecRestoreSubVector(x0, is_diff, &x0_diff)); + PetscCall(VecRestoreSubVector(f, is_diff, &f_diff)); } - break; } case BoutSnesEquationForm::direct_newton: { diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 20de6bd75f..53d39c2123 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -220,6 +220,10 @@ private: bool have_constraints; ///< Are there any constraint variables? Array 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) + bool have_is_maps = false; + PetscLib lib; ///< Handles initialising, finalising PETSc Vec snes_f; ///< Used by SNES to store function Vec snes_x; ///< Result of SNES From afe122196c19041ee0996a2ac4d5ec56386c8e64 Mon Sep 17 00:00:00 2001 From: malamast Date: Wed, 14 Jan 2026 14:29:33 -0800 Subject: [PATCH 3/6] snes: I had accidentally included an exit statement which is now removed. --- src/solver/impls/snes/snes.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 9aa7447db0..62cf3e5cf5 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -616,7 +616,6 @@ int SNESSolver::init() { // Use PETSc fieldsplit PetscCall(PCSetType(pc, PCFIELDSPLIT)); - exit(0); // Give PETSc the index sets PetscCall(PCFieldSplitSetIS(pc, "diff", is_diff)); From c4eb14e40b4cc19d35189be18ed4c0e53e091e0d Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 10 Jun 2026 23:11:23 -0700 Subject: [PATCH 4/6] snes: require backward_euler for constraints --- src/solver/impls/snes/snes.cxx | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 62cf3e5cf5..3d7cb2ad28 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -21,6 +21,7 @@ #include #include #include +#include #include #include "petscerror.h" @@ -400,6 +401,11 @@ int SNESSolver::init() { // Call the Solver function, which sets the array // to one when not a constraint, zero for constraint set_id(std::begin(is_dae)); + + if (equation_form != BoutSnesEquationForm::backward_euler) { + throw BoutException( + "SNES constraints currently require equation_form=backward_euler"); + } } // Initialise PETSc components From 1d1b0aa3aed7f1696178fb4ab81fbf51dcabc97b Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 11 Jun 2026 10:48:00 -0700 Subject: [PATCH 5/6] Small tidying --- src/solver/impls/cvode/cvode.cxx | 1 - src/solver/impls/cvode/cvode.hxx | 2 -- src/solver/impls/snes/snes.cxx | 2 +- 3 files changed, 1 insertion(+), 4 deletions(-) diff --git a/src/solver/impls/cvode/cvode.cxx b/src/solver/impls/cvode/cvode.cxx index 5cc5cae59a..4677ab3c14 100644 --- a/src/solver/impls/cvode/cvode.cxx +++ b/src/solver/impls/cvode/cvode.cxx @@ -1,7 +1,6 @@ /************************************************************************** * Interface to SUNDIALS CVODE * - * ************************************************************************** * Copyright 2010-2026 BOUT++ contributors * diff --git a/src/solver/impls/cvode/cvode.hxx b/src/solver/impls/cvode/cvode.hxx index d1409ad9a3..1160df1cf0 100644 --- a/src/solver/impls/cvode/cvode.hxx +++ b/src/solver/impls/cvode/cvode.hxx @@ -1,8 +1,6 @@ /************************************************************************** * Interface to SUNDIALS CVODE * - * NOTE: Only one solver can currently be compiled in - * ************************************************************************** * Copyright 2010 - 2026 BOUT++ contributors * diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 3d7cb2ad28..90a94d8583 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -459,7 +459,7 @@ int SNESSolver::init() { global_residual = 0.0; if (have_constraints) { - // CreatePETSc-native index sets representing the two parts of your DAE. + // Create PETSc-native index sets representing the two parts of your DAE. PetscInt istart, iend; PetscCall(VecGetOwnershipRange(snes_x, &istart, &iend)); ASSERT2(iend - istart == nlocal); From e202f8e817bb75457d0de5d87c384b2d4ca97b15 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 11 Jun 2026 11:01:47 -0700 Subject: [PATCH 6/6] snes: Update manual to include DAE/constraints --- manual/sphinx/user_docs/time_integration.rst | 63 +++++++++++++++++++- 1 file changed, 61 insertions(+), 2 deletions(-) diff --git a/manual/sphinx/user_docs/time_integration.rst b/manual/sphinx/user_docs/time_integration.rst index 8dd16f9c7b..7f54e40f65 100644 --- a/manual/sphinx/user_docs/time_integration.rst +++ b/manual/sphinx/user_docs/time_integration.rst @@ -419,8 +419,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 useful + for robust convergence to steady state. .. code-block:: ini @@ -429,6 +430,14 @@ The solver supports several timestepping strategies controlled by ``equation_for This method has low accuracy in time but its dissipative properties are helpful when evolving to steady state solutions. +**Backward Euler (DAE/constraints)** + Backward Euler with explicit support for algebraic constraint variables. This form is + required when using BOUT++ ``Solver::constraint(...)`` with the SNES solver. + + .. code-block:: ini + + equation_form = backward_euler + **Direct Newton** Solves the steady-state problem F(u) = 0 directly without timestepping. @@ -450,6 +459,56 @@ 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: advanced with backward Euler +- Algebraic (constraint) variables: solved from the algebraic equations ``G(x) = 0`` at each step + +Current limitations: + +- Constraints are supported only with ``equation_form = backward_euler``. +- Constraint splitting requires ``matrix_free = false`` (``matrix_free_operator`` may still be + used). + +Preconditioner splitting +^^^^^^^^^^^^^^^^^^^^^^^^ + +When constraints are enabled, the SNES solver can optionally split the preconditioner into a +differential part and an algebraic part using PETSc ``fieldsplit``. The split names are: + +- ``diff``: differential variables +- ``alg``: algebraic (constraint) variables + +To enable the split, set ``pc_type = fieldsplit`` in the ``[solver]`` section, then configure PETSc +using standard fieldsplit options (in ``[petsc]`` or on the command line). + +Example: + +.. code-block:: ini + + [solver] + type = snes + equation_form = backward_euler + pc_type = fieldsplit + + [petsc] + # Preconditioner splitting + pc_fieldsplit_type = multiplicative # additive, multiplicative, schur, gkb, ... + + # Differential block + fieldsplit_diff_ksp_type = preonly + fieldsplit_diff_pc_type = hypre + fieldsplit_diff_pc_hypre_type = ilu + + # Algebraic (constraint) block + fieldsplit_alg_ksp_type = preonly + fieldsplit_alg_pc_type = hypre + fieldsplit_alg_pc_hypre_type = boomeramg + Adaptive Timestepping ~~~~~~~~~~~~~~~~~~~~~