I'm working on a code to simulate turbulence in a periodic cube. To do so I start with a velocity vector field of shape [64,64,64,3] which means 64 grid points in each direction (x,y,z) and 3 velocity components at each point in this space (ux,uy,uz). To calculate the next time step the vorticity field must be calculated from the velocity field and here the problem appeared.
Let us for example consider the so called ABC-flow of the form
In this case the rotation of the velocity field stays the same, i.e. velocity field u and vorticity field omega are equal. To calculate the velocity field out of the vorticity field one can use an equation as follows
To calculate the derivatives I used the pseudospectral method such that the equations read like this
.
With u_k and omega_k being the velocity and vorticity field in Fourier space.
I tried to write a code to calculate the vorticity from the ABC-flow velocity field in Fourier space and calculate the velocity field from the vorticity field afterwards . The velocity and vorticity field should look the same in real space but as you can see here they don't.
The picture shows the real part of the first vector component at a single layer out of the cube in real and Fourier space for the velocity at the beginning (pre), the vorticity and the velocity calculated out of the vorticity (post). As you can see the velocity fields look similar as expected while the vorticity field in real space exhibits a stripe like pattern although it should look the same.
Having a look at the imaginary part of the vorticity field in real space shows a stripe like pattern which fits to the one of the real part as you can see here.
My problem now is to find the reason for that pattern building. The vorticity should look like the velocity fields but it doesn't and I cant figure out why. Solving this is very important to avoid a loss of information while calculating the next time step.
My code looks like this:
import numpy as np
import matplotlib.pyplot as plt
xx = np.arange(64) #creating 1D x coordinate
x,y,z = np.meshgrid(xx,xx,xx) #calculating 3D x,y,z coordinates
a = np.zeros(shape=(64,64,64,3)) #creating 3D vector field 'a'
A, B, C = 3., 3., 3. #setting values A,B,C
a[:,:,:,0] = A*np.cos(y/(2.*np.pi)) + B*np.sin(z/(2.*np.pi))
a[:,:,:,1] = B*np.cos(z/(2.*np.pi)) + C*np.sin(x/(2.*np.pi))
a[:,:,:,2] = C*np.cos(x/(2.*np.pi)) + A*np.sin(y/(2.*np.pi))
k = np.zeros(shape=a.shape) # creating k array
k_sqr = np.zeros(shape=a.shape) # creating k^2 array
k_lin = np.fft.fftfreq(64)*2.*np.pi # calculating 1D k coordinate
k[:,:,:,0], k[:,:,:,1], k[:,:,:,2] =
np.meshgrid(k_lin,k_lin,k_lin,indexing='xy') # calculating 3D kx,ky,kz
k_sqr[:,:,:,0] = k[:,:,:,0]**2 + k[:,:,:,1]**2 + k[:,:,:,2]**2
k_sqr[:,:,:,1] = k_sqr[:,:,:,0]
k_sqr[:,:,:,2] = k_sqr[:,:,:,0]
k_sqr[0,0,0,:] += 1. # set zero of k_sqr = 1. to avoide division by 0
a_fft = np.fft.fftn(a,axes=(0,1,2)) # calculate FFT of a (=uk)
wk = np.cross(1j*k,a_fft) # calculating vorticity wk by wk = ik x uk
w = np.fft.ifftn(wk,axes = (0,1,2)) # calculating iFFT von wk (w(x))
uk = np.cross(1j*k,wk)/k_sqr # calculating uk by uk = (ik x wk)/k^2
uk[0,0,0,:] = a_fft[0,0,0,:]
# set value of position of 0-Division to the value it was before
u = np.fft.ifftn(uk,axes = (0,1,2)) # calculating iFFT of uk (u(x))
I didn't find any information about this problem yet. If this question exists already I'm sorry. If not I'm grateful for any ideas.
