Multi-block boundaries leave uninitialized boundary points

[Erik Schnetter]

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:

[Ian Hinder]

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?

[Erik Schnetter]

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.

[Steven R. Brandt]

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?

[Ian Hinder]

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.

[Erik Schnetter]

This https://bitbucket.org/einsteintoolkit/mclachlan/pull-requests/4/correct-scalar-initial-boundary-conditions/diff addresses the second bug. (Not the one that Ian addresses above.)

[Erik Schnetter]

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.

[Erik Schnetter]

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.

[Erik Schnetter]

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

[Steven R. Brandt]

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?

[Erik Schnetter]

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.

[Roland Haas]

Are these bugs still alive or have both pull requests:

• https://bitbucket.org/einsteintoolkit/mclachlan/branch/ianhinder/initgammas
• https://bitbucket.org/einsteintoolkit/mclachlan/pull-requests/4/correct-scalar-initial-boundary-conditions/diff

been applied? If these are ready to be applied, please change the ticket state to "review" to request a review.

[Roland Haas]

Ian?