# 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()
```

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()
```

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()
```