Snes with constrains #3251

malamast · 2026-01-14T22:17:41Z

Work in progress.

Added some functionality to SNES solver to support constraints (solving an differential algebraic equation DAE). I based the implementation from other solvers:
Solvers which can include constraints in BOUT++ are:
IDA, from the SUNDIALS suite
IMEX-BDF2, which uses PETSc SNES
Currently, the constraints only work with equation_form = backward_euler. We will update the rest of the options later.
I added an option to split the pc and apply a different preconditioner for the algebraic equation.
I named the splited vectors as "diff" and "alg"
To enable the split, set in the imput file pc_type = fieldsplit

In the input file, one can do something like this and explore different petsc options:

  [solver]
  nout = 5          # Number of output steps
  output_step = 100.0    # Output timestep, normalised ion cyclotron times [1/Omega_ci]
  
  type = snes                           
  snes_type = newtonls            
  equation_form = backward_euler  
  ksp_type = gmres             
  max_nonlinear_iterations = 26  
  pc_type = fieldsplit            
  lag_jacobian = 12               
  atol = 1e-12                   
  rtol = 1e-5                
  stol = 1e-18
  maxl = 200                   
  diagnose = true
  pid_controller = true
  target_its = 8
  matrix_free_operator = false
  
  [petsc]
  # Preconditioner spliting
  pc_fieldsplit_type = multiplicative # additive,multiplicative,symmetric_multiplicative,schur,gkb
  
  fieldsplit_diff_ksp_type = preonly
  fieldsplit_diff_pc_type = hypre
  fieldsplit_diff_pc_hypre_type = ilu
  fieldsplit_diff_pc_hypre_ilu_type = Block-Jacobi-ILUT
  fieldsplit_diff_pc_hypre_ilu_drop_threshold = 1e-3
  fieldsplit_diff_pc_hypre_ilu_max_nnz_per_row = 70
  fieldsplit_diff_pc_hypre_ilu_local_reordering = true
  fieldsplit_diff_pc_hypre_ilu_tri_solve = true
  
  fieldsplit_alg_ksp_type = preonly
  fieldsplit_alg_pc_type  = hypre
  fieldsplit_alg_pc_hypre_type = boomeramg
  
  fieldsplit_alg_pc_hypre_boomeramg_cycle_type = V
  fieldsplit_alg_pc_hypre_boomeramg_coarsen_type = HMIS
  fieldsplit_alg_pc_hypre_boomeramg_interp_type  = ext+i
  fieldsplit_alg_pc_hypre_boomeramg_relax_type_all =  l1scaled-Jacobi
  fieldsplit_alg_pc_hypre_boomeramg_strong_threshold = 0.5
  
  fieldsplit_alg_pc_hypre_boomeramg_max_iter = 1
  fieldsplit_alg_pc_hypre_boomeramg_tol      = 0.0
  #fieldsplit_alg_pc_hypre_boomeramg_print_statistics = 2

…on of algebraic equations. For now it is only supported by backward_euler option. We will modify the rest of the options later.

…ork 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.

github-actions

clang-tidy made some suggestions

github-actions · 2026-01-14T22:20:04Z

src/solver/impls/snes/snes.cxx

