Snes with constrains updated#3421
Conversation
Supports solution of algebraic differential equations (e.g. for phi) also for the other equation froms like rearranged_backward_euler and the pseudo_transient method as well.
…nl_its and that of the previous steps is equal to the target number of non-linear iterations.
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
|
|
||
| if (have_constraints) { | ||
| // CreatePETSc-native index sets representing the two parts of your DAE. | ||
| PetscInt istart, iend; |
There was a problem hiding this comment.
warning: multiple declarations in a single statement reduces readability [readability-isolate-declaration]
| PetscInt istart, iend; | |
| PetscInt istart; | |
| PetscInt iend; |
|
|
||
| if (have_constraints) { | ||
| // CreatePETSc-native index sets representing the two parts of your DAE. | ||
| PetscInt istart, iend; |
There was a problem hiding this comment.
warning: variable 'iend' is not initialized [cppcoreguidelines-init-variables]
PetscInt istart, iend;
^this fix will not be applied because it overlaps with another fix
|
|
||
| if (have_constraints) { | ||
| // CreatePETSc-native index sets representing the two parts of your DAE. | ||
| PetscInt istart, iend; |
There was a problem hiding this comment.
warning: variable 'istart' is not initialized [cppcoreguidelines-init-variables]
PetscInt istart, iend;
^this fix will not be applied because it overlaps with another fix
| ASSERT2(have_is_maps); | ||
| // Some constraints | ||
|
|
||
| Vec x_diff, x0_diff, delta_x_diff, f_diff; |
There was a problem hiding this comment.
warning: variable 'delta_x_diff' is not initialized [cppcoreguidelines-init-variables]
s
^this fix will not be applied because it overlaps with another fix
| ASSERT2(have_is_maps); | ||
| // Some constraints | ||
|
|
||
| Vec x_diff, x0_diff, delta_x_diff, f_diff; |
There was a problem hiding this comment.
warning: variable 'f_diff' is not initialized [cppcoreguidelines-init-variables]
s
^this fix will not be applied because it overlaps with another fix
| ASSERT2(have_is_maps); | ||
| // Some constraints | ||
|
|
||
| Vec x_diff, x0_diff, f_diff; |
There was a problem hiding this comment.
warning: variable 'f_diff' is not initialized [cppcoreguidelines-init-variables]
ints
^this fix will not be applied because it overlaps with another fix
| ASSERT2(have_is_maps); | ||
| // Some constraints | ||
|
|
||
| Vec x_diff, x0_diff, f_diff; |
There was a problem hiding this comment.
warning: variable 'x0_diff' is not initialized [cppcoreguidelines-init-variables]
ints
^this fix will not be applied because it overlaps with another fix
| ASSERT2(have_is_maps); | ||
| // Some constraints | ||
|
|
||
| Vec x_diff, x0_diff, f_diff; |
There was a problem hiding this comment.
warning: variable 'x_diff' is not initialized [cppcoreguidelines-init-variables]
ints
^this fix will not be applied because it overlaps with another fix
| // Some constraints | ||
|
|
||
| Vec x_diff, x0_diff, f_diff; | ||
| PetscCall(VecGetSubVector(x, is_diff, &x_diff)); |
There was a problem hiding this comment.
warning: multiple declarations in a single statement reduces readability [readability-isolate-declaration]
| PetscCall(VecGetSubVector(x, is_diff, &x_diff)); | |
| ints | |
| diff;Vec x_diff; | |
| Vec x0_diff; | |
| Vec f_diff; |
| 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 |
There was a problem hiding this comment.
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>
cmacmackin
left a comment
cmacmackin
left a comment
There was a problem hiding this comment.
There are some changes required and some things I think need further discussion. The biggest problem, though, is that there is not test for any of this. You should add an integration test.
| // Add slow growth when convergence is good and stable | ||
| if (nl_its <= target_its && nl_its_prev <= target_its) { | ||
| fac *= 1.1; // or 1.05 for more conservative growth | ||
| } | ||
|
|
There was a problem hiding this comment.
I feel like there should be more discussion before adopting this. We are aware that timestep setting tends to be too conservative. So far we've been dealing with it by setting target_its to be much larger than it should be, as seen here. That certainly isn't a good long-term solution, but I think there were some other ideas for how this could be handled. @mikekryjak What are your thoughts?
| n3Dvars(), n2Dvars(), neq, nlocal); | ||
|
|
||
| // Check if there are any constraints | ||
| have_constraints = false; |
There was a problem hiding this comment.
Just give this a default value of false, unless you expect that the init() method might be called more than once.
| have_constraints = false; |
| int nlocal; ///< Number of variables on local processor | ||
| int neq; ///< Number of variables in total | ||
|
|
||
| bool have_constraints; ///< Are there any constraint variables? |
There was a problem hiding this comment.
| bool have_constraints; ///< Are there any constraint variables? | |
| bool have_constraints = false; ///< Are there any constraint variables? |
| n3Dvars(), n2Dvars(), neq, nlocal); | ||
|
|
||
| // Check if there are any constraints | ||
| have_constraints = false; |
There was a problem hiding this comment.
I find have_constraints is too similar a name to has_constraints (I was confused by this during my review). Could you choose another name, such as algebraic_constraints?
| PetscCall(ISCreateGeneral(BoutComm::get(), alg_idx.size(), alg_idx.data(), | ||
| PETSC_COPY_VALUES, &is_alg)); | ||
|
|
||
| have_is_maps = true; |
There was a problem hiding this comment.
What is the purpose of this variable? As far as I can tell, if have_constraints is true then this will always also be true. Is there any situation when the value of have_is_maps and have_constraints could differ?
| // 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 |
There was a problem hiding this comment.
It would be better to put these notes in documentation, rather than as a comment in the code.
| // F = (x0 - x)/Δt + f | ||
| // Algebraic: F = G(x) (already stored in f by rhs_function) | ||
|
|
||
| if (!have_constraints) { |
There was a problem hiding this comment.
There's lots of repetition in this part of the code. I've drafted an alternative approach which should simplify things, but I don't have any problems to test it on. I'll post the diff in a separate comment.
|
Here is a proposed refactor to avoid repeating yourself so much in diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index bdb32ab04..657000ca1 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -1567,9 +1567,27 @@ PetscErrorCode SNESSolver::rhs_function(Vec x, Vec f, bool linear) {
return PETSC_SUCCESS;
}
+template <typename Func, typename... Args>
+PetscErrorCode performCalculation(Func operation, IS indices, Args... args) {
+ if (indices == nullptr) {
+ PetscCall(operation(args...));
+ } else {
+ constexpr std::size_t N = sizeof...(Args);
+ std::array<Vec, N> vectors{args...};
+ std::array<Vec, N> subvectors;
+ for (std::size_t i = 0; i < N; ++i) {
+ PetscCall(VecGetSubVector(vectors[i], indices, &(subvectors[i])));
+ }
+ PetscCall(std::apply(operation, subvectors));
+ for (std::size_t i = 0; i < N; ++i) {
+ PetscCall(VecRestoreSubVector(vectors[i], indices, &(subvectors[i])));
+ }
+ }
+ return PETSC_SUCCESS;
+}
+
// Result in f depends on equation_form
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
@@ -1583,34 +1601,13 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) {
// Rearranged Backward Euler
// F = (x0 - x)/Δt + f
// Algebraic: F = G(x) (already stored in f by rhs_function)
-
- if (!have_constraints) {
-
- // First calculate x - x0 to minimise floating point issues
- VecWAXPY(delta_x, -1.0, x0, x); // delta_x = x - x0
- VecAXPY(f, -1.0 / dt, delta_x); // f <- f - delta_x / dt
-
- } else {
-
- ASSERT2(have_is_maps);
- // Some constraints
-
- Vec x_diff, x0_diff, delta_x_diff, f_diff;
- PetscCall(VecGetSubVector(x, is_diff, &x_diff));
- PetscCall(VecGetSubVector(x0, is_diff, &x0_diff));
- PetscCall(VecGetSubVector(delta_x, is_diff, &delta_x_diff));
- PetscCall(VecGetSubVector(f, is_diff, &f_diff));
-
- PetscCall(VecWAXPY(delta_x_diff, -1.0, x0_diff,
- x_diff)); // delta_x_diff = x_diff - x0_diff
- PetscCall(
- VecAXPY(f_diff, -1.0 / dt, delta_x_diff)); // f_diff <- f_diff - delta_x / dt
-
- PetscCall(VecRestoreSubVector(x, is_diff, &x_diff));
- PetscCall(VecRestoreSubVector(x0, is_diff, &x0_diff));
- PetscCall(VecRestoreSubVector(delta_x, is_diff, &delta_x_diff));
- PetscCall(VecRestoreSubVector(f, is_diff, &f_diff));
- }
+ PetscCall(performCalculation(
+ [this](Vec x, Vec x0, Vec delta_x, Vec f) -> PetscErrorCode {
+ PetscCall(VecWAXPY(delta_x, -1.0, x0, x)); // delta_x = x - x0
+ PetscCall(VecAXPY(f, -1.0 / dt, delta_x)); // f <- f - delta_x / dt
+ return PETSC_SUCCESS;
+ },
+ is_diff, x, x0, delta_x, f));
break;
}
case BoutSnesEquationForm::pseudo_transient: {
@@ -1618,34 +1615,15 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) {
// except that Δt is a vector
// F = (x0 - x)/Δt + f
// Algebraic: F = G(x) (already stored in f by rhs_function)
-
- if (!have_constraints) {
-
- VecWAXPY(delta_x, -1.0, x0, x);
- VecPointwiseDivide(delta_x, delta_x, dt_vec); // delta_x /= dt
- VecAXPY(f, -1.0, delta_x); // f <- f - delta_x
-
- } else {
- ASSERT2(have_is_maps);
-
- Vec x_diff, x0_diff, delta_x_diff, f_diff, dt_vec_diff;
- PetscCall(VecGetSubVector(x, is_diff, &x_diff));
- PetscCall(VecGetSubVector(x0, is_diff, &x0_diff));
- PetscCall(VecGetSubVector(delta_x, is_diff, &delta_x_diff));
- PetscCall(VecGetSubVector(f, is_diff, &f_diff));
- PetscCall(VecGetSubVector(dt_vec, is_diff, &dt_vec_diff));
-
- PetscCall(VecWAXPY(delta_x_diff, -1.0, x0_diff, x_diff));
- PetscCall(
- VecPointwiseDivide(delta_x_diff, delta_x_diff, dt_vec_diff)); // delta_x /= dt
- PetscCall(VecAXPY(f_diff, -1.0, delta_x_diff)); // f <- f - delta_x
-
- PetscCall(VecRestoreSubVector(delta_x, is_diff, &delta_x_diff));
- PetscCall(VecRestoreSubVector(x, is_diff, &x_diff));
- PetscCall(VecRestoreSubVector(x0, is_diff, &x0_diff));
- PetscCall(VecRestoreSubVector(f, is_diff, &f_diff));
- PetscCall(VecRestoreSubVector(dt_vec, is_diff, &dt_vec_diff));
- }
+ PetscCall(performCalculation(
+
+ [](Vec x, Vec x0, Vec delta_x, Vec f, Vec dt_vec) -> PetscErrorCode {
+ PetscCall(VecWAXPY(delta_x, -1.0, x0, x));
+ PetscCall(VecPointwiseDivide(delta_x, delta_x, dt_vec)); // delta_x /= dt
+ PetscCall(VecAXPY(f, -1.0, delta_x)); // f <- f - delta_x
+ return PETSC_SUCCESS;
+ },
+ is_diff, x, x0, delta_x, f, dt_vec));
break;
}
case BoutSnesEquationForm::backward_euler: {
@@ -1653,28 +1631,13 @@ PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) {
// 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));
- }
+ PetscCall(performCalculation(
+ [this](Vec x, Vec x0, Vec f) -> PetscErrorCode {
+ PetscCall(VecAYPX(f, -dt, x)); // f <- x - Δt*f
+ PetscCall(VecAXPY(f, -1.0, x0)); // f <- f - x0
+ return PETSC_SUCCESS;
+ },
+ is_diff, x, x0, f));
break;
}
case BoutSnesEquationForm::direct_newton: { |
2 participants
This is an extension of an older PR #3251 that was later merged in #3386
Probably we will need to merge it with #3386 before it can be merged to next