Description
Hello,
I believe I found a bug with the incompressible random vector field algorithm. Generating a 3d incompressible field A (3 space dimensions, 3-component vectors), I find that the quantity dA_z/dy-dA_y/dz is always zero in the limit of large grid sizes.
A priori there is no reason this should be true for a generic incompressible field, and is unphysical, as it means the curl of A always has x-component equal to zero. This could be an artefact of the projector used to ensure incompressibility. The projector appears to introduce a preferential direction along the x-axis...
So as is, the incompressible random vector field code is not able to simulate an incompressible fluid with non-zero curl along the x-axis. For a fluid with zero mean velocity this is especially problematic.
If this was by design, it would be very useful to have an algorithm which generates a random 3d vector field which is not incompressible, so that an incompressible one may be obtained by then taking the curl. This is as simple as removing the projector in the summate_incompr method. Edit: Actually, I'm wondering if there is any downside with this method? For a Gaussian field, the curl of the field should have the same statistics (e.g correlation length) as the field itself (don't know a proof of this statement, but seems plausible).
I would do this myself, however I would need some help understanding your workflow. I am guessing something like: write cython in the summator.pyx file and then compile summator.c which I think is done automatically in setup.py. However must I rerun 'setup.py' everytime after making changes? Some details on your workflow here so I can try to make some contributions would be great!
Here is the minimal example showing the issue:
from B_field_generation import *
import matplotlib.animation as animation
import matplotlib.pyplot as plt
import numpy as np
import h5py
import gstools
import gstools.covmodel
import time
print("Done!")
############################################################################################
# generate incompressible random vector field in 3D
dim_x=40
x, y, z = np.arange(dim_x)/dim_x, np.arange(dim_x)/dim_x, np.arange(dim_x)/dim_x
cor_x=0.1
model = gs.Gaussian(dim=3, var=1, len_scale=[cor_x, cor_x, cor_x]) # has mean zero
srf = gs.SRF(model, generator='VectorField') # incompressible vector field #, seed=19841203
start = time.time()
A_field = srf((x, y, z), mesh_type='structured') # draw 1 realization for all our time points
end = time.time()
# get field components and derivatives
Az = A_field[2,:,:,:]
dAz_dy = np.gradient(Az, axis=1)
Ay = A_field[1,:,:,:]
dAy_dz = np.gradient(Ay, axis=2)
# plot field components
plt.pcolormesh(Az[:,:,dim_x//2])
plt.title("Az")
plt.figure()
plt.pcolormesh(Ay[:,:,dim_x//2])
plt.title("Ay")
plt.figure()
plt.pcolormesh(dAy_dz[:,:,dim_x//2])
plt.title("dAy_dz")
plt.figure()
plt.pcolormesh(dAz_dy[:,:,dim_x//2])
plt.title("dAz_dy")
plt.figure()
plt.pcolormesh(dAz_dy[:,:,dim_x//2]-dAy_dz[:,:,dim_x//2]) # this approaches zero for large grid sizes
plt.title("dAz_dy-dAy_dz")