Lean's boundary condition setting for BSSN constraints only fills in a single ghostzone
In developing BaikalETK (a NRPy+-based BSSN thorn), I set up a schedule.ccl file that follows almost exactly the same pattern as Lean.
I noticed that when BaikalETK evaluated the constraints without initializing the variables to zero (as is done in Lean), I would get uninitialized on the boundaries **after ApplyBCs was called**. This confused me since BaikalETK, like Lean, calls for "flat" boundary conditions on the constraints. Sure enough when I performed the same run with Lean, most boundary points were filled with zeros (unlike BaikalETK, the constraints in Lean are initialized to zero). I spoke to Roland (Haas) about this, and he was also confused until he noticed that when BaikalETK/Lean calls Boundary_SelectVarForBC/Boundary_SelectGroupForBC a la:
Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, +1, -1, "BaikalETK::somegridfunction", "flat") ,
the "+1" indicates that only one ghost zone point should be filled in. You might then wonder what would happen if we just replaced that "+1" with "cctk_nghostzones[0]". Well, since this function in Lean & BaikalETK is called in "LEVEL" mode, "cctk_nghostzones[0]" does not exist (actually it compiles and gives a garbage value).
Update (Jun 16, 2020): The cctk_nghostzones[0] idea in the previous paragraph (combined with calling the function in LOCAL) mode results in issues when certain numbers of MPI processes are used (thanks to Tanmayee Gupte for pointing this out!) I ended up using Erik’s suggestion (run in LEVEL mode) below. It's currently implemented in BaikalETK (here's the schedule.ccl: https://nbviewer.jupyter.org/github/zachetienne/nrpytutorial/blob/master/Tutorial-BaikalETK.ipynb#scheduleccl) and seems to work great.
Comments (37)
-
-
reporter - marked as bug
- edited description
-
reporter - edited description
-
For the sake of documenting this: the value on
cctk_nghostznones
in LEVEL mode (Carpet only, the issue does not arise in PUGH) is well defined (ie not technically “garbage” as in “unitinialized memory” or “to be garbage collected”) and is 666 (by default) which is chosen to be a recognizable value that is likely to trigger errors when used.See https://bitbucket.org/eschnett/carpet/src/master/Carpet/src/modes.cc#lines-528 and
deadbeef
is a CarpetLib parameter defaulting to 666. -
- marked as major
-
@helvi witek @Miguel Zilhão the way boundary conditions are currently done in Lean_public means that only 1 point is filled in, however for a stencil with of (say) 3 in Lean_public this means that 2 points at the boundary will be left unitialized in the same way they were for Baikal.
This should be fixed in the same way @Zach Etienne fixed it in Baikal, see
#2432. -
i’m trying to fix this in the development version of Lean, but i’m getting very strange warnings. i’ve merely added the following call to the function GetBoundarySizesAndTypes on
Boundaries.F90
CCTK_INT bndsize(6), is_ghostbnd(6), is_symbnd(6), is_physbnd(6) ierr = GetBoundarySizesAndTypes( cctkGH, 6, bndsize, is_ghostbnd, is_symbnd, is_physbnd )
i don’t use any of the output; i merely call the function. when running a sample parameter file, i’m then showered with a bunch of warnings:
WARNING level 1 from host majikthise process 3
while executing schedule bin MoL_PseudoEvolution, routine LeanBSSNMoL::LeanBSSN_Constraints_Bounda
ries
in thorn CoordBase, file /home/mzilhao/dev/ET/Cactus/arrangements/CactusBase/CoordBase/src/Domain.
c:449:
-> The number of symmetry points in direction 0 face 0 is 3, which is different from the number of
ghost points (666); this is probably an error
WARNING level 1 from host majikthise process 3
while executing schedule bin MoL_PseudoEvolution, routine LeanBSSNMoL::LeanBSSN_Constraints_Bounda
ries
in thorn CoordBase, file /home/mzilhao/dev/ET/Cactus/arrangements/CactusBase/CoordBase/src/Domain.
c:449:
-> The number of symmetry points in direction 1 face 0 is 3, which is different from the number of
ghost points (666); this is probably an error
WARNING level 1 from host majikthise process 3
while executing schedule bin MoL_PseudoEvolution, routine LeanBSSNMoL::LeanBSSN_Constraints_Bounda
ries
in thorn CoordBase, file /home/mzilhao/dev/ET/Cactus/arrangements/CactusBase/CoordBase/src/Domain.
c:449:
-> The number of symmetry points in direction 2 face 0 is 3, which is different from the number of
ghost points (666); this is probably an error
WARNING level 1 from host majikthise process 5
while executing schedule bin MoL_PseudoEvolution, routine LeanBSSNMoL::LeanBSSN_Constraints_Bounda
ries
i’ve never seen this before, and it goes away if remove the call to
GetBoundarySizesAndTypes
. again, note that i’m not doing anything with its output -- i just call the function. what could be the issue?
-
- attached LeanBSSN_Schw_hf64.par
- attached BCs.patch
patch for the change that gives problems and sample parfile to run
-
The function to use is actually
GetBoundarySpecification
notGetBoundarySizesAndTypes
. Easiest is most likely to copy the code from ML_ADMConstraint’sKranc.cc
andML_ADMConstraints_evaluate.cc
files, ie:/********************************************************************* * GetBoundaryWidths *********************************************************************/ void GetBoundaryWidths(cGH const * restrict const cctkGH, CCTK_INT nboundaryzones[6]) { CCTK_INT is_internal[6]; CCTK_INT is_staggered[6]; CCTK_INT shiftout[6]; int ierr = -1; if (CCTK_IsFunctionAliased ("MultiPatch_GetBoundarySpecification")) { int const map = MultiPatch_GetMap (cctkGH); /* This doesn't make sense in level mode */ if (map < 0) { static int have_warned = 0; if (!have_warned) { CCTK_WARN(1, "GetBoundaryWidths: Could not determine current map (can be caused by calling in LEVEL mode)"); have_warned = 1; } for (int i = 0; i < 6; i++) nboundaryzones[i] = 0; return; } ierr = MultiPatch_GetBoundarySpecification (map, 6, nboundaryzones, is_internal, is_staggered, shiftout); if (ierr != 0) CCTK_WARN(0, "Could not obtain boundary specification"); } else if (CCTK_IsFunctionAliased ("GetBoundarySpecification")) { ierr = GetBoundarySpecification (6, nboundaryzones, is_internal, is_staggered, shiftout); if (ierr != 0) CCTK_WARN(0, "Could not obtain boundary specification"); } else { CCTK_WARN(0, "Could not obtain boundary specification"); } } /********************************************************************* * GetBoundaryWidth *********************************************************************/ int GetBoundaryWidth(cGH const * restrict const cctkGH) { CCTK_INT nboundaryzones[6]; GetBoundaryWidths(cctkGH, nboundaryzones); int bw = nboundaryzones[0]; for (int i = 1; i < 6; i++) if (nboundaryzones[i] != bw) CCTK_WARN(0, "Number of boundary points is different on different faces"); return bw; }
and
extern "C" void ML_ADMConstraints_evaluate_SelectBCs(CCTK_ARGUMENTS) { #ifdef DECLARE_CCTK_ARGUMENTS_ML_ADMConstraints_evaluate_SelectBCs DECLARE_CCTK_ARGUMENTS_CHECKED(ML_ADMConstraints_evaluate_SelectBCs); #else DECLARE_CCTK_ARGUMENTS; #endif DECLARE_CCTK_PARAMETERS; if (cctk_iteration % ML_ADMConstraints_evaluate_calc_every != ML_ADMConstraints_evaluate_calc_offset) return; CCTK_INT ierr CCTK_ATTRIBUTE_UNUSED = 0; ierr = KrancBdy_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADMConstraints::ML_Ham","flat"); if (ierr < 0) CCTK_WARN(CCTK_WARN_ALERT, "Failed to register flat BC for ML_ADMConstraints::ML_Ham."); ierr = KrancBdy_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADMConstraints::ML_mom","flat"); if (ierr < 0) CCTK_WARN(CCTK_WARN_ALERT, "Failed to register flat BC for ML_ADMConstraints::ML_mom."); return; }
which you could either translate to Fortran (not hard) or keep as C functions and decorate with some
CCTK_FNAME
then call from Fortran (with the usual issues such scGH* cctkGH
turning intocGH ** cctkGH
.Note:
KrancBdy_SelectGroupForBC
is just a wrapper aroundBoundary_SelectGroupForBC
with some extra code for the PreSync functionality in Cactus. -
hmm, but the code on Baikal that i was copying from does use
GetBoundarySizesAndTypes
: https://bitbucket.org/zach_etienne/wvuthorns/src/472cc7f3a8a860040323bb52f5ba1e4e2f4df26d/Baikal/src/BoundaryConditions.c (line 102)i was trying to base the implementation on that, since @Zach Etienne reports it’s working well…
-
Though looking at the code that is used for the
MultiPatch_GetBoundarySpecification
case, this seems incorrect as well. If called in LEVEL mode it will set all boundary widths to 0. And when called in SINGLEMAP mode (so thatmap
is defined) you would end up selecting boundary conditions more than once which makesBoundary
abort.Something similar to what is done when going from
GetBoundaryWidths
toGetBoundaryWidth
seems required, ie.:if (CCTK_IsFunctionAliased ("MultiPatch_GetBoundarySpecification")) { int const maps = MultiPatch_GetMaps (cctkGH); for (int i = 0; i < 6; i++) nboundaryzones[i] = -1; // collect all boundary information from all maps, and get // outer boundary information only for (int m = 0; m < maps; m++) { CCTK_INT nbounds[6]; ierr = MultiPatch_GetBoundarySpecification (m, 6, nboundaryzones, is_internal, is_staggered, shiftout); if (ierr != 0) CCTK_WARN(0, "Could not obtain boundary specification"); for(int i = 0; i < 6; i++) { if(!is_internal[i]) { if(nboundaryzones[i] == -1) { nboundaryzones[i] = nbounds[i]; } else { if(nboundaryzones[i] != nbounds[i]) { CCTK_WARN(0, "Inconsistent number of outer boundary points: %d != %d in the %s %c direction between map %d and an earlier map", nboundaryzones[i], nbounds[i], i % 2 ? "upper" : "lower", "xyz"[i/2], m); } else { // identical boundary sizes, all is fine } } // if nboundaryzones } // if is_internal } // for i } // for m for (int i = 0; i < 6; i++) { assert(nboundaryzones[i] != -1); } } // if multipatch
note that this is utterly untested.
-
Zach may be ignoring the level 1 warning. For this usage case, he can possibly do so (though I am not sure since the results of
CoordBase_GetBoundarySizesAndTypes
will differ on different MPI ranks since it returns the ghost size for patch boundaries that are not physical boundaries but instead inter-processor boundaries.
@Erik Schnetter @Zach Etienne you may need to clarify the suggested code and observed behaviour. -
reporter @Miguel Zilhão : Baikal does not yet officially support symmetries, hence I didn’t see the same warning (regarding symmetry points) you did. In fact I didn’t see any warnings whatsoever.
-
When you ran Cactus did you pass it options to show the warnings? Ie
-L3 -W3
? -
I filed a Kranc ticket about this in https://github.com/ianhinder/Kranc/issues/150
-
shouldn’t there also be a bug with
GetBoundarySizesAndTypes
? at the very least, its name is misleading, since i wouldn’t expect a getter function to change any values, and this example seems to show that it acts non-trivially when called… -
@Miguel Zilhão it does not change any value. It is just its returned values will differ on different MPI ranks since a given boundary (of a grid patch) may be of a different type on different MPI ranks, eg in a situation with 2 ranks and a grid split like this:
| P1 | P2 | x->
then for process 1 the left x-boundary is an outer boundary, but for process 2 it is an interprocessor boundary. Thus
GetBoundarySizesAndTypes
will correctly report that fact on each patch.The point is that when selecting boundary conditions you do not want a per-patch information, but a per-level information (this is the reason for it being n LEVEL mode and not LOCAL mode).
The
GetBoundarySpecification
on the other hand is not concerned with per-patch information and just returns the values that were set in CoordBase’s parameters (repos/cactusbase/CoordBase/src/Domain.c):CCTK_INT CoordBase_GetBoundarySpecification(CCTK_INT const size, CCTK_INT *const nboundaryzones, CCTK_INT *const is_internal, CCTK_INT *const is_staggered, CCTK_INT *const shiftout) { DECLARE_CCTK_PARAMETERS; if (!(size >= 0)) CCTK_ERROR("size is less than zero"); if (!(nboundaryzones)) CCTK_ERROR("nboundaryzones is out of bounds"); if (!(is_internal)) CCTK_ERROR("is_internal is out of bounds"); if (!(is_staggered)) CCTK_ERROR("is_staggered is out of bounds"); if (!(shiftout)) CCTK_ERROR("shiftout is out of bounds"); if (!(size == 6)) CCTK_ERROR("size is out of bounds"); nboundaryzones[0] = boundary_size_x_lower; nboundaryzones[1] = boundary_size_x_upper; nboundaryzones[2] = boundary_size_y_lower; nboundaryzones[3] = boundary_size_y_upper; nboundaryzones[4] = boundary_size_z_lower; nboundaryzones[5] = boundary_size_z_upper;
-
@Zach Etienne I think I understand why you are not seeing a warning. If one looks at the code in CoordBase’s
CoordBase_GetBoundarySizesAndTypes
function then the interesting lines are:421 /* Determine boundary types */ 422 for (int d = 0; d < 2 * dim; ++d) { 423 /* Ghost boundaries (inter-process or mesh refinement boundaries) 424 have cctk_bbox=0 */ 425 is_ghostbnd[d] = !cctkGH->cctk_bbox[d]; 426 /* Symmetry boundaries are not ghost boundaries, are described as 427 symmetry boundaries, and are not multi-patch outer 428 boundaries */ 429 is_symbnd[d] = !is_ghostbnd[d] && symbnd[d] >= 0 && !mp_bbox[d]; 430 /* Physical (outer) boundaries are all other boundaries */ 431 is_physbnd[d] = !is_ghostbnd[d] && !is_symbnd[d]; 432 } 433 434 /* Determine boundary sizes */ 435 for (int dir = 0; dir < dim; ++dir) { 436 for (int face = 0; face < 2; ++face) { 437 438 if (is_ghostbnd[2 * dir + face]) { 439 /* Ghost boundary */ 440 bndsize[2 * dir + face] = cctkGH->cctk_nghostzones[dir]; 441 } else { 442 /* Symmetry or physical boundary */ 443 bndsize[2 * dir + face] = nboundaryzones[2 * dir + face];
where lines 440 and 443 set the elements of the returned array to the number of ghost zones or number of boundary zones depending on the
is_ghostbnd
array. Howeveris_ghostbnd
is set as!cctkGH->cctk_bbox[d]
in line 425. InLEVEL
mode Carpet will setcctk_bbox
to the deadbeef value of 666 which is “true” to C so!cctkGH->cctk_bbox[d]
is 0 and all faces are claimed to be outer boundaries (or symmetry boundaries).While this “works” it is still incorrect use of the function since this depends on the value of
deadbeef
which is really a poison like value and can be set at runtime viatCarpetLib::deadbeef
(though 0 is not a usable value so I cannot demonstrate the issue any more than I already have).Thus, as far as I can tell, using
GetBoundarySizesAndTypes
is incorrect since the function is intended to be used inLOCAL
mode. My guess would be that Erik misspoke and meant to suggestGetBoundarySpecification
for all cases (though note my issues withMultiPatch_GetBoundarySpecification
having to be called inSINGLEMAP
mode which is belowLEVEL
mode ie would cause multiple calls toBoundary_SelectGroupForBCs
). -
reporter @Roland Haas :Thanks for the explanation!
It is still evident that the backported version, despite depending on
deadbeef
within anif()
behaves “more correctly” than the unpatched version. I agree that the backport version isnt “correct” but it’s “correct-er” than the released version. That being said, I’d like to see this fixed properly, and after reading the above, I am left confused as to what is the correct means to fix. Should I useGetBoundarySpecification
, some multipatch version, or write something from scratch? -
Since Baikal does not support Llama (it does not, does it?) all you should have to do is use
GetBoundarySpecification
instead ofGetBoundarySizesAndTypes
. Lean_public which does, or at least there were plans to do so, support Llama has to do more work.Ie this should work (and does for me using repos/wvuthorns/BaikalVacuum/test/BaikalVacuum_EE_O8_sgw3d.par even when increasing resolution to dx=0.0.1 and 12 MPI ranks so that not all MPI ranks have a boundary face all the time) :
diff --git a/BaikalVacuum/src/BoundaryConditions.c b/BaikalVacuum/src/BoundaryConditions.c index 7240e5b..f89a9d2 100644 --- a/BaikalVacuum/src/BoundaryConditions.c +++ b/BaikalVacuum/src/BoundaryConditions.c @@ -96,11 +96,11 @@ void BaikalVacuum_BoundaryConditions_aux_gfs(CCTK_ARGUMENTS) { CCTK_INT ierr CCTK_ATTRIBUTE_UNUSED = 0; CCTK_INT bndsize[6]; - CCTK_INT is_ghostbnd[6]; - CCTK_INT is_symbnd[6]; - CCTK_INT is_physbnd[6]; - GetBoundarySizesAndTypes(cctkGH, 6, bndsize, is_ghostbnd, is_symbnd, is_physbnd); + CCTK_INT is_internal[6]; + CCTK_INT is_staggered[6]; + CCTK_INT shiftout[6]; + GetBoundarySpecification(6, bndsize, is_internal, is_staggered, shiftout); ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, bndsize[0], -1, "BaikalVacuum::HGF", "flat"); if (ierr < 0) CCTK_ERROR("Failed to register BC for BaikalVacuum::HGF!");
-
i’ve checked Roland’s proposal above and it indeed seems to work fine for non-Llama evolutions. the current development version of Lean does support Llama, though, so it’d be good to cover that as well. i see that there is also the function
MultiPatch_GetBoundarySpecification
, which i guess could be used when Llama is used? -
- changed status to open
-
Yes,
MultiPatch_GetBoundarySpecification
is the equivalent function in Llama.However Llama’s function assumes to be called in
SINGLEMAP
mode (or more correctly it assumes being passed a “map” number) which is not available inLEVEL
mode. See https://bitbucket.org/einsteintoolkit/tickets/issues/2288/leans-boundary-condition-setting-for-bssn#comment-58230040 for what I believe to be the correct behaviour. Unfortunately no code that I am aware of currently handles this correctly. Not even McLachlan or the original CTGamma code, so input from the Llama developers would be helpful. -
i see. what would be the best way of handling this issue, in the meantime? what happens if i just call
GetBoundarySpecification
even when Llama is used? or is it preferable to leave things as is until we have input from the Llama developers? -
Calling
GetBoundarySpecification
in Llama code will give you incorrect answers / will fail. The function is provided by the thornCoordBase
which may not be active at all for a Llama run. if it is active the it will report its ow parameter settings and notCoordinates
' (the Llama equivalent thorn toCoordBase
).How to handle this in the developer’s branch is mostly up to you. You could give my proposed solution a try, you could also try and schedule the
SelectBCs
function asSINGLEMAP
then callMultiPatch_GetBoundarySpecifiction
and select boundary conditions on only those faces thatis_internal
is 0.Note that, while this would let you have different boundary widths for different boundary faces, you still cannot have the same face be a boundary on different maps and have different options (eg have 3 point wide boundary in “z” on map 1 but a 1 point wide boundary in “z” on map 2) due to the way theBoundary
thorn’sApplyBCs
group works. None of the current Llama setups requires this though. Eg for the, almost exclusively used, 7patch system the outer boundaries are all on maps 1-6 (the cubed spheres) and they all have the only outer boundary face being in direction of positive third coordinate (“radial”). Ie the coordinate setup is such that having a single “face” be selected for the outer boundary on all maps and that not eg the “radial” direction is a different coordinate number on the different maps.In Lean_public’s case backporting this into the release branch is not really required (though you can if you want to) since, different from Baikal, the current released code does not prevent you from running in certain situations and does give correct results (I think) when using the typical radiation boundary condition (same as McLachlan).
-
i've just committed a possible fix for this in the development version of Lean. i went with an alternative solution of "hardcoding" the boundary size according to which finite difference stencil is used. i hope this is ok.
by the way, this is probably a good time to ask if it'd be ok to merge our current development version onto master? we've been using the development branch for quite some time internally and it has been working fine. i should also mention that this branch has been cross-checked internally by different people.
-
Hell Miguel,
yes, tying the request boundary size to your codes stencil rather than the boundary size in the parfile is fine (done eg by GRHydro). In some sense better (you are actually asking for what you need), in some sense worse (there may not be enough boundary points in the parfile, or if there are too many then some will be left uninitialized).
I am just about to create new point release for ET_2020_05 so if your fixes are in master than I can included them (only the fixes though, no new features).
Since you are hosting Lean_public in your own repository and not one owned by the ET, there is very little that the ET (can) enforce as far as your group handling commits goes. If you are internally developing in in a “development” branch and then, after some internal review, you merge to master then you are already doing fine as far as I am concerned. “master” in the ET is explicitly marked as “this will break your code occasionally” (ie closer to a typical “development” branch).
If there is too much “bad” stuff in master at the time of the next release then the release branch will not update to master and stay stuck at the last release. If master becomes too much broken (or ends up having code that is deemed inappropriate for the ET) then we will either cherry pick from master what we want (lots of work for me so this will only happen if there is manifest interest in the code), not update the release branch (annoying if there are bugfixes), or fork the repo and run a separate ET fork (this would qualify as the “nuclear option” as it is most time consuming) essentially having a “stable” master and a “devel” master (happened to GRHydro vs. Zelmani for a while). Before that happens you can expect emails from the maintainers about “master” being a mess so there should be time to avoid this happening.
-
Pull request: https://bitbucket.org/canuda/lean_public/pull-requests/3 would take care of this I believe.
-
Please review (code is by Miguel and Helvi).
-
@helvi witek , @Miguel Zilhão , @Zach Etienne will review.
-
@helvi witek @Miguel Zilhão
@Miguel Zilhão
@Zach Etienne any progress on this?
-
The pull request https://bitbucket.org/canuda/lean_public/pull-requests/3 should fix this. Would this PR be ok to merge?
-
That is up to the reviewers i.e. Helvi, you, Zach to decide. If you have reached a decision that it can be applied then you should go ahead. I am only poking people but not reviewing.
-
I’ve added a comment about it in the PR above.
-
@helvi witek @Miguel Zilhão @Zach Etienne any news on this? Note that this only needs to be “approved” by your group and then applied. It should however be applied before the feature freeze for the ET_2020_11 release on November 12th. Please see https://docs.einsteintoolkit.org/et-docs/Release_Details#Schedule_for_ET_2020_11 for the planned release timeline for ET_2020_11.
-
I’ve merged PR https://bitbucket.org/canuda/lean_public/pull-requests/3 which takes care of this.
-
- changed status to resolved
Thank you.
- Log in to comment
The proper way to handle this is to inquire with thorn CoordBase how many boundary points there are via its aliased functions. There is e.g. “GetBoundarySizesAndTypes”, or (for multi-block simulations) “MultiPatch_GetBoundarySpecification”.
-erik