+  if (have_constraints) {
+    // CreatePETSc-native index sets representing the two parts of your DAE.
+    PetscInt istart, iend;
+    PetscCall(VecGetOwnershipRange(snes_x, &istart, &iend));


warning: multiple declarations in a single statement reduces readability [readability-isolate-declaration]

Suggested change

PetscCall(VecGetOwnershipRange(snes_x, &istart, &iend));

PetscInt istart;

PetscInt iend;

github-actions · 2026-01-14T22:20:04Z

src/solver/impls/snes/snes.cxx

+  if (have_constraints) {
+    // CreatePETSc-native index sets representing the two parts of your DAE.
+    PetscInt istart, iend;
+    PetscCall(VecGetOwnershipRange(snes_x, &istart, &iend));


warning: variable 'iend' is not initialized [cppcoreguidelines-init-variables]

PetscInt istart, iend; ^

this fix will not be applied because it overlaps with another fix

github-actions · 2026-01-14T22:20:04Z

src/solver/impls/snes/snes.cxx

+  if (have_constraints) {
+    // CreatePETSc-native index sets representing the two parts of your DAE.
+    PetscInt istart, iend;
+    PetscCall(VecGetOwnershipRange(snes_x, &istart, &iend));


warning: variable 'istart' is not initialized [cppcoreguidelines-init-variables]

PetscInt istart, iend; ^

this fix will not be applied because it overlaps with another fix

github-actions · 2026-01-14T22:20:04Z

src/solver/impls/snes/snes.cxx

+    }
+
+    PetscCall(ISCreateGeneral(BoutComm::get(), diff_idx.size(), diff_idx.data(),
+                              PETSC_COPY_VALUES, &is_diff));


warning: no header providing "ISCreateGeneral" is directly included [misc-include-cleaner]

src/solver/impls/snes/snes.cxx:1:

+ #include <petscis.h>

github-actions · 2026-01-14T22:20:04Z

src/solver/impls/snes/snes.cxx

+    // Use PETSc fieldsplit
+    PetscCall(PCSetType(pc, PCFIELDSPLIT));
+    exit(0);
+


warning: no header providing "PCFIELDSPLIT" is directly included [misc-include-cleaner]

PetscCall(PCSetType(pc, PCFIELDSPLIT)); ^

github-actions · 2026-01-14T22:20:05Z

src/solver/impls/snes/snes.cxx

+
+      Vec x_diff, x0_diff, f_diff;
+
+      PetscCall(VecGetSubVector(x, is_diff, &x_diff));


warning: variable 'x0_diff' is not initialized [cppcoreguidelines-init-variables]

ints ^

this fix will not be applied because it overlaps with another fix

github-actions · 2026-01-14T22:20:05Z

src/solver/impls/snes/snes.cxx

+
+      Vec x_diff, x0_diff, f_diff;
+
+      PetscCall(VecGetSubVector(x, is_diff, &x_diff));


warning: variable 'x_diff' is not initialized [cppcoreguidelines-init-variables]

ints ^

this fix will not be applied because it overlaps with another fix

github-actions · 2026-01-14T22:20:05Z

src/solver/impls/snes/snes.cxx

+
+      PetscCall(VecGetSubVector(x, is_diff, &x_diff));
+      PetscCall(VecGetSubVector(x0, is_diff, &x0_diff));
+      PetscCall(VecGetSubVector(f, is_diff, &f_diff));


warning: multiple declarations in a single statement reduces readability [readability-isolate-declaration]

Suggested change

PetscCall(VecGetSubVector(f, is_diff, &f_diff));

ints

diff;

iVec x_diff;

Vec x0_diff;

Vec f_diff;

github-actions · 2026-01-14T22:20:05Z

src/solver/impls/snes/snes.hxx

  int neq;    ///< Number of variables in total

+  bool have_constraints;  ///< Are there any constraint variables?
+  Array<BoutReal> is_dae; ///< If using constraints, 1 -> DAE, 0 -> AE


warning: no header providing "Array" is directly included [misc-include-cleaner]

src/solver/impls/snes/snes.hxx:30:

- #include <bout/build_defines.hxx> + #include "bout/array.hxx" + #include <bout/build_defines.hxx>

github-actions · 2026-01-14T22:20:05Z

src/solver/impls/snes/snes.hxx

+  bool have_constraints;  ///< Are there any constraint variables?
+  Array<BoutReal> is_dae; ///< If using constraints, 1 -> DAE, 0 -> AE
+
+  IS is_diff = nullptr; // is_dae == 1


warning: no header providing "IS" is directly included [misc-include-cleaner]

src/solver/impls/snes/snes.hxx:32:

+ #include <petscistypes.h>

…ved.

github-actions

clang-tidy made some suggestions

github-actions · 2026-01-14T22:33:53Z

src/solver/impls/snes/snes.cxx

+    // Use PETSc fieldsplit
+    PetscCall(PCSetType(pc, PCFIELDSPLIT));
+
+    // Give PETSc the index sets


warning: no header providing "PCFIELDSPLIT" is directly included [misc-include-cleaner]

PetscCall(PCSetType(pc, PCFIELDSPLIT)); ^

github-actions · 2026-01-14T22:33:53Z

src/solver/impls/snes/snes.cxx

+
+      Vec x_diff, x0_diff, f_diff;
+
+      PetscCall(VecGetSubVector(x, is_diff, &x_diff));


warning: variable 'f_diff' is not initialized [cppcoreguidelines-init-variables]

ints ^

this fix will not be applied because it overlaps with another fix

github-actions · 2026-01-14T22:33:54Z

src/solver/impls/snes/snes.cxx

+
+      Vec x_diff, x0_diff, f_diff;
+
+      PetscCall(VecGetSubVector(x, is_diff, &x_diff));


warning: variable 'x0_diff' is not initialized [cppcoreguidelines-init-variables]

ints ^

this fix will not be applied because it overlaps with another fix

github-actions · 2026-01-14T22:33:54Z

src/solver/impls/snes/snes.cxx

+
+      Vec x_diff, x0_diff, f_diff;
+
+      PetscCall(VecGetSubVector(x, is_diff, &x_diff));


warning: variable 'x_diff' is not initialized [cppcoreguidelines-init-variables]

ints ^

this fix will not be applied because it overlaps with another fix

github-actions · 2026-01-14T22:33:54Z

src/solver/impls/snes/snes.cxx

+
+      PetscCall(VecGetSubVector(x, is_diff, &x_diff));
+      PetscCall(VecGetSubVector(x0, is_diff, &x0_diff));
+      PetscCall(VecGetSubVector(f, is_diff, &f_diff));


warning: multiple declarations in a single statement reduces readability [readability-isolate-declaration]

Suggested change

PetscCall(VecGetSubVector(f, is_diff, &f_diff));

ints

diff;

iVec x_diff;

Vec x0_diff;

Vec f_diff;

mikekryjak · 2026-01-15T13:25:44Z

Hi @malamast, this looks exciting! I'd like to make sure I understand it. What are these constraints and what do they do? What do you mean by preconditioner splitting?

malamast · 2026-01-15T17:46:34Z

Hi @malamast, this looks exciting! I'd like to make sure I understand it. What are these constraints and what do they do? What do you mean by preconditioner splitting?

Hi @mikekryjak ,

by constraints here I mean algebraic equations in a differential–algebraic equation (DAE) system. In our case, this is used to include the electrostatic potential φ directly in the SNES state vector and to solve its Poisson-like equation (linking φ to vorticity) implicitly together with the rest of the time-dependent plasma equations.

This is to be used along with this PR: boutproject/hermes-3#470

Previously, φ was computed inside the RHS by explicitly inverting a Laplacian at every residual evaluation (to recover φ from vorticity). That approach is expensive for 2D and also makes the nonlinear solve harder, since the Poisson solve is hidden inside the residual.

With this change, φ is treated as an unknown with an algebraic constraint (its Poisson-like equation), while the other fields remain governed by differential equations. SNES then solves the coupled DAE system in one nonlinear solve, avoiding repeated Laplacian inversions in the RHS and reducing the overall computational cost.

By preconditioner splitting, I mean separating the Jacobian into blocks corresponding to:

the differential equations (transport, momentum, energy, etc.), and
the algebraic constraint for φ.

This allows applying a different (and more appropriate) preconditioner to the φ block (e.g. one tuned for elliptic/Poisson-like operators) while keeping the existing preconditioner for the rest of the system. I've tried to use Hypre/BoomerAMG for phi but with no success so far.

malamast added 2 commits January 13, 2026 16:15

snes: changes to the solver that it can handle constraints for soluti…

5f5b1a0

…on of algebraic equations. For now it is only supported by backward_euler option. We will modify the rest of the options later.

snes: I added petsc index sets IS from is_dae vector so that we can w…

b26bd50

…ork 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.

malamast requested review from ZedThree and bendudson January 14, 2026 22:17

malamast self-assigned this Jan 14, 2026

malamast added enhance request work in progress Not ready for merging labels Jan 14, 2026

github-actions bot reviewed Jan 14, 2026

View reviewed changes

snes: I had accidentally included an exit statement which is now remo…

74a5480

…ved.

malamast force-pushed the snes-with-constrains branch from 221deee to 74a5480 Compare January 14, 2026 22:31

Apply clang-format changes

d0f64f2

github-actions bot reviewed Jan 14, 2026

View reviewed changes

	PetscCall(VecGetOwnershipRange(snes_x, &istart, &iend));
	PetscInt istart;
	PetscInt iend;


		Vec x_diff, x0_diff, f_diff;

		PetscCall(VecGetSubVector(x, is_diff, &x_diff));

-      PetscCall(VecGetSubVector(f, is_diff, &f_diff));
+ints
+diff;
+iVec x_diff;
+      Vec x0_diff;
+      Vec f_diff;

Snes with constrains #3251

Are you sure you want to change the base?

Snes with constrains #3251

Uh oh!

Conversation

malamast commented Jan 14, 2026

Uh oh!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Jan 14, 2026

Choose a reason for hiding this comment

Uh oh!

mikekryjak commented Jan 15, 2026

Uh oh!

malamast commented Jan 15, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants