Fix LaplacePetsc3dAmg #2237
Merged
Fix LaplacePetsc3dAmg #2237
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
These are occasionally needed, and it only adds initialisation cost to calculate them.
The method actually creates the intersection. Fixes #2236.
getUnion() was misleadingly named, and actually gave the intersection. Region::operator+=(Region) updates a Region to be the union of two Regions, which is almost what was wanted, although we have to call mask_region.unique() afterwards too to ensure there are no repeated entries and the entries are sorted.
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
Suggested by clang-tidy.
|
Only the Fedora tests are failing, but they're only failing on the test that this PR is supposed to fix. It says a PETSc matrix element argument was out of range, but I've tested with PETSC-3.14.4 (same version as Fedora) and an identical |
|
It is a bit concerning, and I'd prefer to fix it if at all possible. I'll see if I can reproduce it locally. |
|
I can reproduce it using the fedora packages. petsc is compiled with (mpich) openmpi: |
|
I finally managed to reproduce it by compiling with optimisations, the debug build doesn't fail. Unfortunately, this means the debugger is not much help, as lots of stuff has been inlined and otherwise rearranged. Here's the (abridged) output with some extra debugging stuff put in: The code that's causing the problem is: const int yup = (l.y() == localmesh->yend && upperY.intersects(l.x())) ? -1 : 0,
ydown = (l.y() == localmesh->ystart && lowerY.intersects(l.x())) ? -1 : 0;
operator3D.yup(yup)(l, l.yp()) = 0.0;
operator3D.ydown(ydown)(l, l.ym()) = 0.0;
operator3D.yup(yup)(l, l.xp().yp()) = 0.0;
operator3D.ydown(ydown)(l, l.xp().ym()) = 0.0;
operator3D.yup(yup)(l, l.xm().yp()) = 0.0; // <---- this line
operator3D.ydown(ydown)(l, l.xm().ym()) = 0.0;
operator3D.yup(yup)(l, l.yp().zp()) = 0.0;
operator3D.yup(yup)(l, l.yp().zm()) = 0.0;
operator3D.ydown(ydown)(l, l.ym().zp()) = 0.0;
operator3D.ydown(ydown)(l, l.ym().zm()) = 0.0;
}
operator3D.partialAssemble();which is a bit confusing, because we clearly haven't assembled the matrix yet, so I don't know why it would complain about mallocs. We do call I'm going to try using ASAN/valgrind, see either of them get anywhere |
|
That's bizzare - the preallocation and matrix elements that we set must be the same as in builds that are not raising errors?!? I'd not thought of the debugging mode in PETSc, everything I've tested has been with optimized PETSc builds. @ZedThree what configure command did you use exactly to reproduce the error? |
|
This is my PETSc configure command: and BOUT++: This diff: diff --git a/include/bout/petsc_interface.hxx b/include/bout/petsc_interface.hxx
index 1d1f3dfdf..a44fe5ea1 100644
--- a/include/bout/petsc_interface.hxx
+++ b/include/bout/petsc_interface.hxx
@@ -617,16 +617,19 @@ public:
}
}
Element& operator=(const Element& other) {
+ AUTO_TRACE();
ASSERT3(finite(static_cast<BoutReal>(other)));
return *this = static_cast<BoutReal>(other);
}
Element& operator=(BoutReal val) {
+ AUTO_TRACE();
ASSERT3(finite(val));
value = val;
setValues(val, INSERT_VALUES);
return *this;
}
Element& operator+=(BoutReal val) {
+ AUTO_TRACE();
ASSERT3(finite(val));
auto columnPosition = std::find(positions.begin(), positions.end(), petscCol);
if (columnPosition != positions.end()) {
@@ -641,10 +644,12 @@ public:
private:
void setValues(BoutReal val, InsertMode mode) {
+ TRACE("PetscMatrix setting values at ({}, {})", petscRow, petscCol);
ASSERT3(positions.size() > 0);
std::vector<PetscScalar> values;
std::transform(weights.begin(), weights.end(), std::back_inserter(values),
[&val](BoutReal weight) -> PetscScalar { return weight * val; });
+
int status;
BOUT_OMP(critical)
status = MatSetValues(*petscMatrix, 1, &petscRow, positions.size(),
@@ -724,7 +729,7 @@ public:
BOUT_OMP(critical)
status = MatGetValues(*get(), 1, &global1, 1, &global2, &value);
if (status != 0) {
- throw BoutException("Error when setting elements of a PETSc matrix.");
+ throw BoutException("Error when getting elements of a PETSc matrix.");
}
return value;
}
diff --git a/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx b/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx
index fdb65599a..a8a60387d 100644
--- a/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx
+++ b/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx
@@ -178,6 +178,7 @@ LaplacePetsc3dAmg::~LaplacePetsc3dAmg() {
Field3D LaplacePetsc3dAmg::solve(const Field3D &b_in, const Field3D &x0) {
+ AUTO_TRACE();
// If necessary, update the values in the matrix operator and initialise
// the Krylov solver
if (updateRequired) updateMatrix3D();and using valgrind had a ton of output, but all leaks, rather than invalid read/writes. |
|
AddressSanitizer to the rescue! Fix incoming |
Stencil was storing a lambda which captured some state by reference, which went out of scope before the lambda was called. Fix is to capture by value. `RangeIterator` is not particularly expensive to copy, so this should not affect performance. Possible issue with linked `RangeIterator`s?
|
Fix was basically a single character: modified src/invert/laplace/impls/petsc3damg/petsc3damg.cxx
@@ -511,7 +511,7 @@ OperatorStencil<Ind3D> LaplacePetsc3dAmg::getStencil(Mesh* localmesh,
// If there is a lower y-boundary then create a part of the stencil
// for cells immediately adjacent to it.
if (lowerYBound.max() - lowerYBound.min() > 0) {
- stencil.add([index = localmesh->ystart, &lowerYBound](Ind3D ind) -> bool {
+ stencil.add([index = localmesh->ystart, lowerYBound](Ind3D ind) -> bool {
return index == ind.y() && lowerYBound.intersects(ind.x()); },
lowerEdgeStencilVector);
}
@@ -519,7 +519,7 @@ OperatorStencil<Ind3D> LaplacePetsc3dAmg::getStencil(Mesh* localmesh,
// If there is an upper y-boundary then create a part of the stencil
// for cells immediately adjacent to it.
if (upperYBound.max() - upperYBound.min() > 0) {
- stencil.add([index = localmesh->yend, &upperYBound](Ind3D ind) -> bool {
+ stencil.add([index = localmesh->yend, upperYBound](Ind3D ind) -> bool {
return index == ind.y() && upperYBound.intersects(ind.x()); },
upperEdgeStencilVector);
}Basically the stencil test function was capturing the I'm not sure why this worked in all the other builds! |
ZedThree
approved these changes
Feb 17, 2021
|
Thanks @ZedThree! Fix looks good to me. |
|
|
This was referenced Mar 22, 2021
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
The bug in
test-laplace-petsc3d that @ZedThree found (#2225) was hiding a couple of other bugs inZInterpolationandZHermiteSpline` that made the matrix for 3d Laplace_perp inversion incorrect with non-trivial parallel derivatives. These bugs should be fixed by this PR.I've also added some extra checks (when
CHECK > 2) for finiteness of matrix elements and interpolation weights, and bounds of indices passed togetWeightsForYApproximation(), which helped to find the bugs.Fixes #2225.