Grids

1D, 2D, and 3D grids are supported. We demonstrate here the construction of a one-dimensional grid and how one can use it to perform Fourier transforms and compute spatial derivatives.

A one-dimensional grid with $n_x = 64$ grid points and length $L_x = 2 \pi$ is constructed by

julia> using FourierFlows

julia> nx, Lx = 64, 2π;

julia> grid = OneDGrid(nx, Lx)
OneDimensionalGrid
  ├─────────── Device: CPU
  ├──────── FloatType: Float64
  ├────────── size Lx: 6.283185307179586
  ├──── resolution nx: 64
  ├── grid spacing dx: 0.09817477042468103
  ├─────────── domain: x ∈ [-3.141592653589793, 3.0434178831651124]
  └─ aliased fraction: 0.3333333333333333

The grid domain is, by default, constructed symmetrically around $x = 0$, but this can be altered using the x0 keyword argument of the OneDGrid constructor. The grid spacing is $L_x / n_x$. Note that the last point of the domain is a grid-spacing before $L_x / 2$. This is because periodicity implies that the values of any field at the end-points of the domain are equal and, therefore, grid-point values at both these end-points are reduntant.

We can define an array u that contains the values of a function $u(x)$ on this grid as

u = @. sin(2 * grid.x) + 1/2 * cos(5 * grid.x)

Note that we chose a function that is periodic on our domain. We can visualize u by

using Plots

plot(grid.x, u, label="u", xlabel="x")

Function $u(x)$ can be expanded in Fourier series:

\[u(x) = \sum_{k} \hat{u}(k) \, e^{i k x} ,\]

where $\hat{u}(k)$ is Fourier transform of $u(x)$ and $k$ the discrete set of wavenumbers that fit within our finite domain. We can compute $\hat{u}$ via a Fast Fourier Transform (FFT). Since our u array is real-valued then we should use the real-FFT algorithm. The real-valued FFT transform only saves the Fourier coefficients for $k \ge 0$; the coefficients for negative wavenumbers can be obtained via $\hat{u}(-k) = \hat{u}(k)^{*}$.

The wavenumbers used in FFT are contained in grid.k and they are ordered as:

\[\frac{2\pi}{L_x} \{ 0, 1, \dots, n_x/2-1, -n_x/2, -n_x/2+1, \dots, -1 \} ,\]

while the wavenumbers for real FFT are in grid.kr:

\[\frac{2\pi}{L_x} \{ 0, 1, \dots, n_x/2-1 \} .\]

The grid also includes the FFT plans for both real-valued and complex valued transforms:

grid.fftplan
FFTW forward plan for 64-element array of Complex{Float64}
(dft-direct-64 "n2fv_64_avx2_128")
grid.rfftplan
FFTW real-to-complex plan for 64-element array of Float64
(rdft2-ct-dit/2
  (hc2c-direct-2/2/0 "hc2cfdftv_2_sse2"
    (rdft2-r2hc-direct-2 "r2cf_2")
    (rdft2-r2hc01-direct-2 "r2cfII_2"))
  (dft-direct-32 "n1fv_32_avx2_128"))

We use the convention that variables names with h at the end stand for variable-hat, i.e., $\hat{u}$ is the Fourier transform of $u$ and is stored in array uh. Since u is of size $n_x$, the real-Fourier transform should be of size $n_{kr} = n_x/2+1$.

uh = Complex.(zeros(grid.nkr))

The FFT transform is done as an in-place matrix multiplication using mul!.

using LinearAlgebra: mul!

mul!(uh, grid.rfftplan, u)

The FFT algorithm does not output exactly the Fourier coefficients $\hat{u}(k)$ but rather due to different normalization, FFT outputs something proportional to $\hat{u}(k)$. To obtain $\hat{u}$ we need to divide the FFT output by the length of the original array and by $e^{-i k x_0}$, where $x_0$ is the first point of our domain array.

uhat = @. uh / (nx * exp(- im * grid.kr * grid.x[1])) # due to normalization of FFT

plot(grid.kr, [real.(uhat), imag.(uhat)],
          label = ["real( û )" "imag( û )"],
         xlabel = "k",
          xlims = (-0.5, 10.5),
          ylims = (-0.55, 0.55),
         xticks = 0:10,
         yticks = -0.5:0.25:0.5,
         marker = :auto)

We can compute its derivative via Fourier transforms. To do that we can use the FFTW plans that are constructed with the grid. First we allocate some empty arrays where the values of the derivative in physical and Fourier space will be stored,

∂ₓu  = similar(u)
∂ₓuh = similar(uh)

The grid contains the wavenumbers (both for real-value functions grid.kr and for complex-valued functions grid.k). We populate array ∂ₓuh is with $i k \hat{u}$:

@. ∂ₓuh = im * grid.kr * uh

Then the derivative in physical space, ∂ₓu, is obtained with an inverse Fourier tranform. The latter is obtained again using the FFTW plans but now via ldiv!:

using LinearAlgebra: ldiv!

ldiv!(∂ₓu, grid.rfftplan, ∂ₓuh)

plot(grid.x, [u ∂ₓu], label=["u" "∂u/∂x"], xlabel="x")