diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 31877f8aac..66197b50ed 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/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((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(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(); @@ -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 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); @@ -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: { diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 40412f83b7..a749ef76a9 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/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 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?