Why your energy spectrum is incorrect?

In turbulence simulations, one of the most interesting quantities is the energy spectrum of turbulence. Thus, getting the correct energy spectrum for a given velocity field is crucial for analyzing the turbulence data. However, there are many subtleties in calculating the energy spectrum, which can lead to incorrect results if not handled properly. In this notebook, we will discuss how to calculate the energy spectrum of turbulence correctly using Fourier transforms and binning techniques. We will also provide a code example to illustrate the process.

Turbulence Energy Spectrum

The energy spectrum of turbulence $E(k)$ is a function on the wavenumber $k$, which represents the distribution of kinetic energy across different scales in the turbulent flow. The turbulent kinetic energy (TKE) is defined as “the mean kinetic energy per unit mass in the fluctuating velocity field”[?]. Given a fluctuating velocity field $\mathbf{u}=[u_x, u_y, u_z]$, the TKE can be calculated as:

The energy spectrum $E(k)$ is related to the TKE by the following integral:

Analytical Calculation of Energy Spectrum

In most textbooks, there exists an analytical formula for calculating the energy spectrum from the velocity field. However, this formula is often misused in many postprocessing codes, which can lead to incorrect results. In this section, we will first derive the analytical formula for calculating the energy spectrum and then discuss the common pitfalls in its implementation.

Before we start, let’s start with a discussion on the Fourier transform. Although it is already a well-established mathematical tool, there exist many different conventions for defining the Fourier transform, which can lead to confusion when implementing the energy spectrum calculation. Below, we list three common conventions for the Fourier transform:

1. Fourier transform in ordinary frequency:

1. Fourier transform in angular frequency / wave number, unitary convention: \(\hat{u}(k)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}u(x)e^{-kxj}dx.\)

1. Fourier transform in angular frequency / wave number, non-unitary convention:

To derive the energy spectrum, we will need to use Parseval’s theorem, which states that the total energy in the spatial domain is equal to the total energy in the frequency domain. For different Fourier transform conventions, Parseval’s theorem can be expressed as:

1. Ordinary frequency:

1. Wavenumber and unitary definition:

1. Wavenumber and non-unitary definition:

If we use the unitary definition of the Fourier Transform with the wavenumber, we can transform the TKE into the frequency domain as:

Note that the right-hand side is a function of vector wavenumber $\mathbf{k}=[k_x, k_y, k_z]$, while the energy spectrum is a function of scalar wavenumber $k=|\mathbf{k}|$. To transform the energy from vector wavenumber space to scalar wavenumber space, we can turn the volumetric integration on $\mathbf{k}$ to the integration on $k$:

where $A(k)$ is the unit area of a sphere shell with radius $k$. If the integration space is a unit space, i.e., ${\int\int\int dxdydz}=1$, we will have

I.e.,

This is the most common formula of the TKE spectrum function in textbooks. However, it is important to notice that there are two conditions for the above equation:

1. We use the unitary definition of the Fourier Transform. Using a non-unitary definition will introduce additional coefficients.
2. The volume of the physical space should be equal to 1, i.e., ${\int\int\int dxdydz}=1$

Numerical Calculation of Energy Spectrum

Although we have the analytical formula for calculating the energy spectrum, it is not straightforward to implement it numerically. In this section, we will discuss how to calculate the energy spectrum numerically using the Discrete Fourier Transform (DFT).

In most framework like numpy and PyTorch, the DFT is defined as:

It is also possible to use other definitions of DFT by setting the normalization coefficient and the sign of the exponent (the norm parameter is available in numpy to change the normalization coefficient). If you don’t set the norm parameter, which is the most common case, the DFT and its inverse will be defined as above. In this case, Parseval’s theorem can be expressed as:

If we want to use the wavenumber $k=2\pi f$ rather than the frequency $f$, we simply replace $f$ with $k/2\pi$ and sum over all $k$s, i.e.,

You may notice that the additional coefficient $\frac{1}{N}$ is introduced in Parseval’s theorem here, which already indicates that the energy spectrum calculated using DFT will be different from the one calculated using the continuous Fourier transform.

In DFT scenarios, the TKE is calculated as

Here, we can split the sum into two parts just like what we have done in the continuous form of TKE spectrums. We can define a new value called “energy in shell $E_s$” where we sum all the values which have the same $k=\sqrt{k_x^2+k_y^2+k_z^2}=\sqrt{(2\pi f_x)^2+(2\pi f_y)^2+(2\pi f_z)^2}$. Thus, we will have:

