Support for reflective boundary conditions?

[antoniusf]

Right now, the tricubic interpolation supports only periodic boundaries. If reflective boundaries are configured, it uses a weird mix where the index of the grid cell is determined using reflective boundaries, but the retrieving the neighbouring points still uses periodic boundaries.

Do we want to fix this before merging?

[LeanderSchlegel]

Thank you :) should now be fixed with a reflectiveBoundary function analog to the periodic one, that uses the same algorithm for reflective boundaries as the closestValue function.

[antoniusf]

Awesome, thanks for fixing! While you're at it, could you also fix the reflectiveClamp (x > n) condition? I think this should be (x >= n), just like the fix we did for reflectiveBoundary, because otherwise the case x == n wouldn't be treated correctly.

[LeanderSchlegel]

Oh thanks! I wanted to do that, but I noticed that there seems to be a problem related to the first calculation that is done in the interpolation functions: Vector3d r = (position - gridOrigin) / spacing; For a position of Vector3d(0,0,0) and an origin of Vector3d(0,0,0) this evaluates to components of -0.5 and the fix would map this forever to -0.5 due to the -1. Probably here two conventions, on the side of the clamping and the side of the coordinate system, collide?

[antoniusf]

Hmm, I don't think I understand yet. Which values for x and n would you have to feed into reflectiveClamp, and what would happen if you did that?

[LeanderSchlegel]

If I see it correctly, for every n the while loop in the fixed reflectiveClamp gets infinite in the case of x=-0.5 due to the -1 that is needed in the case of x==n.

[antoniusf]

Hmm, so this is the code with my suggestion applied (right?):

while ((x < 0) or (x >= n))
  x = 2 * n * (x > n) -x;

Let's take n = 2 and x = -0.5.

First iteration: (x < 0) is true, so the loop body is executed.

x_new = 2 * n * (x > n) - x = 2 * n * 0 - x = -x = 0.5

Second iteration: x is 0.5. (x < 0) <=> (0.5 < 0) <=> false (x >= n) <=> (0.5 >= 2) <=> false

both loop conditions are false, so the loop terminates here.

(this works analogously for all positive integer values of n.)

Have I missed something?

[LeanderSchlegel]

Oh perhaps I am wrong, I used for the fix while ((x < 0) or (x >= n)) x = 2 n (x >= n) -x-1; as in closestValue to make sure the indices are zero based.

[antoniusf]

Oh, I see – hold on, why the minus one?

[antoniusf]

Oh sorry, I see what you mean now. Huh, give me some more time to think this through?

[LeanderSchlegel]

Of course thank you for your help! I also have to think about it for some time, I am not sure right now if the -1 is only due to zero based indices or due to the selection of the mirror axis of the reflection, for example for n=100 i thought of 0 ... 98, 99, 99, 98, ... and the -1 should map the 100 correctly to 99.

[antoniusf]

Soo, I though I'd analyze the current code to figure out exactly which kind of boundary conditions we're going for, and I do think I have an answer, but I also think the current code is kind of wrong.

For reference, I'm working off the most current commit as of this writing; the code we're looking at starts here: https://github.com/CRPropa/CRPropa3/blob/89dacac530b59e009401865847dad2fa75df324e/include/crpropa/Grid.h#L259.

First, in line 261, we perform a coordinate system transformation. The coordinate system of r has the grid point with index (0, 0, 0) at r = (0, 0, 0); and the grid point with index (n-1, n-1, n-1) at r = (n-1, n-1, n-1). (The spacing in r-space is trivially one, since we divide the normal grid spacing out. Also, r = (0, 0, 0) maps to ix = iy = iz = 0 and fX = fY = fZ = 1 independent of boundary condition, which means that b = get(ix, iy, iz) = get(0, 0, 0) – all other terms are zero.)

All further operations are performed in this coordinate system; the original position is not used anymore.

We then determine the indices of the lower and upper grid cells in each axis using reflectiveClamp. I'm going to ignore the case of x=n (and all cases that reduce to this), since this ends up returning lo = floor(x) = floor(n) = n, which is not a valid grid index. The rest of the code works by mirroring values < 0 around 0, and values > n around n, repeatedly if necessary. This creates a (partially reflected) tiling of the origin cube from (0, 0, 0) to (n, n, n).

Recall, though, that the coordinates passed into that function are in r-space, where the lowest grid point is at (0, 0, 0), and the highest is at (n-1, n-1, n-1), meaning that the actual grid is not centered within this tiling. Note that while this mirrored/tiled x isn't directly returned, it is used to compute lo = floor(x), meaning that this weird behavior is represented in the final output.

