Add feature for zero mean velocity incompressible field

[joshua-benabou]

The IncomprRandMeth cannot deal with zero mean velocity fluids, since the projector used to ensure incompressibility has a preferential direction along the x-axis, which ensures the generated vector field has curl with zero x-component. For a zero-mean velocity fluid, there is no preferential direction, so we can instead use the general method introduced by Kraichnan (1970).

The main changes here are the addition of summate_incompr_zero_vel and the corresponding class IncomprRandZeroVelMeth . the coefficients of the cosine and sine Fourier-space basis functions are given by the cross-product of a given wave vector of dimension vec_dim and a vector drawn from a multivariate normal distribution of dimension vec_dim.

Please use my script included in Issue #258 to check that field realizations indeed have zero-divergence on average, in the limit of large grid-sizes. I explicitly verified that the curl of the generated field is non-zero in the x-direction (can attach test script if desired).

My changes here are just to show what should be done and probably don't respect the spirit of your code. For example there are clearly a few issues I did not have time to work out: since the super method of generator.py only draws z_1 and z_2 for a univariate Gaussian, I overwrite them inside IncomprRandZeroVelMeth, which doesn't respect the seed business you guys have set up to ensure reproducibility.

Further, I attempted to use a np.cross to perform the cross product in the summate_incompr_zero_vel lines 112,113, but this led to the method running indefinitely when called. I am guessing this is due to type instability or something C-related, but since I'm not used to Cython, I'm not sure. In the meantime I commented out those lines and hard-coded the cross-product for the case vec_dim = 3only.

Further, I paid no attention to the possible normalization of the field (if any), i.e some constant prefactor to ensure the inputted mean velocity of IncomprRandZeroVelMeth is correct. For my personal application, I just rescale the velocity anyway, but I can do this properly once I have a bit more time.

I also added the function summate_generic_vector_field and correspondingly GenericRandVectorFieldMeth . These are supposed to be for generating vector fields which are not necessarily incompressible. The point is that, one may then take the curl to get a field with zero divergence. These functions are not correct yet, so please ignore them for the moment.

Please excuse my surely incomplete pull request, it's my first time contributing to open source code, but I'm happy to learn from mistakes here.

[joshua-benabou]

A quick update: I added a feature for periodic boundary conditions in the spatial dimension (first N coordinates where N=vec_dim), but only for the generator dealing with zero-velocity incompressible random fields. However one could trivially copy paste my code for periodic b.c here to use for scalar and generic incompressible random fields.

All it does is, for a given dimension, only use Fourier modes which are multiples of 2pi/box_length, where box_length is the length of the simulation box in that dimension (which the user needs to input for now). I don't touch the spectral density of the generated modes, but just round the generated Fourier modes to the nearest of these multiples.

I explicitly checked the periodic b.c setup works, can share code for this check later.

[LSchueler]

Please delete this line :-)

[LSchueler]

Hi Joshua,

I think you forked off the main branch or already merged/ rebased changes from there to your local main-branch, which are not on the nd-vector-fields branch 😅
The diff looks really messy. I merged all the commits from the main branch into the nd-vector-fields branch, but that didn't help at all.
Would it be a lot of work on your side to make a new fork off the nd-vector-fields branch and recommit your changes?

[LSchueler]

Are you sure that the seed does not work? - You are using the RNG class for your z1 and z2 so I guess you should be able to create reproducible fields with a given seed.

[LSchueler]

I think these are nice debug prints, but for production code, they should be removed or at least you should use the verbose keyword, as is done in elsewhere in this file.

[LSchueler]

Using numpy-functions in Cython code is very very slow. I don't think that you ran into numerical problems, you just didn't wait long enough ;-)
I think it would be general enough for the time being if you could provide the curl operation in 2d and 3d.

[LSchueler]

Could you add your checks of the divergence and so on to the tests folder? - Preferably very small test cases which run fast.
That would increase the reliability of your code and any future changes tremendously.

[LSchueler]

I'm sorry, I'm running out of time. I will continue with the review tomorrow. I hope that with my few hints and suggestions, you can continue to improve your code!

[joshua-benabou]

Dear Lennart,

Sorry for the long delay in my reply. I finally had time to come back to this and aim to make these changes shortly. However, in the process of testing things I encountered an issue relating to the OpenMP functionality.

I've thus far used pip install -e . in the main package dir for my local changes. Now, I set the env variable export GSTOOLS_BUILD_PARALLEL=1 . Then doing pip install -e . I get among other things in the terminal ...

Installing collected packages: gstools
Running setup.py develop for gstools
Successfully installed gstools-1.3.6.dev37

The issue is nothing actually changed because, setup.py (shown below) is supposed to print "OpenMP=True" if the env variable is set to GSTOOLS_BUILD_PARALLEL=1 in the linux terminal , and print something else if its not set to 1.

I tried instead just python setup.py install but that gives

UNKNOWN 0.0.0 is already the active version in easy-install.pth

Installed /global/u1/b/benabou/.conda/envs/healpy_conda_gstools_dev/lib/python3.8/site-packages/UNKNOWN-0.0.0-py3.8-linux-x86_64.egg
Processing dependencies for UNKNOWN==0.0.0
Finished processing dependencies for UNKNOWN==0.0.0

and import gstools

no longer works correctly.

So how can I install my edited version of the package with OpenMP support?

[LSchueler]

Hi Joshua,

it's great to hear from you again! I think pip surpresses all output from within setup.py. You can try passing the verbose flag like pip install -e . -v and then you should see that OpenMP is actually probably being used with export GSTOOLS_BUILD_PARALLEL=1.
You can always run a simple but expensive calculation, like generating a big 3d-field and watch your systems monitor to see how many cores are being used.
I hope that works 🤞