Lean's boundary condition setting for BSSN constraints only fills in a single ghostzone

[Zach Etienne]

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.

[Erik Schnetter]

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

[Roland Haas]

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.

[Roland Haas]

@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.

[Miguel Zilhão]

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?

‌

[Miguel Zilhão]

patch for the change that gives problems and sample parfile to run

[Roland Haas]

The function to use is actually GetBoundarySpecification not GetBoundarySizesAndTypes. Easiest is most likely to copy the code from ML_ADMConstraint’s Kranc.cc and ML_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 s cGH* cctkGH turning into cGH ** cctkGH.

Note: KrancBdy_SelectGroupForBC is just a wrapper around Boundary_SelectGroupForBC with some extra code for the PreSync functionality in Cactus.

[Miguel Zilhão]

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…

[Roland Haas]

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 that map is defined) you would end up selecting boundary conditions more than once which makes Boundary abort.

Something similar to what is done when going from GetBoundaryWidths to GetBoundaryWidth 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.

‌

[Roland Haas]

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.

[Zach Etienne]

@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.

[Roland Haas]

When you ran Cactus did you pass it options to show the warnings? Ie -L3 -W3?

[Roland Haas]

I filed a Kranc ticket about this in https://github.com/ianhinder/Kranc/issues/150

[Miguel Zilhão]

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…

[Roland Haas]

@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;

‌

[Roland Haas]

@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. However is_ghostbnd is set as !cctkGH->cctk_bbox[d] in line 425. In LEVEL mode Carpet will set cctk_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 viat CarpetLib::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 in LOCAL mode. My guess would be that Erik misspoke and meant to suggest GetBoundarySpecification for all cases (though note my issues with MultiPatch_GetBoundarySpecification having to be called in SINGLEMAP mode which is below LEVEL mode ie would cause multiple calls to Boundary_SelectGroupForBCs).

[Zach Etienne]

@Roland Haas :Thanks for the explanation!

It is still evident that the backported version, despite depending on deadbeef within an if() 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 use GetBoundarySpecification, some multipatch version, or write something from scratch?

[Roland Haas]

Since Baikal does not support Llama (it does not, does it?) all you should have to do is use GetBoundarySpecification instead of GetBoundarySizesAndTypes. 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!");

‌

[Miguel Zilhão]

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 functionMultiPatch_GetBoundarySpecification, which i guess could be used when Llama is used?

[Roland Haas]

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 in LEVEL 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.

[Miguel Zilhão]

i see. what would be the best way of handling this issue, in the meantime? what happens if i just callGetBoundarySpecification even when Llama is used? or is it preferable to leave things as is until we have input from the Llama developers?

[Roland Haas]

Calling GetBoundarySpecification in Llama code will give you incorrect answers / will fail. The function is provided by the thorn CoordBase which may not be active at all for a Llama run. if it is active the it will report its ow parameter settings and not Coordinates' (the Llama equivalent thorn to CoordBase).

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 as SINGLEMAP then call MultiPatch_GetBoundarySpecifiction and select boundary conditions on only those faces that is_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 the Boundary thorn’s ApplyBCs 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).

[Miguel Zilhão]

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.

[Roland Haas]

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.

[Roland Haas]

Pull request: https://bitbucket.org/canuda/lean_public/pull-requests/3 would take care of this I believe.

[Roland Haas]

Please review (code is by Miguel and Helvi).

[Roland Haas]

@helvi witek , @Miguel Zilhão , @Zach Etienne will review.

[Roland Haas]

@helvi witek @Miguel Zilhão

@Miguel Zilhão

@Zach Etienne any progress on this?

‌

[Miguel Zilhão]

The pull request https://bitbucket.org/canuda/lean_public/pull-requests/3 should fix this. Would this PR be ok to merge?

[Roland Haas]

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.

[Miguel Zilhão]

I’ve added a comment about it in the PR above.

[Roland Haas]

@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.

[Miguel Zilhão]

I’ve merged PR https://bitbucket.org/canuda/lean_public/pull-requests/3 which takes care of this.

[Roland Haas]

Thank you.