Interestingly, not even the fractional values fx, fy and fz are computed from this tiled x, as might be expected; rather, they are derived from the unmodified r. This kind of makes sense, considering that the tiling mainly affects the grid points chosen as neighbours and should keep the position relative to the surrounding grid cell unmodified; however, since some of the tiled repetitions of the origin cube are mirrored, but these fractional values aren't, this has the fun consequence of making the interpolation appear to go "backwards" relative to the mirrored parts of the grid.

To visualize these effects, consider the following 3x3(x3) grid, with these values at each grid point:

0.5 0   0
2   4   0
1   3   0

The grid was created with offset 0, and spacing 1. In the image below, the grid was sampled starting at (0.5, 0.5, 0.5) for the lower left corner, going up to (9.5, 9.5, 0.5) in the upper right. White lines highlight the boundaries of the tiling cube, while the black-and-white dots mark the actual, original grid.

Thus:

• note how the mirrored repetitions reproduce the rough structure of the original grid (you can see how they are indeed mirrored to the other side of their tiling cube, due to the aforementioned offset of the grid), while the gradients inside each grid cell follow the same direction as the unmirrored ones instead of being mirrored.
• note also the bright blips that align with the grid lines, where an interpolation query hit a cell boundary exactly and thus, due to floor(x) = x, reproduced somewhat more accurate values
• in addition to this, the points of the top and right edges of the original grid effectively occur three times in series at the tiling boundaries; in contrast, the points of the bottom and left edges are conflated together and thus only occur once. If instead the grid was centered within its tiling cube, the points of an edge would always occur twice at the boundary, regardless of orientation.
• unfortunately hidden by the white lines: memory overflow artifacts from the erroneous x=n condition in reflectiveClamp. I'll post the full image once I figure out how, without cluttering this comment.

The image was created using commit 1aad73b4. This is a bit old because I forgot to check before making the image, and I'm not going to put in all this work again. However, no changes to Grid.h were made since that version, so the file I analyzed in the text above is still the exact same as was used in making this image.

I'm not super sure how to proceed from here – I think the smartest way would be to center the grid inside the tiling cube, such that each border point occurs twice at the boundary. Another option might be to collate neighbouring border points, as is currently done for the bottom and left edges. But this is a physics question, and I don't think I have the insight necessary to answer it. On the technical side, computing the fractional values directly from the reduced/mirrored values should solve the concrete bug as well as simplifying the code; it might also be advisable to add high-level tests that specifically test the symmetry invariants provided by each type of boundary condition (initializing a grid, retrieving a set of random points from the origin cube, and then checking that equivalent points match).

[LeanderSchlegel]

Wow, that's amazing! Thank you for your detailed analysis! I worked with our old jupyter notebook and reproduced your plot :) Then, I modified the reflectiveClamp in order to also give back the mirrored/tiled x directly to calculate therefrom the fractional values to the lower neighbours, it fixes the bug of the interpolation going backwards :) The difference to the upper neighbour however, was still calculated with 1-fX0, would that be correct? To solve the problem of the shifted coordinate systems, I currently try to mirror around 0 and (n-1) with while ((x < 0) or (x > (n-1))) { x = 2 * (n-1) * (x >(n-1)) -x;}

[LeanderSchlegel]

Today @JulienDoerner and I worked on the code and found that the fix I proposed using (n-1) is probably not a good idea, because it would make the volume that is repeated by one spacing smaller, laying the mirroraxis onto a gridpoint, and instead your plot already was the right solution. While the repetition of a gridpoint, like ..98,99,99,98,..., which also Ilja Jaroschewski, Marcel Schroller and I discussed, could be artificial because particles traveling through the boundary would see it exactly twice (the influence should decrease with the gridsize), we think it would make sense based on the definition of the volume that is represented by the grid in CRPropa that extends half a spacing more than the last gridpoint into each direction. We now replaced the (n-1) again with n and fixed the x==n bug with an if statement, mapping the n to n-1.

[antoniusf]

First off, thanks for your work on the code! I'm glad you could reproduce the results, too.

Could you show the solution that you currently have? I agree that having the repetition contain only an (N-1)-size cube might be a problematic somehow, which is why I'd also suggested that we might keep the N-size tiling cube but center the grid within. (In other words, mirror around –0.5 and (N – 0.5) instead of 0 and N.) I can't really tell if this is where you're at right now though.

