Random field on a structured grid with periodic boundaries

[joshua-benabou]

Does gstools have support for random field generation on a structured grid with periodic boundaries? This is very useful for e.g reducing finite size effects in simulations.

Consider for example a scalar field in 2D in a NxN grid, such that the field values are, in the limit of large N, continuous across the left-right and top-bottom edges of the grid.

This should amount to a redefinition of the isotropic spatial covariance function C(x1,x2) = C(|x1-x2|) between the random field variables at grid points x1, x2., via C(x1,x2)=C(|x1-x2|_periodic) with |x1-x2|_periodic the distance between these points using periodic boundaries.

Equivalently, at the level of the field's Fourier decomposition, one would just need to sum Fourier modes \vec{k} whose components are integer multiples of 2pi/N.

I suspect this would be either a minor change to summate and summate_incompr in gstools/field/summator.pyx , or to the part of the code that samples from the spectral density distribution of the covariance model, to impose that the sampled wave-vectors k have components integer multiples of 2pi/N.

A more complex use case is spatiotemporal field in 3 space + 1 time dimensions, with periodic boundary conditions in space, but not time. For this all that is needed is the ability to specify which wavevector components should have periodic boundaries (e.g by including an optional boolean array parameter periodic_bc with length given by the model dimension and default value [False,....,False].

I am happy to think more about how to make this work if not already implemented.

[LSchueler]

Random fields on a donut, what a brilliant idea :-)
I think for many geospatial applications this is not of concern. Especially for Lagrangian particle simulations, you can calculate the field values at the particle positions very efficiently, because all the heavy sampling is done in a pre-processing step.
But I would definitely like periodic boundaries.
In case you want this to be your own contribution, please feel free to come up with a solution, put your name into the CONTRIBUTIONS.md, if you like and open a PR. Otherwise, I will come back to this when I have a bit more time.

If these changes introduce many if-branches or computational overhead to the already existing methods, maybe this would be a good first use case of creating a new Generator class?

[joshua-benabou]

Note: I added this feature in the PR # 260.