GitHub

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

using FourierFlows

nx, Lx = 64, 2π

grid = OneDGrid(; nx, Lx, x0 = -Lx/3)

# output

OneDimensionalGrid
  ├─────────── Device: CPU
  ├──────── FloatType: Float64
  ├────────── size Lx: 6.283185307179586
  ├──── resolution nx: 64
  ├── grid spacing dx: 0.09817477042468103
  ├─────────── domain: x ∈ [-2.0943951023931953, 4.090615434361711]
  └─ 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 redundant.

We can define an array u that contains the values of a function $u(x)$ on this grid. If, for example, $u$ is

\[u(x) = \sin(2x) + \frac{1}{2} \cos(5x) ,\]

then

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

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

using CairoMakie

lines(grid.x, u, label="u", axis = (xlabel="x",))
axislegend()
Example block output

Any periodic 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 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 ComplexF64
(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_avx2_128"
    (rdft2-r2hc-direct-2 "r2cf_2")
    (rdft2-r2hc01-direct-2 "r2cfII_2"))
  (dft-direct-32 "n2fv_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

fig = Figure()
ax = Axis(fig[1, 1];
          xlabel = "k",
          limits = ((-0.5, 10.5), (-0.55, 0.55)),
          xticks = 0:10,
          yticks = -0.5:0.25:0.5)

scatterlines!(ax, grid.kr, real.(uhat);
              markersize = 16, label = "real(û)")

scatterlines!(ax, grid.kr, imag.(uhat);
              markersize = 22, marker = :diamond, label = "imag(û)")

axislegend()
Example block output

Note that with the normalization we get what we expect since

\[u(x) = \sin(2x) + \frac1{2} \cos(5x) = -\frac{i}{2} (e^{2ix} - e^{-2ix}) + \frac{1}{4} (e^{5ix} + e^{-5ix}) .\]

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

∂ₓ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 the array ∂ₓuh with $\mathrm{i} k \hat{u}$:

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

Then the derivative in physical space, ∂ₓu, is obtained through an inverse Fourier transform. For this again we use the FFTW plans but now via ldiv!:

using LinearAlgebra: ldiv!

ldiv!(∂ₓu, grid.rfftplan, ∂ₓuh)
64-element Vector{Float64}:
 -3.1650635094610906
 -2.6388612568260568
 -1.750214503657824
 -0.7619779085543852
  0.05153169814995007
  0.4708904837349843
  0.3844356009901506
 -0.18487950393303665
 -1.0850031948125751
 -2.0713702725438923
  ⋮
  3.718520440734795
  4.146865373291539
  3.9383720631513
  3.0971223767148777
  1.7652025002410927
  0.19132260965489634
 -1.3253879837092064
 -2.5032342382718373
 -3.1415403992972397

We can plot ∂ₓu and compare it with the analytic derivative. 

fig = Figure()
ax = Axis(fig[1, 1], xlabel = "x")

lines!(ax, grid.x, u;
       label = "u")
lines!(ax, grid.x, ∂ₓu;
       label = "∂u/∂x")
lines!(ax, grid.x, @. 2 * cos(2 * grid.x) - 5/2 * sin(5 * grid.x);
       linewidth=1.5, linestyle=:dash, label = "analytical ∂u/∂x")

axislegend()
Example block output