I also agree that repetition of grid points is probably okay; however, with the old tiling (where the grid isn't centered), the upper right set of grid points is repeated twice (each grid point appears three times), while the lower right set isn't repeated at all (each grid point appears only once). This was my problem with the old tiling solution.

[LeanderSchlegel]

Thank you for your explanation! I think I see what you mean and would also think a centered cube would probably be best, I will try to mirror around -0.5 and (N-0.5) as you proposed :) That's how it looks right now with the latest commit for the reflective boundary and for the periodic boundary. I noticed however, that nearest neighbour is still buggy, I will take a look on that.

[antoniusf]

Thanks for posting your current state, and for taking a look at shifting the mirroring axes.

I'm not sure if this helps with your nearest neighbour problem, but eventually I'd like to do a bit of consistency checking to make sure that all of the function that interact with the boundary (Clamp, Boundary, and closestValue) are mathematically sound and produce compatible results. Yesterday evening I was thinking about how to avoid weird corner cases for integer values of x in reflectiveClamp, because I was worried that there might be discrepancies in the computation of res vs the computation of hi and lo (because the mirroring condition is reproduced there, right?) So I thought about whether we could prove this correct somehow and I came up with this set of conditions that reflectiveClamp should fulfill:

reflectiveClamp(x, N).lo = reflectiveClamp(floor(x), N).res reflectiveClamp(x, N).hi = reflectiveClamp(floor(x) + 1, N).res

Do these make sense? Basically, without boundary conditions we'd just compute lo = floor(x) and hi = floor(x) + 1, so with boundary conditions we kind of want to ensure that the transformed lo, hi, and res follow the same rules.

I started trying to prove the current implementation correct but then realized that it wasn't and stopped. I would like to do this to the implementation eventually to make sure it's correct, but I can also try to work on this a bit and see if I can get a correct implementation this way. (Though for nearestNeighbour, you'd probably have to look at closestValue right? I can definitely recommend the formalization approach though, everything suddenly gets a lot complicated. You just try to work through your proof and if something fails you go to fix it, rather than trying to keep this whole complicated model in your head at once.)

[antoniusf]

Oh, you know what, you were right with the fact that this suggestion of mine:

Awesome, thanks for fixing! While you're at it, could you also fix the reflectiveClamp (x > n) condition? I think this should be (x >= n), just like the fix we did for reflectiveBoundary, because otherwise the case x == n wouldn't be treated correctly.

would lead to an infinite loop, and that it'd do so for exactly the case that I wanted to fix: x == n. (This is a fixpoint in the loop iteration, since x_new = 2 * n * (x > n) - x = 2 * n - n = n.) I assume you knew this already, which is why you didn't use the fix, but I just realized and wanted to share the insight that I'd been wrong.

[JulienDoerner]

Hey,

today @LeanderSchlegel and I discussed the repition of values and we found a problem which may relate to the problem @antoniusf mentioned before

the upper right set of grid points is repeated twice (each grid point appears three times), while the lower right set isn't repeated at all (each grid point appears only once)

If we asume a 1D grid with values 1, 0, 2 we would expect to repeat them as 1, 0, 2, 2, 0, 1, 1, 0, 2 ... Till now the repetition in the plots from Leander have to many 2 and the 1 is not repeated. Therefore the big gap occures.

The problem comes from the definition of the upper and lower neighbour in the reflectiveClamp function. Wouldn't it be easier to use a function which clamps the local r vector to a position in the regular grid and then calculate the neighbours with floor(x) and floor(x)+(x<n).

[LeanderSchlegel]

Hey, thank you for your ideas! I thought about how to fix the problem of the fixpoints of the loop (@antoniusf Thank you for your reply, that got me onto this :)) and tried to avoid the loop by determining the boxnumber the point is in, via a division and afterwards determining the value inside the box where the point should get mapped to, via modulo, after seeing that there are only two different cases that have to be considered depending on whether the boxnumber is even or odd, because for the mirroraxis we chose, in one case the box indices look (for n=3) like 0 1 2 and in the other 2 1 0. -5 -4 -3 -2 -1 0 1 2 3 4 5 x / Indices 1 2 2 1 0 0 1 2 2 1 0 f(x) / Mapped Indices I discussed it with @reichherzerp today, and we found the problem that even without using a loop the problem of a fixpoint seems unavoidable in the case of a mirroraxis that lies between two gridpoints, because the axis itself would have to be mapped onto itself, like -0.5 -> -0.5 and also every point in the interval (-1,0) (for an axis of -0.5) would infinitely get mirrored back and forth around -0.5. Do we see it correctly?

