AEILocalInterp should not off-centre the interpolation stencil by default

[Ian Hinder]

AEILocalInterp by default will off-centre the interpolation stencil if there are insufficient points to perform the interpolation. In a parallel setting, this could happen due to there being insufficient ghost points for the interpolator chosen. The off-centering leads to an interpolation error which is of the correct order but larger than for a centered stencil. More importantly, it leads to different results on different numbers of processes. This violates a basic design principle of Cactus, and the expectation of users, that changing the number of processors should not change the results of a simulation.

I propose that instead of silently off-centering the stencil, AEILocalInterp should abort with an error indicating that there are insufficient ghost-zones. The interpolator options corresponding to this are:

boundary_off_centering_tolerance={0.0 0.0 0.0 0.0 0.0 0.0} boundary_extrapolation_tolerance={0.0 0.0 0.0 0.0 0.0 0.0}

Keyword: AEILocalInterp

[Erik Schnetter]

Thornburg suggests (in the code) to use 1.0e-10 instead of 0.0 as the default values to avoid round-off problems. Comments?

[Roland Haas]

I would vote for making AEILocalInterp abort. The respective changes would be in src/InterpLocalUniform.h where LAGRANGE_BNDRY_OFF_CNTR_TOL_DEF etc are defined. The documentation always said that sufficient ghost points must be present for the interpolator stencil, so no well behaved user code should see a difference.

[Erik Schnetter]

Agreed.

This should wait until after the release.

[Ian Hinder]

Increasing priority to major, since this can lead to reproducibility issues when the number of processors is changed, and can cause people to spend a lot of time looking for the problem. It would be good to implement this before the next release.

[Ian Hinder]

We just ran into this issue again. However, I don't understand in our case why there is a problem. We are using order=4 Lagrange interpolation, which has a 5-point interpolation stencil, and we use 3 ghost and symmetry points. Hence, the interpolator should never have to off-centre the stencil, because the interpolation point will always be in the 'default-centering-region', as defined in the documentation (https://einsteintoolkit.org/documentation/ThornDoc/Numerical/AEILocalInterp/). In fact, 2 ghost points should be sufficient. We are using periodic boundary conditions, and the interpolation point is exactly on the boundary in one (maybe two) dimensions. If CarpetInterp always gives the point to the component that "owns" it, and the symmetry interpolator the same, then setting these options (boundary_off_centering_tolerance= boundary_extrapolation_tolerance={0...}) should have no effect. Instead, we see that without these parameters set, there is a symmetry violation (a geodesic which should evolve along the cell edge by symmetry diverges from it), and with them set, everything is fine. Further, my understanding was that if AEILocalInterp previously decided to off-centre, and you set these options to 0, then it would return an error code, which we do not see, rather than giving "better" results. If it could give "good" results without off-centering, then why is it deciding to off-center?

[Erik Schnetter]

With periodic boundaries, which grid points are boundary points and which ones are evolved? Where are the grid points placed with respect to the boundary? There are basically three options: (1) an evolved grid point on the lower boundary, (2) an evolved grid point on the upper boundary, and (3) grid points staggered about the boundary. The two former go with vertex centred refinement, the latter with cell centred refinement.

Where is the boundary located according to AEILocalInterp? If this thorn assumes that the boundary is directly at the outermost evolved grid point, then this notion clashes with what CoordBase chooses (see above), as boundary according to CoordBase will necessarily be located further outwards. Thus you might need to all off-centring the stencil, simply because AEILocalInterp's notion of what off-centring means is different from CoordBase's notion.

Thorn CoordBase has some nice sketches explaining things that I tried to explain above.

[Eloisa Bentivegna]

The boundary Ian refers to is of type (1): its points are evolved. I am not sure how AEILocalInterp sees this, though.

[Ian Hinder]

