Fix edge case where meshblocks can be improperly de-refined#1248
Conversation
|
@pgrete any idea what's going on with the CI here? This looks like a local issue with the CI machine? |
| @@ -1,3 +1,4 @@ | |||
| #!/usr/bin/env python | |||
There was a problem hiding this comment.
This shebang was just missing before.
|
Thanks for the additional cleanup @lroberts36 ! Looks like builds are now failing? |
For whatever reason some lock file was not cleaned up. A manual removal seems to have fixed the issue. |
Yurlungur
left a comment
There was a problem hiding this comment.
Thanks @lroberts36 for the additional cleanup. Original bug is below. Luke's changes are non-functional cleanup that better isolate the tree machinery to prevent something like this from being re-introduced.
| if (r < tnderef) { | ||
| if ((lderef[n].lx1() + i) == lderef[r].lx1() && | ||
| (lderef[n].lx2() + j) == lderef[r].lx2() && | ||
| (lderef[n].lx3() + k) == lderef[r].lx3() && | ||
| lderef[n].level() == lderef[r].level()) | ||
| if (lderef[n].GetSameLevelNeighbor(oi, oj, ok) == lderef[r]) { | ||
| rr++; | ||
| } |
There was a problem hiding this comment.
this change here was the location of the bug.
|
@pgrete this error looks like the CI machine didn't see any GPUs? |
pgrete
left a comment
There was a problem hiding this comment.
LGTM (as far as I was able to follow).
| for (std::int64_t ok = 0; ok <= lk; ok++) { | ||
| for (std::int64_t oj = 0; oj <= lj; oj++) { | ||
| for (std::int64_t oi = 0; oi <= 1; oi++) { |
There was a problem hiding this comment.
I cannot say that I have fully internalized the updated mesh structure yet, but double checking for my own understanding:
Why are lk and lj tied to symmetry but i is not?
I see that our grids are always at least 1D (which I'd understand overall a different behavior), but what determines the use of symmetry?
IIRC I had a similar question in the past and it was tied to being periodic but I'm not that sure any more.
There was a problem hiding this comment.
@lroberts36 should jump in but... as I understand it, symmetry is totally synonymous with dimension. It's essentially, "are we plane-symmetric in this direction?"
There was a problem hiding this comment.
Yes, symmetry(dir) being true means that we have translational symmetry in direction dir. Previously in many places in the code we just checked if the size of a block was one in a given direction, but this created some issues with multigrid when multiple block sizes were allowed. I think we could also probably remove the symmetry stuff and replace it with ndim checks.
There was a problem hiding this comment.
I prefer checking symmetry. It's more self-descriptive.
There was a problem hiding this comment.
Maybe it should be changed to TranslationalSymmetry or something because this isn't the first time someone has been confused about what that call implies.
There was a problem hiding this comment.
A reasonable suggestion for a later MR.
There was a problem hiding this comment.
yeah, when I wrote "Maybe it should be changed..." I was thinking the change should occur at some (possibly infinite) time in the future
I think the machine is currently under heavy use. I just retriggerd the job (expecting that it'll go through). |
Co-authored-by: Philipp Grete <pgrete@hs.uni-hamburg.de>
Co-authored-by: Philipp Grete <pgrete@hs.uni-hamburg.de>
PR Summary
The last week or so, I've been chasing down an AMR bug in riot. With a Rayleigh-Taylor instability, I was seeing the solution become asymmetric. I chased down the problem to the mesh refining incorrectly. Here's a plot of asymmetry in the fluid solution, with the meshblocks labeled by gid right before the problem appears:



On the next time step, a collection of blocks at the top derefines, while the symmetric block at the bottom does not. In particular, blocks 87, 88, 89, and 90 get coarsened, while blocks 2, 3, 4, and 5 do not.
One step later, the mesh re-refines:
At this point, the top blocks have experienced restriction and prolongation, while the bottom blocks have not. This introduces asymmetric truncation error.
The problem turned out to be in
mesh-amr_loadbalance.cpp. The code should only de-refine if all children of a parent block want to de-refine. However, the code that checks for this is quite brittle. It compares logical locations, rather than actual parent associativity. When we moved to a forest of octrees, rather than a single tree, we needed to add a comparison for tree as well as location within a tree, and that check was missing.In this MR, I move to a slightly less brittle approach suggested by @lroberts36 which checks for "same level neighbor" rather than computing by hand. That said, this code is still quite brittle and we may wish to move away from the current logic, which assumes morton ordering, to logic that just tracks which blocks want to de-refine via, e.g., an unordered set.
While I was at it, I also couldn't resist cleaning up the memory management, but replacing the
newanddeletecalls with astd::vectorand adding some explanatory comments.I also added a little documentation and added the script I used to check for symmetry.
PR Checklist