1. Einstein Toolkit Avatar Einstein Toolkit
  2. Einstein Toolkit
  3. tickets
  4. Issues

Multi-block boundaries leave uninitialized boundary points

Issue #2030 resolved
Erik Schnetter created an issue

Applying multi-block boundaries can leave uninitialized boundary points. This happens when multi-block boundaries are used in combination with another boundary condition that looks at the interior of the grid, such as e.g. a symmetry condition or radiative / extrapolating boundaries.

There is another bug currently in McLachlan that applies its "scalar" boundary conditions to too few grid points (it uses "BoundaryNoSync" instead of "Boundary"), which means there are uninitialized boundary points for initial_boundary_condition = "scalar" as well in multi-block systems.

Keyword:

Comments (24)

  1. Ian Hinder
    • removed comment

    The first bug mentioned is one that we have discussed before, many years ago. One solution was to apply the BCs twice, which would probably need to be handled in Kranc. (I just looked, and cannot find a Kranc issue about this, unfortunately.) This was the reason that for a long time, my McLachlan Llama BBH parameter files used poison = 0, because setting the uninitialised points to zero was better than having them as NaN, especially when the boundary was causally disconnected. However, recently, I was not able to reproduce the problem. I thought that maybe something had been fixed/changed, though this was perhaps optimistic. Do you have a parameter file and a method for identifying that this problem has occurred?

    Presumably the long-term solution is for Cactus to keep track of variables read and written in different parts of the grid. Is there a short-term solution that doesn't involve "clever" scheduling, and applying the BCs twice?

    As I understand it, this problem occurs only on initial data, not during evolution. Is that your understanding as well?

  2. Erik Schnetter reporter
    • removed comment

    Yes, this is the same problem.

    You can avoid it by using both initial_boundary_condition="scalar" and rhs_boundary_condition="scalar", and by not using any symmetries. Any other multi-block case, i.e. either using bitant mode, or using extrapolation for the initial boundary condition or radiative for the rhs boundary condition, should trigger the issue.

    The second bug compounds the issue, because the initial boundary condition is not handled correctly even in the scalar case. I will propose a patch soon.

    As you say, the best current solution for this problem is to apply the boundary conditions twice, both before and after multi-block interpolation. The true solution, which is not supported by Cactus's boundary conditions, is to distinguish between faces, edges, and corners, and to apply some of these before and the others after multi-block interpolation. While this still requires boundary conditions twice in the schedule, they will be executed only once per boundary point.

  3. Steven R. Brandt
    • removed comment

    It's not obvious to me what is different about synching between multiblock grids and regular ghostzone points. Aren't these operations performed at the same time?

  4. Ian Hinder
    • removed comment

    When syncing ghost zones (for interprocess communication), each ghost point can be updated in precisely one way, which is to set its value to the value held by the process which "owns" the real point. i.e., we have only one "boundary condition" which can be applied to that point. A "standard" interpatch boundary point is the same; in the usual case, there is a single component on a single process which owns the point, and the boundary point is filled by interpolation from the neighbouring points in that component. Similarly, a "standard" outer boundary point can be filled by applying the physical boundary condition. However, a problem arises when a point is BOTH an outer boundary point and an interpatch boundary point. This would happen for a point at r = r_max and phi = phi_max.

    I think this problem only happens during initialisation, not evolution. The steps are these:

    1. CCTK_INITIAL/ADMBase_InitialData: TwoPunctures initialises the ADMBase metric variables pointwise at every point on the grid.
    2. CCTK_INITIAL/ML_BSSN::ML_BSSN_InitialADMBase2Interior: McLachlan computes the BSSN variables from the ADMBase variables. The BSSN Gamma variables are initialised using finite differencing, since they are computed from derivatives of the metric. This means that they are set only on "interior" points.
    3. CCTK_INITIAL/ML_BSSN_Helper::ML_BSSN_ExtrapolateGammas: Fill the outer boundary of the Gammas by extrapolation (radial) from the interior.
    4. CCTK_POSTINITIAL/MoL_PostStep/Interpolate2::Interpolate2ApplyBC: The interpatch boundaries are applied, but ONLY to points which have a wide enough stencil. Specifically, the points on both the outer boundary and the interpatch boundary are NOT set.

    The extrapolation (step 3) extrapolates uninitialised points into the "corners" (those points which are both r = r_max and phi = phi_max) because the phi interpatch boundary points haven't been set yet. The interpatch boundary condition doesn't touch these points, because the interpolation stencil cannot be applied so close to the outer boundary. Hence, the points remain uninitialised.

    If the interpatch BC was applied first, its stencil would include the uninitialised outer boundary points, and corrupt the interior. The extrapolation into the outer boundary would then use this corrupted value, so you would have uninitialised data in both the interior and the exterior.

    A working solution is to follow the above steps, which leaves uninitialised corner points, and then apply the radial extrapolation a second time. We could probably implement this by manually scheduling ML_BSSN_Helper::ML_BSSN_ExtrapolateGammas in CCTK_POSTINITIAL after ML_BSSN_ADMBase_ApplyBCs. I think the reason I didn't do this originally was that I wanted a general solution to this problem in Kranc, and that was a harder problem. Here, we seem to be using the helper thorn to do this extrapolation anyway, so writing manual code to fix it seems like the best option. ... I have implemented this and put it on a branch of McLachlan. See https://bitbucket.org/einsteintoolkit/mclachlan/branch/ianhinder/initgammas. It compiles and passes the tests, but there is no test for this problem. If someone has an example which currently shows uninitialised points, please try it with this version and see if it is fixed.

  7. Erik Schnetter reporter
    • removed comment

    Ian, you analysis is correct. I have not yet tested your patch. However, the problem also exists during evolution. There we use radiative boundaries, which are similar to extrapolation in their requirements. I will attach a test case for your patch.

  8. Erik Schnetter reporter
    • removed comment

    The test case above activates the NaNChecker. It will abort if the error is present, i.e. if there is a nan. If should succeed if there is no nan.

  10. Steven R. Brandt
    • removed comment

    I'm still a little confused. So you have a portion of the grid which is both a physical boundary and a synchronization zone... that should lead to zones being filled in multiple times, not being filled in at all, shouldn't it?

  11. Erik Schnetter reporter
    • removed comment

    Synchronization does not play a role here. (Let's be exact about terminology to avoid confusion.) Yes, there is part of the grid -- an edge -- that is both an inter-patch boundary and a symmetry boundary.

    Let's also assume that we look at the lower left corner of a 2D grid for simplicity. Let's assume that the x<0 boundary is a multi-block interpolation boundary, and that the y<0 boundary is a symmetry boundary. Initially, the interior of the grid (x>0, y>0) is defined.

    The multi-block interpolation requires an interpolation stencil. It can never fill the y<0 part of the grid, since the stencil doesn't fit. It can only fill the (x<0, y>0) boundary, and for this, it requires all (x>0) points to be defined so that the stencil can be evaluated. That means it will look at (y<0) points, so it needs to be applied after the symmetry boundary.

    The symmetry boundary can define all (y<0) points, and requires the respective (y>0) points to be defined for the same x coordinate. So initially, since all boundaries are undefined, we can only set the (x>0, y<0) part of the boundary via the symmetry condition. Since we don't have any values for (x<0, y>0), we cannot define the lower left corner (x<0, y<0) with a symmetry condition (yet).

    However, after applying the symmetry condition and defining the (x>0, y<0) boundary, now all points with (x>0) are defined. We can now apply the multi-block interpolation, which fills part of the (y<0) boundary. Since the stencil has a finite size, this defines the (x<0, y>0) points, leaving the lower left corner still undefined.

    And now we can, in the final step, apply the symmetry condition again, this time to define the missing (x<0, y<0) corner.

  15. Roland Haas

    @Ian Hinder , @Erik Schnetter any updates on this? The last “real” activity is from 2017-10. After that is is just me pinging you.

  17. Roland Haas

    @Ian Hinder @Erik Schnetter any updates on this? The last “real” activity is from 2017-10. After that is is just me pinging you.

  18. Ian Hinder

    @Roland Haas I don’t have any updates on this. I’m not working actively with Llama, and very likely won’t have time to follow this up in the near future.

  20. Roland Haas

    Since this has been bugging me for a while and I periodically get hit by it myself, reassigning to myself.

    I did a bit of digging into this and it turns out that the “usual” 7patch system used does indeed do fill in all combined physical boundary / interpatch boundary points in Interpolate2’s ApplyBCs. However this happens only if the domain in Coordinates::symmetry is set to full ie no symmetry is used.

    This is ultimately due to a check in Interpolate/init.cc that looks at the value of the grid scalar Coordinates::interpolate_boundary_points which is set to 1 for the 7patch system only if symmetry is full (see https://bitbucket.org/llamacode/llama/src/839808fc05fdfdc261cc1da24c4038f170057bde/Coordinates/src/thornburg04.cc?at=master#lines-472):

    *interpolate_boundary_points = CCTK_EQUALS(symmetry,"full") ? 1 : 0;

    So I would content that in the example parfile (https://bitbucket.org/einsteintoolkit/tickets/issues/attachments/2030/einsteintoolkit/tickets/1549060822.21/2030/ML_BSSN_MP_extrapolation.par) the issue is not so much scalar vs extrapolate boundary condition but more the use of z-symmetry. Namely the parfile works fine with extrapolate b/c if one disables z-symmetry.

    Having said that, for cases where interpolate_boundary_points==0 the issue certainly exists since Interpolate2 will in that case not fill in interpatch boundary points that are also outer boundary points (physical or symmetry). Ian’s proposed fix in https://bitbucket.org/einsteintoolkit/mclachlan/branch/ianhinder/initgammas should fix that indeed at t=0.

    For t>0 the NewRad boundary condition is encoded in the RHS and no physical is selected for any of the evolved variables (eg ML_Gamma_bound = “none”) but symmetry (Llama’s actual symmetry and Llama’s use of Cactus symmetry for interpatch boundary interpolation) conditions are still applied. Since for interpolate_boundary_points==0 the interpatch interpolation does not happen in outer boundary points the RHS computed must actually be valid in the points that are both outer-boundary and interpatch boundary points ie. not just in the interior which is unusual and not provided by NewRad (certainly not for z_is_radial = true since it wants to do 1d extrapolation in the 3rd direction, and also not for regular NewRad since it will not extrapolate into regions marked by Cactus as symmetry points ie the interpatch boundaries since git hash 080c20b "NewRad: do not extrapolate into symmetry boundaries" of einsteinevolve).

    My opinion would be that Llama really should not require science thorns to compute a valid RHS anywhere other than the interior of the grid and should itself take care that it can fill in the regions that it claims to be symmetry points, if this means that it needs to shift its interpolation stencil near the outer boundary then so be it. I would argue that complexity, if unavoidable, is for the infrastructure code (ie Cactus, Carpet, Llama) to handle and not to be handed down to the science code .

    For the special case of z-symmetry in the 7patch Thornburg04 patch system, branch rhaas/symmetry_ip_interpolate (https://bitbucket.org/llamacode/llama/src/8c7c70abadebadc4541793944501e634b919b1e3/?at=rhaas%2Fsymmetry_ip_interpolate) re-enables interpatch interpolation in (physical) outer boundaries even when symmetry != “full” since the symmetry code already has special handling in there to handle points that are both symmetry and interpatch boundary points. With the branch the example parfile from the top of the ticket runs without creating NaNs.

  22. Roland Haas

    Approved by @Erik Schnetter in pull request.

    Applied as git hash 14f6053 "Coordinates: avoid compiler warning" of llama

    Will keep open due to issue still remaining an this being a partial fix only.

  23. Roland Haas

    The remaining issue does not seem to have been due to unitialized variables. Namely the effect in my BBH test run seems to be an exponential growth of oscillations in the conformal factor in all the physical boundary points (with or without the z-symmetry). At least it does not trigger any of the unitialized data checks in Llama itself that already fills the multipatch which exists in Interpolate2 all the time:

    738:     // interpolation source
739:     if (verbose) CCTK_INFO ("Poisoning inter-patch boundaries");
740:
741:     {
742:       list<scatter_setup_t>::const_iterator scatter_setup_it =
743:         interp2_setup.scatter_setups.begin();
744:       BEGIN_LOCAL_MAP_LOOP (cctkGH, CCTK_GF) {
745:         BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) {
746:           assert (scatter_setup_it != interp2_setup.scatter_setups.end());
747:           scatter_setup_t const & scatter_setup = *scatter_setup_it;
748:           assert (scatter_setup.m == Carpet::map);
749:           assert (scatter_setup.c == Carpet::component);
750:
751:           DECLARE_CCTK_ARGUMENTS_CHECKED(Interpolate2ApplyBC);
752:
753:           for (int n = 0; n < nvars; ++ n) {
754:
755:             CCTK_REAL * restrict const var =
756:               static_cast <CCTK_REAL *>
757:               (CCTK_VarDataPtrI (cctkGH, 0, indices.AT(n)));
758:             assert (var);
759:
760: #pragma omp parallel for
761:             for (int i=0; i<int(scatter_setup.indices.size()); ++i) {
762:               var[scatter_setup.indices.AT(i)] = NAN; // poison;
763:             }

    which ensures that none of the target points' data is accidentally used as data for the interpolation itself. This is not quite the same as poisoning all interpatch boundary points just before boundary conditions get applied but is close already.

  25. Log in to comment
Assignee
Type
bug
Priority
major
Status
resolved
Component
EinsteinToolkit thorn
Milestone
–
Version
development version
Votes
0
Watchers
0
Jira Software: the preferred issue tracker for Bitbucket. Join the team!