Enhance user experience with aliasing + (fix bug in getaliasedwaveumbers)

Project.toml

@@ -6,7 +6,7 @@ authors = ["Gregory L. Wagner <wagner.greg@gmail.com>", "Navid C. Constantinou <
 description = "Tools for building fast, hackable, pseudospectral partial differential equation solvers on periodic domains."
 documentation = "https://fourierflows.github.io/FourierFlowsDocumentation/stable/"
 repository = "https://github.com/FourierFlows/FourierFlows.jl"
-version = "0.6.19"
+version = "0.7.0"
 
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"

docs/make.jl

@@ -50,6 +50,7 @@ pages = [
     "Installation Instructions" => "installation_instructions.md",
     "Code Basics" => "basics.md",
     "Grids" => "grids.md",
+    "Aliasing" => "aliasing.md",
     "Problem" => "problem.md",
     "GPU" => "gpu.md",
     "Examples" => [ 

docs/src/aliasing.md

@@ -0,0 +1,77 @@
+# Aliasing
+
+
+```@setup 1
+using FourierFlows
+using Plots
+Plots.scalefontsizes(1.25)
+Plots.default(lw=3)
+```
+
+In pseudospectral methods, often aliasing errors come into play. These errors originate from
+the discrete version of the grid. A grid discretized with ``n_x`` points can only resolve a 
+total of ``n_x`` wavenumbers in Fourier space. 
+
+On a grid with a total of ``n_x`` wavenumbers, both harmonics ``e^{2\pi i k x / L_x}`` and 
+``e^{2\pi i (k+n_x) x / L_x}``, with ``k`` an integer, are indistinguishable when evaluated
+on the discrete grid-points of this grid. When we compute nonlinear terms in physical space, 
+we may end up with terms that project on higher wavenumbers, beyond those that our grid can 
+represent. In that case, those wavenumbers will be erroneously projected onto some lower 
+wavenumber mode that fits our domain.
+
+Take, for example, functions ``\cos(4x)`` and ``\cos(6x)`` and let's see how they are represented 
+on a grid ``x \in [-π, π)`` with ``n_x = 10``.
+
+```@example 1
+nx, Lx = 10, 2π
+grid = OneDGrid(nx, Lx)
+
+f1(x) = cos(4x)
+f2(x) = cos(6x)
+
+p = plot(grid.x, f1.(grid.x), lw=0, marker=:circle, c=:red, ms=8, ylims=(-1.6, 1.6), label="cos(4x)")
+plot!(p, f1, lw=3, alpha=0.2, c=:red, xlims=(-Lx/2, Lx/2), label="")
+plot!(p, grid.x, f2.(grid.x), lw=0, marker=:star5, ms=8.5, color=:blue, alpha=0.8, label="cos(6x)")
+plot!(p, f2, lw=3, alpha=0.2, c=:blue, xlims=(-Lx/2, Lx/2), label="")
+
+plot(p, xlabel="x", xlims=(-3.3, 3.3))
+
+savefig("assets/plot4.svg"); nothing # hide
+```
+
+![](assets/plot4.svg)
+
+The take home message is that on this grid we cannot distinguish harmonics of wavenumbers 4 and 6
+and attempting to represent harmonics with wavenumber 6 on this grid will lead to aliasing errors.
+For example, say that we are solving an equation on this grid and at some point we compute the product 
+``\cos(2x) \cos(4x)``. The result is ``\frac1{2} \cos(2x) + \frac1{2} \cos(6x)``, but on this 
+grid ``\cos(6x)`` is indistinguishable from ``\cos(4x)`` and, therefore, we get an answer 
+which is the sum of ``\frac1{2} \cos(2x) + \frac1{2} \cos(4x)``!
+
+There are two ways to avoid aliasing errors we either *(i)* need to discard some of the wavenumber 
+components in Fourier space before we tranform to physical space, or *(ii)* pad our Fourier 
+represresentation with more wavenumbers that will have zero power. In FourierFlows.jl the former
+is implemented
+
+!!! info "De-aliasing scheme"
+    FourierFlows.jl curently implement dealiasing by zeroing out the top-`aliased_fraction` 
+    wavenumber components on a `grid`.
+
+How many wavenumber components we need to discard depends on the order of the nonlinearity. For
+quadradic nonlinearities, one would intuitively say that we need to discard the top-1/2 of the 
+wavenumber components. However, Orszag (1972) pointed out that simply only discarding the 
+top-1/3 of wavenumber components is enough. Actally, with Orszag's so-called 2/3-rule for dealiasing, 
+still some aliasing errors occur, but only into wavenumbers that will be zero-ed out next time 
+we dealias.
+
+When constructing a `grid` we can specify the `aliased_fraction` parameter. By default, this is 
+set to ``1/3``, appropriate for quadratic nonlinearities. Then `dealias!(fh, grid)` will zero-out 
+the top-`aliased_fraction` wavenumber components of `fh`. 
+
+If we construct a grid with `aliased_fraction=0`, e.g.,
+
+```@example 1
+grid_nodealias = OneDGrid(nx, Lx; aliased_fraction=0)
+```
+
+then `dealias!(fh, grid_nodealias)` will have _no effect_ whatsoever on `fh`.

docs/src/gpu.md

@@ -26,7 +26,9 @@ OneDimensionalGrid
   ├────────── size Lx: 2.0
   ├──── resolution nx: 16
   ├── grid spacing dx: 0.125
-  └─────────── domain: x ∈ [-1.0, 0.875]
+  ├─────────── domain: x ∈ [-1.0, 0.875]
+  └─ aliased fraction: 0.3333333333333333
+
 ```
 
 gives out a grid whose arrays are `CuArrays`. (Calling `OneDGrid(n, L)` defaults to CPU, i.e., 

docs/src/grids.md

@@ -26,8 +26,9 @@ OneDimensionalGrid
   ├────────── size Lx: 6.283185307179586
   ├──── resolution nx: 64
   ├── grid spacing dx: 0.09817477042468103
-  └─────────── domain: x ∈ [-3.141592653589793, 3.0434178831651124]
-```
+  ├─────────── 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 `OneDGrid` constructor. The grid 

docs/src/problem.md

@@ -185,10 +185,10 @@ plot!(x -> cos(π * x) * exp(-prob.params.α * 2), -1, 1, label = "analytical")
 
 plot!(x -> cos(π * x), -1, 1, linestyle=:dash, color=:gray, label = "initial condition")
 
-savefig("assets/plot4.svg"); nothing # hide
+savefig("assets/plot5.svg"); nothing # hide
 ```
 
-![](assets/plot4.svg)
+![](assets/plot5.svg)
 
 A good practice is to encompass all functions and type definitions related with a PDE under 
 a single module, e.g.,

examples/OneDShallowWaterGeostrophicAdjustment.jl

@@ -120,7 +120,7 @@ function calcN!(N, sol, t, clock, vars, params, grid)
   @. N[:, 2] = - params.f * vars.uh                                        #  - f u
   @. N[:, 3] = - im * grid.kr * params.H * vars.uh                         #  - H ∂u/∂x
 
-  dealias!(N, grid, grid.kralias)
+  dealias!(N, grid)
 
   return nothing
 end

src/FourierFlows.jl

@@ -77,7 +77,7 @@ using Base: fieldnames
 using FFTW: fftfreq, rfftfreq
 
 "Abstract supertype for grids."
-abstract type AbstractGrid{T, A} end
+abstract type AbstractGrid{T, A, Alias} end
 
 "Abstract supertype for timesteppers."
 abstract type AbstractTimeStepper{T} end