Snes with constrains

src/solver/impls/snes/snes.cxx

@@ -28,10 +28,10 @@
 
 class ColoringStencil {
 private:
-  bool static isInSquare(int const i, int const j, int const n_square) {
+  bool static isInSquare(const int i, const int j, const int n_square) {
     return std::abs(i) <= n_square && std::abs(j) <= n_square;
   }
-  bool static isInCross(int const i, int const j, int const n_cross) {
+  bool static isInCross(const int i, const int j, const int n_cross) {
     if (i == 0) {
       return std::abs(j) <= n_cross;
     }
@@ -40,7 +40,7 @@ class ColoringStencil {
     }
     return false;
   }
-  bool static isInTaxi(int const i, int const j, int const n_taxi) {
+  bool static isInTaxi(const int i, const int j, const int n_taxi) {
     return std::abs(i) + std::abs(j) <= n_taxi;
   }
 
@@ -154,7 +154,7 @@ PetscErrorCode SNESSolver::FDJinitialise() {
                       .doc("Extent of stencil (taxi-cab norm)")
                       .withDefault<int>((n_square == 0 && n_cross == 0) ? 2 : 0);
 
-    auto const xy_offsets = ColoringStencil::getOffsets(n_square, n_taxi, n_cross);
+    const auto xy_offsets = ColoringStencil::getOffsets(n_square, n_taxi, n_cross);
     {
       // This is ugly but can't think of a better and robust way to
       // count the non-zeros for some arbitrary stencil
@@ -624,7 +624,9 @@ SNESSolver::SNESSolver(Options* opts)
                      .withDefault<bool>(false)),
       asinh_vars((*options)["asinh_vars"]
                      .doc("Apply asinh() to all variables?")
-                     .withDefault<bool>(false)) {}
+                     .withDefault<bool>(false)) {
+  has_constraints = true; ///< This solver can handle constraints
+}
 
 int SNESSolver::init() {
   Solver::init();
@@ -645,6 +647,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
@@ -695,6 +720,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<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));
+
+    have_is_maps = true;
+  }
+
   // Nonlinear solver interface (SNES)
   output_info.write("Create SNES\n");
   SNESCreate(BoutComm::get(), &snes);
@@ -771,6 +824,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);
@@ -823,6 +879,24 @@ 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));
+
+    // 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);
 
@@ -1667,7 +1741,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;
   }
@@ -1691,10 +1771,33 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) {
     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)
+
+    if (!have_constraints) {
+
+      VecAYPX(f, -dt, x);   // f <- x - Δt*f
+      VecAXPY(f, -1.0, x0); // f <- f - x0
+
+    } else {
+
+      ASSERT2(have_is_maps);
+      // Some constraints
+
+      Vec x_diff, x0_diff, f_diff;
+
+      PetscCall(VecGetSubVector(x, is_diff, &x_diff));
+      PetscCall(VecGetSubVector(x0, is_diff, &x0_diff));
+      PetscCall(VecGetSubVector(f, is_diff, &f_diff));
+
+      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(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: {

src/solver/impls/snes/snes.hxx

@@ -201,8 +201,8 @@ private:
   BoutReal kI; ///< (0.2 - 0.4) Integral parameter (smooths history of changes)
   BoutReal kD; ///< (0.1 - 0.3) Derivative (dampens oscillation - optional)
   bool pid_consider_failures; ///< Reduce timestep increases if recent solves have failed
-  BoutReal recent_failure_rate;            ///< Rolling average of recent failure rate
-  BoutReal last_failure_weight;            ///< 1 / number of recent solves
+  BoutReal recent_failure_rate; ///< Rolling average of recent failure rate
+  BoutReal last_failure_weight; ///< 1 / number of recent solves
 
   BoutReal nl_its_prev;
   BoutReal nl_its_prev2;
@@ -215,6 +215,13 @@ 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<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)
+  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
@@ -244,7 +251,7 @@ private:
   bool matrix_free_operator; ///< Use matrix free Jacobian in the operator?
   int lag_jacobian;          ///< Re-use Jacobian
   bool jacobian_persists; ///< Re-use Jacobian and preconditioner across nonlinear solves
-  bool use_coloring;         ///< Use matrix coloring
+  bool use_coloring;      ///< Use matrix coloring
 
   bool jacobian_recalculated; ///< Flag set when Jacobian is recalculated
   bool prune_jacobian;        ///< Remove small elements in the Jacobian?