Erik: This is unigrid, so I think there is no difference between vertex and cell centring. Does your question relate to whether the upper or lower boundary is open or closed? Usually, the boundaries are specified as [0, L), so we have an evolved point at 0, and three symmetry points < 0, and the last evolved point would be L - h, followed by three symmetry points. The lower boundary would have a shiftout of 1, and the upper 0. As far as I can tell, AEILocalInterp doesn't look at anything other than the options you pass in; there is provision for passing a number of boundary points, but we don't use it, and it defaults to 0. So AEILocalInterp thinks that all points can be interpolated to and from.

I think that this is fine, because the symmetry thorn and CarpetInterp should only ever be asking for points for which there are enough ghost zones, in our case.

PeriodicCarpet folds the interpolation coordinates using the physical_min and physical_max, which it gets from CoordBase's GetDomainSpecification. Is this logic implemented correctly? Are these the first and last evolved points? If PeriodicCarpet is leaving the coordinate "unfolded", so that it is actually outside the domain for which we have enough points, then this might explain some of what is happening.

But still, even in that case, if the interpolator was previously having to off-centre, then setting the options to tell it not to should result in an error, not in the "right" answer. So I am still confused.

[Erik Schnetter]

Let's assume you have 10 points on your grid, and your domain is [0; 1). Your points' coordinates are thus [0.0, ..., 0.9]. If you try to interpolate at 0.95, then this is fine, as it's inside the domain. However, all AEILocalInterp sees is that it has grid points in the range [0.0; 0.9], and it will thus complain.

[Frank Löffler]

Replying to [comment:9 eschnett]:

Let's assume you have 10 points on your grid, and your domain is [0; 1). Your points' coordinates are thus [0.0, ..., 0.9]. If you try to interpolate at 0.95, then this is fine, as it's inside the domain. However, all AEILocalInterp sees is that it has grid points in the range [0.0; 0.9], and it will thus complain.

Why should it complain? With 3 symmetry points, the actual grid should go from -0.3 ... 1.2 (1.3?), and with an interpolation stencil of 2 (0.2) in each direction, shouldn't 0.95 be fine? Why should AEILocalInterp care about the size of the evolved grid, and not only on the size of the grid it has to interpolate on?

[Erik Schnetter]

Can you remind me what the actual case is where you receive the complaint?

[Ian Hinder]

This is the problem! There is no complaint. When we don't tell it not to off-centre (the default), there is some sort of symmetry violation, suggesting that it is off-centering. When we tell it not to off-centre, this symmetry violation goes away, so it looks like it is no longer off-centering. This doesn't make any sense. Telling it not to off-centre should cause it to give an error if it had been off-centering before. The case we have is the one in comment:5.

Let's assume you have 10 points on your grid, and your domain is [0; 1). Your points' coordinates are thus [0.0, ..., 0.9]. If you try to interpolate at 0.95, then this is fine, as it's inside the domain. However, all AEILocalInterp sees is that it has grid points in the range [0.0; 0.9], and it will thus complain.

Are these just the evolved points? We would also have three symmetry points at [-0.3,-0.2,-0.1] and [1.0,1.1,1.2]. Suppose we ask for interpolation at 0.95. This is outside the 'logical' domain, so PeriodicCarpet should fold it back. Aha! Now I see something odd. This point is equivalent to 0.95 - 1, which is -0.05, which is STILL outside the logical domain. So there is actually no interpolation point we can use which would lie inside [0,1). Is this to do with using vertex centring? What would PeriodicCarpet do in this case? Even so, I still don't understand the problem. Suppose PeriodicCarpet does fold it back to -0.05. AEILocalInterp will see the points [-0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, ...]. It needs 5 points to do the interpolation, so it should use -0.3, -0.2, -0.1, 0.0, 0.1 (or shifted to the right by one; I find it hard to parse the open/closed diagran in the AEILocalInterp docs.) So it should still work, and should never have to off-centre. We get away with this because we have three ghost zones, and this interpolator only needs two (modulo the PeriodicCarpet issue above).

[Roland Haas]

Ian, since the point in question is just at the edge of the grid (and a cell) is it possible that you are being hit by some sort of roundoff issue?