We then discussed if instead the idea of a mirroraxis that lies on a gridpoint (like in the second last commit) could nevertheless be a good idea, because next to the avoided problem of an artificial repetition, the not reproduced volume of the cube, which @JulienDoerner argued, (also in comparison to the periodic boundary, where the volume has to be reproduced), would make just a very little influence with increasing gridsize, going with 1/(n^3).

[reichherzerp]

An additional argument is that having the mirroraxis that lies between gridpoints results in the artificial repetition of values. This artificially amplifies local minima and maxima at borders. This is avoided with the method of mirroring on grid points.

[antoniusf]

Hi, so I'm not quite sure what the problem with the fixpoint is? When mirroring, the mirror axis is always a fixpoint, no matter where it lies. I'd just been wrong about how to handle this. I now think the correct way to do it is to include both mirror axes in the final (reduced) volume. For mirroring around -0.5 and (N – 0.5) the final volume is thus [-0.5, N – 0.5]³. If the loop breaks once a point lies within this volume, there shouldn't be an infinite loop? Let's take -0.8 \in (-1, 0) as an example: The first iteration will map this to -1 – (-0.8) = -0.2. This is in [-0.5, N – 0.5], so the loop will terminate; the same goes for all other values within (-1, -0.5). Values in [-0.5, 0) are already in the reduced volume and thus don't need to be reduced further.

However, @reichherzerp's note still stands. As I said above,

But this is a physics question, and I don't think I have the insight necessary to answer it.

I agree though that the repeated points do feel a bit icky.

[JulienDoerner]

Hey, the problem I see with the reduced volume is not for the usage with magnetic fields on grids but ruther more for gas densities and photon fields. In both cases a repition of values should not be a problem, but the volume representation is important at the edge. For this case we would also need a way where the grid is repeated with the last gridpoint outside the original volume. An example would be a galactic gas density which drops down at edge of the Galaxy to the density of the IGM and would be (nearly) constant from that on. For magnetic fields i would just recomend to use the periodic repetition of the grid to avoid these artifitial fields.

[LeanderSchlegel]

Hi, @antoniusf: Sorry, I was wrong, you were completely right that the mirror axis always has to be a fixpoint and that no infinite loop is produced. I tried to implement your solution, and my code reads:

` while ((x < -0.5) or (x > (n-0.5)))

    x = 2 * n * (x > (n-0.5)) -x-1;
res = x;
lo = floor(x);
hi = lo + (lo < n-1);
`

For the trilinear interpolation and sampling starting at (0.5,0.5,0.5), it looks like this: however, it looks weird starting at (0,0,0), perhaps due to negative values that are given back when x is in [-0.5,0), which i tried to handle with if (x<0) {lo=0;hi=0;res=0.0;} , however, that statement changed the picture making it more symmetric than it should be regarding the gradients given by the grid values.

For the tricubic interpolation I noticed artefacts that I am not sure about, and the picture seems also to be symmetrized. The artefacts reach negative values in some regions.

[LeanderSchlegel]

I tested that negative x values are rounded down by floor(x) to lo=-1 which is given to the get() function and think a handling for these values is needed. For all values in [-0.5,0) I think the neighbours should both be 0, however, res should probably not be changed to maintain the real distances to these points, perhaps: if (x<0) {lo=0;hi=0;} after the while-loop.

[antoniusf]

I agree that res shouldn't be changed. How about clamping lo to be ≥ 0, i.e. lo = max(floor(x), 0) or something? That way it's nice and symmetric with hi, which is clamped to ≤ n – 1.

[antoniusf]

@JulienDoerner I take what you're saying as an argument that we should keep the repeated points for reflective boundary conditions, is that right? Also, what's the second thing you meant to say? That the last gridpoint in each direction just gets repeated out to infinity?

[antoniusf]

@LeanderSchlegel Also, thanks so much for all the work you've been doing. I've been stuck on trying to find a nice correctness proof while I could have been putting a lot more effort into an actual implementation.

[LeanderSchlegel]

@antoniusf Oh thanks, that's so nice :) I also have to thank you and everybody for your help!

[JulienDoerner]

