2

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.

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.

Does anyone know where the mistake is?

Thanks for your help!

Wil
  • 55
  • 5
  • Is it close to zero, like accumulated round-off error close? – Tim Roberts Aug 25 '23 at 23:01
  • @TimRoberts unfortunately not, because it's even larger than the divergence of the potential component.. – Wil Aug 26 '23 at 12:08
  • Look, honestly I'm out of my depth here, but I noticed some things. The proof requires you to take an integral of the FFT on the third derivative and multiply by 1/2pi^(3/2). I don't see any of those steps. – Tim Roberts Aug 26 '23 at 17:29
  • @TimRoberts the integrals are performed by the numpy functions fft and ifft: these are the Fourier and the inverse Fourier transforms respectively. Why the third derivative? Thanks for your help – Wil Aug 26 '23 at 21:13
  • You're quite right about the integration. The formula in the proof mentions `d3r` (with the 3 as an exponent), which usually means "third derivative of r". Again, however, I'm out of my depth. – Tim Roberts Aug 27 '23 at 00:17
  • @TimRoberts the 3 as exponent means an integration over the 3D space, since r has 3 components. Thank you anyway for your time!! – Wil Aug 27 '23 at 10:54

0 Answers0