3

I have a 3D vector field and I want to extract its divergence-free part (also called transverse component), using the Helmholtz decomposition. In principle, this can be done in the Fourier space, as illustrated here in the proof (and also reported here).

I wrote this Python function, implementing this straightforward calculation, where the input is a vector field of shape (3,N,N,N), since it has 3 components and is defined in a 3D square grid with NxNxN nodes and size equal to boxsize:

def helmholtz_decomposition(vector_field, boxsize):
    kx, ky, kz = np.meshgrid(2 * np.pi * np.fft.fftfreq(vector_field.shape[1], boxsize / (vector_field.shape[1])),
                             2 * np.pi * np.fft.fftfreq(vector_field.shape[2], boxsize / (vector_field.shape[2])),
                             2 * np.pi * np.fft.fftfreq(vector_field.shape[3], boxsize / (vector_field.shape[3])),
                             indexing='ij')

k_vector = np.array([kx, ky, kz])
k_squared = kx ** 2 + ky ** 2 + kz ** 2
k_squared[k_squared == 0] = 1e-12  # Avoid division by zero

vector_field_fourier = np.fft.fftn(vector_field)

dot_product = kx * vector_field_fourier[0] + \
              ky * vector_field_fourier[1] + \
              kz * vector_field_fourier[2]

potential_part_fourier = (dot_product / k_squared) * k_vector
solenoidal_part_fourier = vector_field_fourier - potential_part_fourier

potential_part = np.real(np.fft.ifftn(potential_part_fourier))
solenoidal_part = np.real(np.fft.ifftn(solenoidal_part_fourier))

return potential_part, solenoidal_part

Here, the potential and solenoidal parts should be curl-free and divergence-free respectively.

The problem is that the divergence of the solenoidal part is not zero. I have already checked that the solenoidal part in the Fourier space is perpendicular to the wave vector k_vector, as it should be.

For example: the divergence of the original total vector field is 0.0111, while the divergences of its solenoidal and potential parts are 0.00584 and 0.00555 respectively.

Does anyone know where the mistake is?

Thanks a lot for your help!

Wil
  • 31
  • 2
  • Welcome to scicomp! Is the number of real space mesh points equal to the number of cells? Or do you have n+1 DOF for your n cells? – MPIchael Aug 29 '23 at 11:55
  • Thank you @MPIchael . The number of grid points is equal to the length of the array, if that's what you mean. – Wil Aug 29 '23 at 11:58
  • I am not familiar with the details in python/numpy, but many errors in fourier algorithms stem from missing or doubled normalization. In your example you scale by 2*pi? See for example: https://dsp.stackexchange.com/questions/66058/am-i-supposed-to-normalize-fft-in-python – MPIchael Aug 29 '23 at 14:04
  • Unfortunately, the 2*pi factor doesn't change the result. But I will think about it again. Thanks a lot @MPIchael. – Wil Aug 29 '23 at 14:35

0 Answers0