@antoniusf Yes that are cases where the repetition of gridpoints is neccesarry. And for the second point I mentioned is the restriction of a grid to a given volume an set all other values to zero (or the last value given at the corner, I don't know which case is better). I think this could also be done with the RestrictToRegion Module (line). But this is only for modules not for magnetic fields.

[LeanderSchlegel]

@JulienDoerner and I discussed the case of negative values we saw when using the tricubic interpolation and he had the idea, that even when all grid values are positive, the polynom could be bend towards the negative halfspace in case of two following zeros but non-zero slopes in these points. I tested it in 1D by plotting the cubic spline (see https://www.paulinternet.nl/?page=bicubic) for the values 1,0,0,1.

[LeanderSchlegel]

There are also neighbouring cases like 1,0.1,0,1 for which the curve gets negative.

[LeanderSchlegel]

We discussed that this is therefore in principle expectable (I'm not sure if for all cubic polynoms, at least for our polynom) and I would think it is okay to leave it like that at this point, however, for scalar fields like photon densities, negative values could make problems and the trilinear interpolation which would avoid that, could be the better choice.

[LeanderSchlegel]

I thought about using lo=max(floor(x),0) and stayed now with if (x<0) {lo=0;hi=0;} after the while-loop, because otherwise hi would get 1. I hope this way all cases that are left after the clamping loop do get covered correctly.

Here is how it looks for the trilinear interpolation, does the picture look correct to you? I was not sure about the symmetric looking inner part of the four cubes, however, before it probably was not correct with the dark blue/zero lines, because the zeros should be on the outside. Interestingly this fix does not produce comparable results for the tricubic interpolation, where it instead seems to look comparable without the fix. I tested it by effectively leave the fix out/ set the lower neighbours explicitly to -1 for res in [-0.5,0), perhaps this is related to the definition of the used polynom and its sampling points of -1,0,1,2 see https://www.paulinternet.nl/?page=bicubic) and the picture looks like this: This is how it looks for nearest neighbour interpolation: I just pushed this current state to our branch.

[antoniusf]

Yay, yes, I think the trilinear picture looks good now. The edges and center of each of the squares have these small bands of constant values, which are what we would expect when two successive points are at the same value right?

As for the tricubic interpolation, I think that your explanation, that the negative values come from undershooting of the cubic polynomial, is sufficient. And the fix for lo and hi only applies to trilinear interpolation anyways, right? Tricubic just uses reflectiveBoundary to determine all of its gridpoints and thus doesn't need that logic.

Ohhhh wait hold on – tricubic currently uses reflectiveClamp? I'm a bit confused, I don't think this should be the case? I've pushed a change to enhanced_interpolation_afrie, could you either test that or send me the jupyter notebook that you've been using to make these plots?

Reasons that I think this is wrong: Basically, we always have to choose whether we do operations in the original coordinate system (r-space as I've called it above), or in the reduced coordinate system. reflectiveClamp computes lo and hi in the reduced coordinate system, which is why we have to use res from the reduced coordinate system as well. The reason is that since reflective boundary conditions mirror the coordinate space, lo and hi can swap places. Eg: if floor(x) ≤ x ≤ floor(x) + 1 in r-space, then it is very well possible for reduce(floor(x) + 1) ≤ reduce(x) ≤ reduce(floor(x)) after the coordinates have been reduced into the origin cube.

However, for the tricubic interpolation, computing the neighbours in the reduced space becomes unnecessarily complicated and error prone, which is why we introduced the reflectiveGet function. Using this function, we compute the neighbours in the (non-reduced) r-space, just by adding and subtracting, and then perform the reduction using reflectiveGet/reflectiveBoundary. This not only means that we don't need reflectiveClamp (since we'll perform the reduction later anyways), but it also means that using the res returned from reflectiveClamp to perform the computation of (fX, fY, fZ) should result in incorrect values... (though I don't quite understand why the tricubic plot you've made wouldn't show more apparent artifacts if my line of thinking was correct, but the reason might be that we're working off of the lo returned by reflectiveClamp. In either case, I think it's worth trying the other approach and seeing what changes.)

[LeanderSchlegel]

Thank you very much for your improved implementation!

I am now performing tests using real turbulent fields and thought about a high-level test scenario using random points and checking if the boundary conditions are fulfilled, however, I think then one would have to implement the same conditions in the test script and code would get checked against code. I came up with another idea, in which a grid is taken and sampled with an offset in the sampler of n*spacing/2. in the x and y direction, therefore four quadrants of reflection are covered. To cover the symmetry in the z-axis as well, @reichherzerp suggested to check the same for an offset in x and z-direction.

I sampled on the gridpoints itself and it looks good! XY: XZ:

The doubled lines at the mirror axis are visible and I checked the whole volume numerically by comparing the opposite lines (middle and middle, middle-1 and middle+1,...) in one plane for all planes for left to right and up to down, comparing each vector component with each other, a single failure would produce an exception and the test would fail. For TL and NN they are exactly equal, for TC at least within an absolute difference of 2e-6.

[LeanderSchlegel]

We now tested also sampling between gridpoints for all interpolations and it works as expected :)