Skip to content

Commit 7a963c1

Browse files
authored
Merge pull request #285 from FourierFlows/ncc/aliased_fraction-part-of-grid
Enhance user experience with aliasing + (fix bug in `getaliasedwaveumbers`)
2 parents a4bcb46 + 2939929 commit 7a963c1

File tree

11 files changed

+415
-182
lines changed

11 files changed

+415
-182
lines changed

Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ authors = ["Gregory L. Wagner <wagner.greg@gmail.com>", "Navid C. Constantinou <
66
description = "Tools for building fast, hackable, pseudospectral partial differential equation solvers on periodic domains."
77
documentation = "https://fourierflows.github.io/FourierFlowsDocumentation/stable/"
88
repository = "https://github.com/FourierFlows/FourierFlows.jl"
9-
version = "0.6.19"
9+
version = "0.7.0"
1010

1111
[deps]
1212
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"

docs/make.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,7 @@ pages = [
5050
"Installation Instructions" => "installation_instructions.md",
5151
"Code Basics" => "basics.md",
5252
"Grids" => "grids.md",
53+
"Aliasing" => "aliasing.md",
5354
"Problem" => "problem.md",
5455
"GPU" => "gpu.md",
5556
"Examples" => [

docs/src/aliasing.md

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
# Aliasing
2+
3+
4+
```@setup 1
5+
using FourierFlows
6+
using Plots
7+
Plots.scalefontsizes(1.25)
8+
Plots.default(lw=3)
9+
```
10+
11+
In pseudospectral methods, often aliasing errors come into play. These errors originate from
12+
the discrete version of the grid. A grid discretized with ``n_x`` points can only resolve a
13+
total of ``n_x`` wavenumbers in Fourier space.
14+
15+
On a grid with a total of ``n_x`` wavenumbers, both harmonics ``e^{2\pi i k x / L_x}`` and
16+
``e^{2\pi i (k+n_x) x / L_x}``, with ``k`` an integer, are indistinguishable when evaluated
17+
on the discrete grid-points of this grid. When we compute nonlinear terms in physical space,
18+
we may end up with terms that project on higher wavenumbers, beyond those that our grid can
19+
represent. In that case, those wavenumbers will be erroneously projected onto some lower
20+
wavenumber mode that fits our domain.
21+
22+
Take, for example, functions ``\cos(4x)`` and ``\cos(6x)`` and let's see how they are represented
23+
on a grid ``x \in [-π, π)`` with ``n_x = 10``.
24+
25+
```@example 1
26+
nx, Lx = 10, 2π
27+
grid = OneDGrid(nx, Lx)
28+
29+
f1(x) = cos(4x)
30+
f2(x) = cos(6x)
31+
32+
p = plot(grid.x, f1.(grid.x), lw=0, marker=:circle, c=:red, ms=8, ylims=(-1.6, 1.6), label="cos(4x)")
33+
plot!(p, f1, lw=3, alpha=0.2, c=:red, xlims=(-Lx/2, Lx/2), label="")
34+
plot!(p, grid.x, f2.(grid.x), lw=0, marker=:star5, ms=8.5, color=:blue, alpha=0.8, label="cos(6x)")
35+
plot!(p, f2, lw=3, alpha=0.2, c=:blue, xlims=(-Lx/2, Lx/2), label="")
36+
37+
plot(p, xlabel="x", xlims=(-3.3, 3.3))
38+
39+
savefig("assets/plot4.svg"); nothing # hide
40+
```
41+
42+
![](assets/plot4.svg)
43+
44+
The take home message is that on this grid we cannot distinguish harmonics of wavenumbers 4 and 6
45+
and attempting to represent harmonics with wavenumber 6 on this grid will lead to aliasing errors.
46+
For example, say that we are solving an equation on this grid and at some point we compute the product
47+
``\cos(2x) \cos(4x)``. The result is ``\frac1{2} \cos(2x) + \frac1{2} \cos(6x)``, but on this
48+
grid ``\cos(6x)`` is indistinguishable from ``\cos(4x)`` and, therefore, we get an answer
49+
which is the sum of ``\frac1{2} \cos(2x) + \frac1{2} \cos(4x)``!
50+
51+
There are two ways to avoid aliasing errors we either *(i)* need to discard some of the wavenumber
52+
components in Fourier space before we tranform to physical space, or *(ii)* pad our Fourier
53+
represresentation with more wavenumbers that will have zero power. In FourierFlows.jl the former
54+
is implemented
55+
56+
!!! info "De-aliasing scheme"
57+
FourierFlows.jl curently implement dealiasing by zeroing out the top-`aliased_fraction`
58+
wavenumber components on a `grid`.
59+
60+
How many wavenumber components we need to discard depends on the order of the nonlinearity. For
61+
quadradic nonlinearities, one would intuitively say that we need to discard the top-1/2 of the
62+
wavenumber components. However, Orszag (1972) pointed out that simply only discarding the
63+
top-1/3 of wavenumber components is enough. Actally, with Orszag's so-called 2/3-rule for dealiasing,
64+
still some aliasing errors occur, but only into wavenumbers that will be zero-ed out next time
65+
we dealias.
66+
67+
When constructing a `grid` we can specify the `aliased_fraction` parameter. By default, this is
68+
set to ``1/3``, appropriate for quadratic nonlinearities. Then `dealias!(fh, grid)` will zero-out
69+
the top-`aliased_fraction` wavenumber components of `fh`.
70+
71+
If we construct a grid with `aliased_fraction=0`, e.g.,
72+
73+
```@example 1
74+
grid_nodealias = OneDGrid(nx, Lx; aliased_fraction=0)
75+
```
76+
77+
then `dealias!(fh, grid_nodealias)` will have _no effect_ whatsoever on `fh`.

docs/src/gpu.md

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,9 @@ OneDimensionalGrid
2626
├────────── size Lx: 2.0
2727
├──── resolution nx: 16
2828
├── grid spacing dx: 0.125
29-
└─────────── domain: x [-1.0, 0.875]
29+
├─────────── domain: x [-1.0, 0.875]
30+
└─ aliased fraction: 0.3333333333333333
31+
3032
```
3133

3234
gives out a grid whose arrays are `CuArrays`. (Calling `OneDGrid(n, L)` defaults to CPU, i.e.,

docs/src/grids.md

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -26,8 +26,9 @@ OneDimensionalGrid
2626
├────────── size Lx: 6.283185307179586
2727
├──── resolution nx: 64
2828
├── grid spacing dx: 0.09817477042468103
29-
└─────────── domain: x ∈ [-3.141592653589793, 3.0434178831651124]
30-
```
29+
├─────────── domain: x ∈ [-3.141592653589793, 3.0434178831651124]
30+
└─ aliased fraction: 0.3333333333333333
31+
```
3132

3233
The grid domain is, by default, constructed symmetrically around ``x = 0``, but this
3334
can be altered using the `x0` keyword argument of `OneDGrid` constructor. The grid

docs/src/problem.md

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -185,10 +185,10 @@ plot!(x -> cos(π * x) * exp(-prob.params.α * 2), -1, 1, label = "analytical")
185185
186186
plot!(x -> cos(π * x), -1, 1, linestyle=:dash, color=:gray, label = "initial condition")
187187
188-
savefig("assets/plot4.svg"); nothing # hide
188+
savefig("assets/plot5.svg"); nothing # hide
189189
```
190190

191-
![](assets/plot4.svg)
191+
![](assets/plot5.svg)
192192

193193
A good practice is to encompass all functions and type definitions related with a PDE under
194194
a single module, e.g.,

examples/OneDShallowWaterGeostrophicAdjustment.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -120,7 +120,7 @@ function calcN!(N, sol, t, clock, vars, params, grid)
120120
@. N[:, 2] = - params.f * vars.uh # - f u
121121
@. N[:, 3] = - im * grid.kr * params.H * vars.uh # - H ∂u/∂x
122122

123-
dealias!(N, grid, grid.kralias)
123+
dealias!(N, grid)
124124

125125
return nothing
126126
end

src/FourierFlows.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,7 @@ using Base: fieldnames
7777
using FFTW: fftfreq, rfftfreq
7878

7979
"Abstract supertype for grids."
80-
abstract type AbstractGrid{T, A} end
80+
abstract type AbstractGrid{T, A, Alias} end
8181

8282
"Abstract supertype for timesteppers."
8383
abstract type AbstractTimeStepper{T} end

0 commit comments

Comments
 (0)