 removed comment
Applying multiblock boundaries can leave uninitialized boundary points. This happens when multiblock 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 multiblock systems.
Keyword:
Comments (23)


reporter  removed comment
Yes, this is the same problem.
You can avoid it by using both
initial_boundary_condition="scalar"
andrhs_boundary_condition="scalar"
, and by not using any symmetries. Any other multiblock 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 multiblock 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 multiblock interpolation. While this still requires boundary conditions twice in the schedule, they will be executed only once per boundary point.

 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?

 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:
 CCTK_INITIAL/ADMBase_InitialData: TwoPunctures initialises the ADMBase metric variables pointwise at every point on the grid.
 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.
 CCTK_INITIAL/ML_BSSN_Helper::ML_BSSN_ExtrapolateGammas: Fill the outer boundary of the Gammas by extrapolation (radial) from the interior.
 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.

 changed status to open
 removed comment

reporter  removed comment
This https://bitbucket.org/einsteintoolkit/mclachlan/pullrequests/4/correctscalarinitialboundaryconditions/diff addresses the second bug. (Not the one that Ian addresses above.)

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.

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.

reporter  removed comment
Ian, I tested your patch. Unfortunately, it does not work.

 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?

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 interpatch 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 multiblock 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 multiblock 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 multiblock 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.

 removed comment
Are these bugs still alive or have both pull requests:
 https://bitbucket.org/einsteintoolkit/mclachlan/branch/ianhinder/initgammas
 https://bitbucket.org/einsteintoolkit/mclachlan/pullrequests/4/correctscalarinitialboundaryconditions/diff
been applied? If these are ready to be applied, please change the ticket state to "review" to request a review.

 removed comment
Ian?

@Ian Hinder @Erik Schnetter are these bugs still alive or have both pull requests:
 https://bitbucket.org/einsteintoolkit/mclachlan/branch/ianhinder/initgammas
 https://bitbucket.org/einsteintoolkit/mclachlan/pullrequests/4/correctscalarinitialboundaryconditions/diff
been applied?

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


assigned issue to
 edited description

assigned issue to

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

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


assigned issue to
 marked as major

assigned issue to

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 tofull
ie no symmetry is used.This is ultimately due to a check in
Interpolate/init.cc
that looks at the value of the grid scalarCoordinates::interpolate_boundary_points
which is set to 1 for the 7patch system only ifsymmetry
isfull
(see https://bitbucket.org/llamacode/llama/src/839808fc05fdfdc261cc1da24c4038f170057bde/Coordinates/src/thornburg04.cc?at=master#lines472):*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
vsextrapolate
boundary condition but more the use of zsymmetry. Namely the parfile works fine withextrapolate
b/c if one disables zsymmetry.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
theNewRad
boundary condition is encoded in the RHS and no physical is selected for any of the evolved variables (egML_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 forinterpolate_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 outerboundary and interpatch boundary points ie. not just in the interior which is unusual and not provided by NewRad (certainly not forz_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 zsymmetry in the 7patch Thornburg04 patch system, branch rhaas/symmetry_ip_interpolate (https://bitbucket.org/llamacode/llama/src/8c7c70abadebadc4541793944501e634b919b1e3/?at=rhaas%2Fsymmetry_ip_interpolate) reenables 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. 
Pull request for z symmetry is here: https://bitbucket.org/llamacode/llama/pullrequests/7/rhaassymmetryipinterpolate


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 zsymmetry). 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 interpatch 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.
 Log in to 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 longterm solution is for Cactus to keep track of variables read and written in different parts of the grid. Is there a shortterm 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?