Merge pull request #246 from FourierFlows/ncc/proper-dealias

Dealias `sol` before computing nonlinear terms
FourierFlows · Jun 4, 2021 · bf62e73 · bf62e73 · navidcy · Jun 4, 2021
2 parents c30c5a3 + 6e44b4d
commit bf62e73
Show file tree

Hide file tree

Showing 10 changed files with 63 additions and 5 deletions.
diff --git a/Project.toml b/Project.toml
@@ -2,7 +2,7 @@ name = "GeophysicalFlows"
 uuid = "44ee3b1c-bc02-53fa-8355-8e347616e15e"
 license = "MIT"
 authors = ["Navid C. Constantinou <navidcy@gmail.com>", "Gregory L. Wagner <wagner.greg@gmail.com>", "and co-contributors"]
-version = "0.12.3"
+version = "0.13.0"
 
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
@@ -21,7 +21,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 CUDA = "^1, ^2.4.2, ^3.0"
 DocStringExtensions = "^0.8"
 FFTW = "^1"
-FourierFlows = "^0.6.19"
+FourierFlows = "^0.7.1"
 JLD2 = "^0.1, ^0.2, ^0.3, ^0.4"
 Reexport = "^0.2, ^1"
 SpecialFunctions = "^0.10, ^1"

diff --git a/docs/Project.toml b/docs/Project.toml
@@ -8,5 +8,5 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
 [compat]
-FourierFlows = "≥ 0.6.17"
+FourierFlows = "≥ 0.7.1"
 Plots = "≥ 1.10.1"
