Dealias sol before computing nonlinear terms

Project.toml

@@ -2,7 +2,7 @@ name = "GeophysicalFlows"
 uuid = "44ee3b1c-bc02-53fa-8355-8e347616e15e"
 license = "MIT"
 authors = ["Navid C. Constantinou <[email protected]>", "Gregory L. Wagner <[email protected]>", "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"

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"

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[

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`.

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"

src/barotropicqgql.jl

@@ -297,7 +297,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
@@ -340,6 +343,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

src/multilayerqg.jl

@@ -496,8 +496,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
@@ -643,6 +647,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

src/singlelayerqg.jl

@@ -299,7 +299,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
@@ -355,6 +358,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

src/surfaceqg.jl

@@ -230,7 +230,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
@@ -275,6 +278,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

src/twodnavierstokes.jl

@@ -244,6 +244,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)
@@ -289,6 +291,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