[WIP] Splice multilevel

[cphyc]

This is a work in progress to implement the splicing operation for multi-level setup.

[cphyc]

As of 856f1b0, this is what the result of splicing a field a into another field b, obtaining a spliced field f. Shown here is a 1D splice through the center of the spliced region, which is a sphere of radius 9 Mpc in a zoom of size 25 Mpc for a box size of size 50 Mpc. The primary grid has 128³ pixels, the zoom has the same number of resolution elements and a physical span that's half that.

[cphyc]

Note, in 8ba68c4 and 5e136b5 I had to update the answer test, as the new splicing is slightly different from the old one.

Old Answer, diff between spliced field and inner field	New answer	Difference old/new answer

[cphyc]

Note for the future: when converting from delta space to whitenoise and vice-versa, we loose the fundamental mode (which is set to 0). This result in the spliced field to always have a mean value of 0 (for each grid) and causes it to differ by a constant from the target field. A workaround would be to prevent the conversion to happen in the first place, i.e. to do the splicing at the location where the conversion from whitenoise to delta is done.

[apontzen]

Note for the future: when converting from delta space to whitenoise and vice-versa, we loose the fundamental mode (which is set to 0). This result in the spliced field to always have a mean value of 0 (for each grid) and causes it to differ by a constant from the target field. A workaround would be to prevent the conversion to happen in the first place, i.e. to do the splicing at the location where the conversion from whitenoise to delta is done.

Shouldn't the filtering make this irrelevant? All long-wavelength modes should be coming from the parent grid, not the zoom

[cphyc]

My understanding is that, in non-spliced situation, the fields are in whitenoise space until a call to write or done is performed. At this stage, the fields are converted into delta space (by applying applyPowerSpec) and then the different levels are recombined (through a cascade of calls starting from ensureParticleGeneratorInitialised all the way to initialiseParticleGeneratorBasedOnTemplate).

For splicing, we solve for alpha and then combine all the levels in delta space to form b, M (a - b) and Mbar alpha. The solution to the splicing problem is then f = b + M(a-b) + Mbar alpha, where f is now a combined field in delta space.
As it currently stands, the splicing is done before the call to either write or done, so I short-circuited the recombination step and converted back into whitenoise (but with combined fields, https://github.com/cphyc/genetIC-1/blob/6a4da19fa6c8fd09e280ce139f5441146199223a/genetIC/src/simulation/modifications/splice.hpp#L343). This almost does the trick, but the issue is that this drops the k=0 mode (which was non-null at the finest level because of the combination operation).

Overall, this leads to a small offset between the spliced field and the target field in the zoomed region, as can be seen below. The difference should be of the order of 1e-14, but is instead exactly the mean value of the field f if I compute it before converting to whitenoise.

Note that if I force the k=0 mode to be conserved (by setting it to 1 instead of 0), then the results are perfect down to machine precision.

[apontzen]

My point is that the k=0 should not be contributing, because it should be filtered in from the higher level grids. So I don't quite understand this issue.

I am also nervous about doing splicing after everything else is complete. I have a sense that this will lead to confusing behaviour -- e.g. what if you try to modify something that you spliced?

Ideally we just rewrite the splicing operation to apply in noise space, as I did in the original one-level case.