diff --git a/docs/make.jl b/docs/make.jl
@@ -63,6 +63,7 @@ sitename = "GeophysicalFlows.jl",
    pages = Any[
             "Home"    => "index.md",
             "Installation instructions" => "installation_instructions.md",
+            "Aliasing" => "aliasing.md",
             "GPU" => "gpu.md",
             "Examples" => [
               "TwoDNavierStokes" => Any[

diff --git a/docs/src/aliasing.md b/docs/src/aliasing.md
@@ -0,0 +1,33 @@
+# Aliasing
+
+
+In pseudospectral methods, computing nonlinear terms results in aliasing errors. (Read more about
+aliasing errors in the [FourierFlows.jl Documentation](https://fourierflows.github.io/FourierFlowsDocumentation/stable/aliasing/).) To avoid aliasing errors, we need to apply some dealiasing to our fields 
+in Fourier space before transforming to physical space to compute nonlinear terms.
+
+!!! info "De-aliasing scheme"
+    FourierFlows.jl curently implements dealiasing by zeroing out the highest-`aliased_fraction` 
+    wavenumber components on a `grid`. By default in FourierFlows.jl, `aliased_fraction=1/3`.
+    Users can construct a grid with different `aliased_fraction`, e.g.,
+
+    ```julia
+    julia> grid = OneDGrid(64, 2π, aliased_fraction=1/2)
+
+    julia> 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.5
+    ```
+
+Currently, all nonlinearities in all modules included in GeophysicalFlows.jl modules are quadratic 
+nonlinearities. Therefore, the default `aliased_fraction` of 1/3 is appropriate.
+
+All modules apply de-aliasing by calling, e.g., `dealias!(prob.sol, prob.grid)` both before
+computing any nonlinear terms and also during updating all variable, i.e., within `updatevars!`.
+
+To disable de-aliasing you need to create a problem with a grid that has been constructed with 
+the keyword `aliased_fraction=0`.
diff --git a/examples/Project.toml b/examples/Project.toml
@@ -1,7 +1,6 @@
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-FourierFlows = "2aec4490-903f-5c70-9b11-9bed06a700e1"
 GeophysicalFlows = "44ee3b1c-bc02-53fa-8355-8e347616e15e"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -10,5 +9,5 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-FourierFlows = "≥ 0.6.17"
+FourierFlows = "≥ 0.7.1"
 Plots = "≥ 1.10.1"
diff --git a/src/barotropicqgql.jl b/src/barotropicqgql.jl
@@ -330,7 +330,10 @@ N = - \\widehat{𝖩(ψ, q+η)} + F̂ .
 ```
 """
 function calcN!(N, sol, t, clock, vars, params, grid)
+  dealias!(sol, grid)
+
   calcN_advection!(N, sol, t, clock, vars, params, grid)
+
   addforcing!(N, sol, t, clock, vars, params, grid)
 
   return nothing
@@ -373,6 +376,8 @@ end
 Update the `vars` of a problem `prob` that has `grid` and `params` with the solution in `sol`.
 """
 function updatevars!(sol, vars, params, grid)
+  dealias!(sol, grid)
+
   CUDA.@allowscalar sol[1, 1] = 0
   @. vars.zetah = sol
   CUDA.@allowscalar @. vars.zetah[1, :] = 0

diff --git a/src/multilayerqg.jl b/src/multilayerqg.jl
@@ -541,8 +541,12 @@ N_j = - \\widehat{𝖩(ψ_j, q_j)} - \\widehat{U_j ∂_x Q_j} - \\widehat{U_j
 function calcN!(N, sol, t, clock, vars, params, grid)
   nlayers = numberoflayers(params)
 
+  dealias!(sol, grid)
+
   calcN_advection!(N, sol, vars, params, grid)
+
   @views @. N[:, :, nlayers] += params.μ * grid.Krsq * vars.ψh[:, :, nlayers]   # bottom linear drag
+
   addforcing!(N, sol, t, clock, vars, params, grid)
 
   return nothing
@@ -688,6 +692,8 @@ end
 Update all problem variables using `sol`.
 """
 function updatevars!(vars, params, grid, sol)
+  dealias!(sol, grid)
+
   @. vars.qh = sol
   streamfunctionfrompv!(vars.ψh, vars.qh, params, grid)
   @. vars.uh = -im * grid.l  * vars.ψh

diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl
@@ -333,7 +333,10 @@ N = - \\widehat{𝖩(ψ, q+η)} + F̂ .
 ```
 """
 function calcN!(N, sol, t, clock, vars, params, grid)
+  dealias!(sol, grid)
+
   calcN_advection!(N, sol, t, clock, vars, params, grid)
+
   addforcing!(N, sol, t, clock, vars, params, grid)
 
   return nothing
@@ -389,6 +392,8 @@ end
 Update the variables in `vars` with the solution in `sol`.
 """
 function updatevars!(sol, vars, params, grid)
+  dealias!(sol, grid)
+
   @. vars.qh = sol
   streamfunctionfrompv!(vars.ψh, vars.qh, params, grid)
   @. vars.uh = -im * grid.l  * vars.ψh

diff --git a/src/surfaceqg.jl b/src/surfaceqg.jl
@@ -256,7 +256,10 @@ N = - \\widehat{𝖩(ψ, b)} + F̂ .
 ```
 """
 function calcN!(N, sol, t, clock, vars, params, grid)
+  dealias!(sol, grid)
+
   calcN_advection!(N, sol, t, clock, vars, params, grid)
+
   addforcing!(N, sol, t, clock, vars, params, grid)
 
   return nothing
@@ -301,6 +304,8 @@ Update variables in `vars` with solution in `sol`.
 function updatevars!(prob)
   vars, grid, sol = prob.vars, prob.grid, prob.sol
 
+  dealias!(sol, grid)
+
   @. vars.bh = sol
   @. vars.uh =   im * grid.l  * sqrt(grid.invKrsq) * sol
   @. vars.vh = - im * grid.kr * sqrt(grid.invKrsq) * sol

diff --git a/src/twodnavierstokes.jl b/src/twodnavierstokes.jl
@@ -274,6 +274,8 @@ N = - \\widehat{𝖩(ψ, ζ)} + F̂ .
 ```
 """
 function calcN!(N, sol, t, clock, vars, params, grid)
+  dealias!(sol, grid)
+
   calcN_advection!(N, sol, t, clock, vars, params, grid)
 
   addforcing!(N, sol, t, clock, vars, params, grid)
@@ -319,6 +321,8 @@ Update variables in `vars` with solution in `sol`.
 function updatevars!(prob)
   vars, grid, sol = prob.vars, prob.grid, prob.sol
 
+  dealias!(sol, grid)
+
   @. vars.ζh = sol
   @. vars.uh =   im * grid.l  * grid.invKrsq * sol
   @. vars.vh = - im * grid.kr * grid.invKrsq * sol