Repository navigation
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.
|
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: { |
- Restored Jacobian-domain-error handling. - Removed the PID growth tweak. - Added cleanup for the constraint IS objects in the SNES destructor. - Updated manual page on SNES constraints - Added a new integrated test at `tests/integrated/test-snes-constraints` that tests convergence to a known solution with a constraint.
| } | ||
| }; | ||
|
|
||
| BOUTMAIN(TestSnesConstraints); |
There was a problem hiding this comment.
warning: variable 'init_err' of type 'int' can be declared 'const' [misc-const-correctness]
BOUTMAIN(TestSnesConstraints);
^Additional context
include/bout/physicsmodel.hxx:430: expanded from macro 'BOUTMAIN'
int init_err = BoutInitialise(argc, argv); \
^- Remove redundant `have_is_maps` flag. - Refactor common code suggested by @cmacmackin. - Move settings suggestions from comments into manual.
| PetscErrorCode ierr = PETSC_SUCCESS; | ||
| std::size_t acquired = 0; | ||
| for (; acquired < N; ++acquired) { | ||
| ierr = VecGetSubVector(vectors[acquired], indices, &subvectors[acquired]); |
There was a problem hiding this comment.
warning: do not use array subscript when the index is not an integer constant expression [cppcoreguidelines-pro-bounds-constant-array-index]
ierr = VecGetSubVector(vectors[acquired], indices, &subvectors[acquired]);
^| PetscErrorCode ierr = PETSC_SUCCESS; | ||
| std::size_t acquired = 0; | ||
| for (; acquired < N; ++acquired) { | ||
| ierr = VecGetSubVector(vectors[acquired], indices, &subvectors[acquired]); |
There was a problem hiding this comment.
warning: do not use array subscript when the index is not an integer constant expression [cppcoreguidelines-pro-bounds-constant-array-index]
ierr = VecGetSubVector(vectors[acquired], indices, &subvectors[acquired]);
^| } | ||
|
|
||
| for (std::size_t i = acquired; i > 0; --i) { | ||
| PetscErrorCode restore_ierr = |
There was a problem hiding this comment.
warning: variable 'restore_ierr' of type 'PetscErrorCode' (aka 'int') can be declared 'const' [misc-const-correctness]
| PetscErrorCode restore_ierr = | |
| PetscErrorCode const restore_ierr = |
|
|
||
| for (std::size_t i = acquired; i > 0; --i) { | ||
| PetscErrorCode restore_ierr = | ||
| VecRestoreSubVector(vectors[i - 1], indices, &subvectors[i - 1]); |
There was a problem hiding this comment.
warning: do not use array subscript when the index is not an integer constant expression [cppcoreguidelines-pro-bounds-constant-array-index]
VecRestoreSubVector(vectors[i - 1], indices, &subvectors[i - 1]);
^|
|
||
| for (std::size_t i = acquired; i > 0; --i) { | ||
| PetscErrorCode restore_ierr = | ||
| VecRestoreSubVector(vectors[i - 1], indices, &subvectors[i - 1]); |
There was a problem hiding this comment.
warning: do not use array subscript when the index is not an integer constant expression [cppcoreguidelines-pro-bounds-constant-array-index]
VecRestoreSubVector(vectors[i - 1], indices, &subvectors[i - 1]);
^
|
Thanks for your review and suggestions @cmacmackin ! I think I have addressed them. I also merged the changes in #3386 into this PR. I think it's ready to go. |
bendudson
added
feature
Comments to explain the expected solution. Use boututils launch to run example
For constraints it may be clearer to write `residual(var)` rather than `ddt(var)`
Protected variable in Solver, renamed for clarity. If a model has algebraic constraints then the solver sets a variable `has_constraint_variables`.
|
Thanks @cmacmackin ! I've renamed |
Use residual(phi) rather than ddt(phi) in example.
|
|
||
| /// Can this solver handle constraints? Set to true if so. | ||
| bool has_constraints{false}; | ||
| bool supports_constraints{false}; |
There was a problem hiding this comment.
warning: member variable 'supports_constraints' has protected visibility [cppcoreguidelines-non-private-member-variables-in-classes]
bool supports_constraints{false};
^
|
Looks good, except that I see you didn't add a |
Fixes unit tests
|
Hi @cmacmackin ! |
Ah, that would be a good reason! |
3 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