Improved interpolation routines and bugfix of reflective boundary

[LeanderSchlegel]

Dear all,

with this pull request we want to contribute an additional tricubic- and nearest neighbour interpolation for gridbased data to CRPropa, as reviewed regarding the technical properties and physical effects in (Schlegel et al 2020 ApJ 889 123 / https://arxiv.org/abs/1907.09934). Also we want to fix the bug in the reflective boundary condition that was described in issue #361

New interpolations

First, we introduced a new grid property named "interpolationType" allowing the user to switch between the three now available interpolation routines for scalar and vectorgrids, setting it to "TRILINEAR" , "TRICUBIC", "NEAREST_NEIGHBOUR" ,respectively, while trilinear interpolation remains the default case.

The new nearest neighbour interpolation could simply be implemented by using the existing closestValue function.

The new tricubic interpolation is implemented according to https://www.paulinternet.nl/?page=bicubic and http://graphics.cs.cmu.edu/nsp/course/15-462/Fall04/assts/catmullRom.pdf.

Optimization

The tricubic interpolation was SIMD optimized on vectorgrids to gain faster execution times. However, there is no option to use it without the SIMD extension, therefore, we introduced a new CMAKE flag HAVE_SIMD depending on the SIMD_EXTENSION flag given in the installation. Scalargrids work still without SIMD. In the tricubic interpolation Grid3d are downcasted to float-precision and a corresponding warning is given.

Test coverage

The unit tests in testCore.cpp were extended to cover next to the trilinear interpolation also the new interpolations and boundary conditions analogue to the existing tests.

Bug fix

The reflective boundary condition had a bug as described in issue #361 Here, the grid was not centered inside of the tiling cube that was mirrored by the clamping and the fractional distances used as weights in the trilinear interpolation were calculated with unreduced/unmirrored coordinates. Sometimes indices were out of the grid range. We fixed it on the way by shifting the mirroraxis and using the reduced/mirrored coordinates to calculate the fractional distances to the neighbours.

In the following plot all interpolations with both correct boundaries are shown with the example given in #361

Thank you in advance! Leander, @antoniusf, @JulienDoerner, @reichherzerp, @eicheB

[lukasmerten]

After a quick look this makes totally sense to me. I will compile this version tomorrow and report back after running you tests.

You could already now resolve the conflict in the CMakeLists file.

[antoniusf]

I'm sorry for the delay, but as you've requested in the last meeting, I have now taken a look at the simd implementation specifically and checked that the vectorized and unvectorized implementations of tricubicInterpolate compute the same results, up to associativity and commutativity. (Note that this means that the actual numerical results returned by the two implementations may still differ by small amounts, since associativity is not strictly upheld by IEEE-754 floating point operations.)

In addition, of course, the vectorized implementation currently only operates on single-precision numbers, meaning that Vector3ds are downcast to Vector3fs, and all interpolation operations are performed at reduced precision.

Overall, the vectorized implementation looks good to me; according to my analysis, it's correct if the non-vectorized implementation is correct.

[reichherzerp]

@lukasmerten, do you think that the current version of the code looks ready?

[lukasmerten]

I am currently on vacation until November 1st. So I am not sure if I will find time to review things before the first week of November.

[lukasmerten]

@LeanderSchlegel I compiled and tested this version locally. I could not find any strange behavior. Warning messages, when compiling without SIMD extensions are meaningful. Downcasting warning when using Grid3d works, too. There is one last request:

• [ ] Please add two short sentences in the CHANGELOG file under new features and bug fixes.

[LeanderSchlegel]

@lukasmerten Thank you for your testing! Just added the short changelog :)