The energy in shell function $E_s(k)$ is different than the energy spectrum function $E(k)$ as $TKE=\sum_k E_s(k)=\sum_kE(k)\Delta k$.

Since the $\Delta k$ is not uniformly distributed in 3D space, we usually need to define the $Delta k$ by ourself, i.e., choose the bin size for $k$ and sum all the values in the same bin together to get the energy in shell function $E_s(k)$ and devide it by the bin size to get the energy spectrum function $E(k)$.

Now, let’s summarize the key points for calculating the energy spectrum of turbulence in numerical ways:

1. You shouldn’t use the analytical formula for calculating the energy spectrum directly, as it is derived based on the continuous Fourier transform and may not hold in discrete scenarios.
2. You should use the correct definition of DFT and its corresponding Parseval’s theorem to calculate the TKE in the frequency domain. For the most common implementation of DFT, an additional coefficient of $\frac{1}{(N_xN_yN_Z)^2}$ will be introduced in the TKE calculation.
3. You should define the bin size for $k$ by yourself, i.e., choose the bin size for $k$ and sum all the values in the same bin together to get the energy in the shell function $E_s(k)$ and divide it by the bin size to get the energy spectrum function $E(k)$.

Code Example

Based on the above discussion, we can implement the energy spectrum calculation in code as follows:

from torchfsm import FourierMesh
from torchfsm._type import SpatialTensor
import torch,math

def get_energy_spectrum(u:SpatialTensor,energy_coeff=1.0):
    """
    Calculate the energy spectrum E(k) given a velocity field $\mathbf{u}$ which satisfies
    $$\int E(k) dk =\int c|\mathbf{u}|^2 d\mathbf{x}$$,
    where c is the energy coefficient.
    The bin size of the wavenumber k is 1, and the center of the first bin is 0.5. 
    This also makes the returned Energy spectrum E(k) the same as the energy in shell $E_s(k)$, 
    satisfing $\frac{\sum c|\mathbf{u}|^2 dV}{\sum dV} = \sum E_s(k) = \sum E(k) \Delta k$.

    Args:
        u (SpatialTensor): The input velocity field, with shape [1, C, Nx, Ny, Nz].
        energy_coeff (float, optional): The energy coefficient c. Defaults to 1.0.
    
    Returns:
        k_bins (torch.Tensor): The center of the wavenumber bins.
        Ek (torch.Tensor): The energy spectrum E(k) for each bin.


    """
    assert u.shape[0] == 1, "The batch size must be 1 for energy spectrum calculation."
    spatial_shapes=u.shape[2:]
    f_mesh = FourierMesh([(0, 1, spatial_shape_i) for spatial_shape_i in spatial_shapes])
    u_fft = f_mesh.fft(u)
    energy_fft = torch.sum(energy_coeff*torch.abs(u_fft)**2, dim=1)/(math.prod(spatial_shapes)**2)
    # Use the `long` type will auto convert the float values to the corresponding integer bin index with bin size 1.
    k_long=torch.norm(f_mesh.bf_vector*2*torch.pi,dim=1).view(-1).long()
    # `torch.bincount` will sum the energy values in the same bin together, and return the energy spectrum with bin size 1.
    Ek=torch.bincount(k_long, weights=energy_fft.view(-1))

    # We define the actual wavenumber of each bin to be the center of the bin, which is 0.5, 1.5, 2.5, ...
    k_bins = (torch.arange(Ek.shape[0]) + 0.5)
    return k_bins, Ek

We use TorchFSM, a PyTorch-based Fourier Spectral Method library for PDE simulations, to perform the Fourier transform and calculate the energy spectrum. You can also replace it with some other libraries like numpy. In the above code, we set the bin size for $k$ to be 1, which means the final returned spectrum values can be both treated as the energy in the shell function $E_s(k)$ and the energy spectrum function $E(k)$, as the $\Delta k$ is 1. If you want to get the energy spectrum function $E(k)$, you can simply divide the returned values by the bin size.

We can verify the correctness of the above code by comparing the energy calculated in the spatial domain and the energy calculated in the frequency domain. If they are close enough, we can conclude that the energy spectrum calculation is correct.

u=torch.randn(1,3,128,128,128)
k_bins, Ek = get_energy_spectrum(u)
Ek.sum(),torch.mean(torch.sum(u**2, dim=1))

(tensor(2.9969), tensor(2.9969))

References

1: Pope, S. B. (2000). Turbulent Flows, p88. Cambridge: Cambridge University Press.