[WIP] Add support for nd vectors on md field

[LSchueler]

This is a first suggestion of how to implement spatio-temporal vector fields. It is hardly tested and more a foundation for a discussion.
We want to implement 2d vectors (let's call this vec_dim=2) in 2 spatial and 1 temporal dimensions (let's call this field_dim=3)and 3d vectors (vec_dim=3) in 3 spatial and 1 temporal dimension (field_dim=4). This means, that the incompressibility of the vector field only has to be ensured in the spatial plane, which is field_dim - 1 and always equal to vec_dim. This means that I didn't have to change much in the summate_incompr function in summator.pyx. The phase is now computed for all field dimensions, but the summed_modes are computed only for the vec_dim.
As I don't have a lot of time at the moment, I haven't thought this completely through and I'm not sure, if this makes 100% sense.
I'm sure there are also many (corner?) cases which we have to check, e.g. I haven't checked unstructured grids at all.

This is the only script I have used to test this PR:

import numpy as np
import matplotlib.pyplot as pt
import gstools as gs

x, y, z, u = np.arange(10), np.arange(10), np.arange(10), np.arange(10)

# model = gs.Gaussian(dim=4, var=1, len_scale=10)
# srf = gs.SRF(model, generator='VectorField', vec_dim=3, seed=19841203)
# field = srf((x, y, z, u), mesh_type='structured')

model = gs.Gaussian(dim=3, var=1, len_scale=[0.5, 1, 2])
srf = gs.SRF(model, generator='VectorField', vec_dim=2, seed=19841203)
field = srf((x, y, z), mesh_type='structured')


pt.streamplot(x, y, field[0,:,:,0], field[1,:,:,0])
pt.streamplot(x, y, field[0,:,:,1], field[1,:,:,1])
pt.streamplot(x, y, field[0,:,:,2], field[1,:,:,2])
pt.streamplot(x, y, field[0,:,:,3], field[1,:,:,3])
pt.show()

I'm very much looking forward to your thoughts and criticism! @joshua-benabou please join in on this PR and discussion, which arose from discussion #245 .

[joshua-benabou]

Hi Lennart,

thanks much for your code changes. I will pull this branch and do the obvious tests that (i) the field realizations remain incompressible for the spatial components, (ii) the spatial correlation length is always the one specified at each time-step (iii) the correlation time-scale is the one specified.

