To use this, set seed_density2 for your initial conditions

TODO: think about cylindrical case (can we have a divergence free background
field?)

To make sure it does not affect performance

And check whether the output folder is writable

make sure 3D variant is also deterministic (with same random numbers)

But don't use the geometry to restrict the solution.

And re-use previous solution

This should lead to a smaller error in the gradient of the computed solution.
The new approach probably does not work well for cases with jumps in a
coefficient; the old method is still available.