Merge SNES constraints into refactor#3386
Conversation
…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[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
|
|
||
| if (have_constraints) { | ||
| // Create PETSc-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) { | ||
| // Create PETSc-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) { | ||
| // Create PETSc-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, 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
|
|
||
| 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; | |
| iVec 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>
| .doc("Apply asinh() to all variables?") | ||
| .withDefault<bool>(false)) {} | ||
| .withDefault<bool>(false)) { | ||
| has_constraints = true; ///< This solver can handle constraints |
There was a problem hiding this comment.
Should this not default to false? Then L384 could be removed?
Also, maybe better to set the default in the declaration in the header
| 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; | ||
| } | ||
| } |
There was a problem hiding this comment.
We can use the shiny ranges algorithms now!
| 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; | |
| } | |
| } | |
| using std::ranges::any_of; | |
| has_constraints = any_of(f2d, &VarStr<Field2D>::constraint) | |
| or any_of(f3d, &VarStr<Field3D>::constraint); |
I think that pointer-to-member function syntax is right, you might need a lambda:
auto f_has_constraint = []<class T>(const VarStr<T>& f) { return f.constraint; };as an intermediate. This has the added advantage of using all the valid bracket types in one expression :)
| 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)); | ||
|
|
||
| if (equation_form != BoutSnesEquationForm::backward_euler) { | ||
| throw BoutException( | ||
| "SNES constraints currently require equation_form=backward_euler"); | ||
| } |
There was a problem hiding this comment.
Can we do the check first?
| 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)); | |
| if (equation_form != BoutSnesEquationForm::backward_euler) { | |
| throw BoutException( | |
| "SNES constraints currently require equation_form=backward_euler"); | |
| } | |
| if (equation_form != BoutSnesEquationForm::backward_euler) { | |
| throw BoutException( | |
| "SNES constraints currently require equation_form=backward_euler"); | |
| } | |
| 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)); |
3 participants
Merges the changes in #3251 into the refactor-petsc-precon branch.