What corner cases do you have in mind besides unstructured grids (which, for my personal application, I don't need to bother with)?

[joshua-benabou]

Dear Lennart,

I cloned the branch a tested it quickly. Some quick tests by eye show that the correlation length in time and in space match the specified ones, and the 3D spatial field realizations are incompressible. This is really great!

However, there seems to be one issue:

• in your example, trying to impose a mean vector of the field does not work, e.g replacing

srf = gs.SRF(model, generator='VectorField', vec_dim=2, seed=19841203)
field = srf((x, y, z), mesh_type='structured')

by

srf = gs.SRF(model, mean=(-1, 0), generator='VectorField', vec_dim=2, seed=19841203)
field = srf((x, y, z), mesh_type='structured')

gives the error

ValueError: Mean/Trend: Wrong size ([-1. 0.])

This error comes from line 54 of gstools/field/base.py, in the function _set_mean_trend.

I think it is expecting a mean vector of size 3 here (even though that is not what we want as the spatial vectors are 2D here). Trying instead

srf = gs.SRF(model, mean=(-1, 0, 0), generator='VectorField', vec_dim=2, seed=19841203)
field = srf((x, y, z), mesh_type='structured')

gives the error

ValueError: operands could not be broadcast together with shapes (2,10,10,10) (3,10,10,10) (2,10,10,10)

which comes from line 100 of gstools/src/gstools/normalizer/tools.py (function apply_mean_norm_trend)

I think a fix for this is all that is needed to get this working!

[LSchueler]

Hi Joshua,

that is great to hear, thanks for testing the code! I hope I'll find some time this weekend to fix the problem you found.

[joshua-benabou]

Dear Lennart,

After further testing, I believe the incompressibility does not work here :(

Attached is a minimal example showing that the divergence of the field is zero (rather, approaches zero for large grid sizes) for 2d and 3d vector fields, but not for a (3 space + 1 time) field with 3d vectors.

I wrote two functions to calculate the divergence of a 2d and 3d field, respectively, and one function test_incompressible which tests each of the 3 cases I listed above, calculating and plotting the divergence of the field on a suitable spatial slice. The grid size can easily be increased here, which shows that the divergence approaches zero for the 2d and 3d case, but not for the (3+1)d case.

Note that the divergence is always inaccurate at the edges for all cases, but that is expected.

Please let me know if my test is faulty, or if there is actually an issue.

import numpy as np
import gstools as gs
import matplotlib.pyplot as plt


    
def div_3d(F, grid_spacing):
    Fx = F[:,:,:,0] 
    Fy = F[:,:,:,1] 
    Fz = F[:,:,:,2] 

    dFx_dx= np.gradient(Fx, grid_spacing, axis=[0])
    dFy_dy= np.gradient(Fy, grid_spacing, axis=[1])
    dFz_dz= np.gradient(Fz, grid_spacing, axis=[2])
    
    div = dFx_dx + dFy_dy + dFz_dz
    mean_div=np.sum(div)/np.shape(F)[0]**3
    
    print("mean div: ", mean_div)
    
    return div

def div_2d(F, grid_spacing):
    Fx = F[:,:,0] 
    Fy = F[:,:,1] 
 
    dFx_dx= np.gradient(Fx, grid_spacing, axis=[0])
    dFy_dy= np.gradient(Fy, grid_spacing, axis=[1])
    
    div = dFx_dx + dFy_dy
    mean_div=np.sum(div)/np.shape(F)[1]**2
    
    print("mean div: ", mean_div)
    
    return div

def test_incompressible():
    dim_x = 40
    dim_t = 5
    cor_x = 0.1
    cor_t = 0.3
    
    x = np.arange(dim_x)/dim_x
    y = np.arange(dim_x)/dim_x
    z = np.arange(dim_x)/dim_x
    t = np.arange(dim_t)/dim_t

    grid_spacing = 1/dim_x

    #############################################################################################
    # 2d, it works
    print("\n2d incompressible field...")

    model_2d = gs.Gaussian(dim=2, var=1, len_scale=cor_x)
    srf_2d = gs.SRF(model_2d, generator='VectorField',mean=(-1, 0))
    srf_2d((x, y), mesh_type='structured')
    F_2d = srf_2d.field

    div_arr_2d=div_2d(np.transpose(F_2d, (1,2,0)),grid_spacing)
    
    print("\nplotting divergence")

    plt.figure()
    plt.pcolormesh(div_arr_2d)
    plt.colorbar()
    ###########################################################################################
    # 3d, it works

    print("\n3d incompressible field...")
    model_3d = gs.Gaussian(dim=3, var=1, len_scale=cor_x)
    srf_3d = gs.SRF(model_3d, generator='VectorField',mean=(-1, 0,0))
    srf_3d((x, y,z), mesh_type='structured')
    F_3d = srf_3d.field

    div_arr_3d=div_3d(np.transpose(F_3d, (1,2,3,0)),grid_spacing)

    print("\nplotting divergence of z-slice")
    plt.figure()
    plt.pcolormesh(div_arr_3d[:,:,dim_x//2])
    plt.colorbar()
    
    ##############################################################################################
    # 4d, it fails
    
    print("\n4d field of 3d vectors which is supposed to be (spatially) incompressible ...")

    model_4d = gs.Gaussian(dim=4, var=1, len_scale=[cor_x,cor_x,cor_x,cor_t])
    srf_4d = gs.SRF(model_4d, generator='VectorField', vec_dim = 3) # !!! note that I don't specifiy mean here since it doesn't work
    srf_4d((x, y, z, t), mesh_type='structured')
    F_4d = srf_4d.field
    
    print(np.shape(F_4d))
    div_arr_3d=div_3d(np.transpose(F_4d[:,:,:,:,2], (1,2,3,0)), grid_spacing) # select 3rd time-slice

    print("\nplotting divergence of (z,t)-slice")
    plt.figure()
    plt.pcolormesh(div_arr_3d[:,:,dim_x//2])
    plt.colorbar()

    
    
#######################################################
# run test
    
test_incompressible()

[LSchueler]

That's a pity! But thank you so much for your script.

Regarding the mean value, actually I'm not sure at the moment how to interpret it when applied to velocity fields. A lot of additions came to the code with scalar fields in mind. The way to set the mean velocity in x-direction (as you can always rotate your problem in such a way that we do not loose on generality... except of course you want 0 mean velocity ;-) ), is by setting srf.generator.mean_u.

I'll need some time to see, if I can find a solution.

[LSchueler]

Ah, by the way, in case you haven't done this already, you can significantly decrease computational time by using this line before compiling/ installing GSTools:
export GSTOOLS_BUILD_PARALLEL=1

[LSchueler]

I just calculated the divergence of the equation given in the notes field here and the divergence is actually only zero, if the dimension of the cov. model coincides with the dimension of the vectors :-( I guess I should have done that calculation before hacking away ;-)

[LSchueler]

Maybe I just need more sleep or more coffee...
I just realised that I can make a tiny edit to the equations and it works...
I think it should work now. Can you still somehow use this, with non-zero mean velocity? What about subtracting the vector (-1, 0, 0, ...) after the random field generation?

[joshua-benabou]

Dear Lennart,

Thank you for your excellent and rapid communication on this issue.

Let me pull this again and check that everything works; your fix seems quite logical! As far as the mean velocity, this is not too problematic; previously i have been commenting out line 448 in generator.py:

#self.mean_u * e1

Since the default is mean_u=1, this gives me a zero-velocity fluid.

[joshua-benabou]

Dear Lennart,

I confirm that your fix worked: when generating 3D vector fields in (3 space + 1 time) dimension, the field realizations are now (spatially) incompressible (awesome!), although with zero curl in the x-direction.

For my application, I desire a fluid of zero mean velocity. The zero velocity is easy to get by commenting out the default mean velocity in the x-direction that your code automatically imposes, but the incompressibility projector also needs to be modified.

I plan to make two minor code contributions to treat the zero-velocity case:

(i) simply remove the projector for the incompressible vector field generator to get realizations which are not incompressible, but which can be transformed to an incompressible field by then taking the curl

(ii) i will try to implement the method suggested by (Kraichnan, 1970) Eq.19, 20 who appears to have treated the zero-velocity case with a more general projector then what you are using. It appears to require more computation though, because each wave coefficient must be drawn for a 2 or 3 dimensional multivariate Gaussian. I suppose your projector is a sort of special case? If you know of any code implementations (available on github) of the zero-velocity case, I am happy to use those as a start point for contributing to your code.
kraichnan1970.pdf

[joshua-benabou]

Dear Lennart,

I filed a new pull request (#260) for this issue.

[LSchueler]

Dear Joshua,

I'm sorry for taking so long to come back the 4d fields. The last 10 days or so I just had too much things to do. Tomorrow I will have time to take a look at the PR (I live in central European time), I'm very much looking forward to it!