Add Timestepping sections in Docs + info about filtering

.github/workflows/Documenter.yml

@@ -11,10 +11,10 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@latest
         with:
-          version: 1.6
+          version: 1.9
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy

docs/Project.toml

@@ -2,7 +2,6 @@
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

docs/make.jl

@@ -2,7 +2,6 @@ using
   Documenter,
   Literate,
   CairoMakie,  # so that Literate.jl does not capture precompilation output
-  Glob,
   FourierFlows
 
 #####
@@ -29,14 +28,10 @@ end
 ##### Build and deploy docs
 #####
 
-# Set up a timer to print a space ' ' every 240 seconds. This is to avoid CI
-# timing out when building demanding Literate.jl examples.
-Timer(t -> println(" "), 0, interval=240)
-
 format = Documenter.HTML(
     collapselevel = 2,
        prettyurls = get(ENV, "CI", nothing) == "true",
-        canonical = "https://fourierflows.github.io/FourierFlowsDocumentation/dev/",
+        canonical = "https://fourierflows.github.io/FourierFlowsDocumentation/stable",
 )
 
 pages = [
@@ -45,6 +40,7 @@ pages = [
     "Code Basics" => "basics.md",
     "Grids" => "grids.md",
     "Aliasing" => "aliasing.md",
+    "Time stepping" => "timestepping.md",
     "Problem" => "problem.md",
     "Diagnostics" => "diagnostics.md",
     "Output" => "output.md",
@@ -73,17 +69,33 @@ makedocs(;
   checkdocs = :exports
 )
 
-@info "Cleaning up temporary .jld2 and .nc files created by doctests..."
-for file in vcat(glob("docs/*.jld2"))
+@info "Cleaning up temporary .jld2 and .nc files created by doctests or literated examples..."
+
+"""
+    recursive_find(directory, pattern)
+
+Return list of filepaths within `directory` that contains the `pattern::Regex`.
+"""
+recursive_find(directory, pattern) =
+    mapreduce(vcat, walkdir(directory)) do (root, dirs, files)
+        joinpath.(root, filter(contains(pattern), files))
+    end
+
+files = []
+for pattern in [r"\.jld2", r"\.nc"]
+  global files = vcat(files, recursive_find(@__DIR__, pattern))
+end
+
+for file in files
   rm(file)
 end
 
 withenv("GITHUB_REPOSITORY" => "FourierFlows/FourierFlowsDocumentation") do
   deploydocs(
-            repo = "github.com/FourierFlows/FourierFlowsDocumentation.git",
-        versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"],
-       forcepush = true,
-    push_preview = false,
-       devbranch = "main"
+             repo = "github.com/FourierFlows/FourierFlowsDocumentation.git",
+         versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"],
+        forcepush = true,
+     push_preview = true,
+        devbranch = "main"
   )
 end

docs/src/aliasing.md

@@ -1,4 +1,4 @@
-# Aliasing
+# [Aliasing](@id aliasing)
 
 ```@setup 1
 using FourierFlows
@@ -96,7 +96,7 @@ the highest-`aliased_fraction` wavenumber components of `fh`.
 FourierFlows.dealias!
 ```
 
-If we construct a grid with `aliased_fraction=0`, e.g.,
+If we construct a grid with `aliased_fraction = 0`, e.g.,
 
 ```@example 1
 grid_nodealias = OneDGrid(; nx, Lx, aliased_fraction=0)

docs/src/problem.md

@@ -110,20 +110,18 @@ equation = FourierFlows.Equation(L, calcN!, grid)
 ```
 
 Last, we have to pick a time-stepper and a time-step `dt` and gather everything 
-a FourierFlows's [`Problem`](@ref FourierFlows.Problem):
+a FourierFlows's [`Problem`](@ref FourierFlows.Problem). Time-steppers are prescribed
+via a string that corresponds to the name of the implemented time-steppers _without_
+the `TimeStepper` ending. (See [Time-stepping section](@ref timestepping) for a list
+of implemented time-stepping schemes.) For example, a problem that uses a
+[`ForwardEulerTimeStepper`](@ref) with time step of 0.02 is constructed via:
 
 ```@example 2
 stepper, dt = "ForwardEuler", 0.02
 
 prob = FourierFlows.Problem(equation, stepper, dt, grid, vars, params)
 ```
 
-Currently, the implemented time-steppers are [`ForwardEuler`](@ref ForwardEulerTimeStepper),
-[`AB3`](@ref AB3TimeStepper) (Adams-Basmforth 3rd order), [`RK4`](@ref RK4TimeStepper) (Runge-Kutta
-4th order), and [`ETDRK4`](@ref ETDRK4TimeStepper) (Exponential Time Differencing Runge-Kutta
-4th order). Also, there exist the `Filtered` versions of all the above, in which a high-wavenumber
-filter is applied after every time-step.
-
 By default, the `Problem` constructor takes `sol` a complex valued array filled with zeros
 with same size as `L`.
 

docs/src/timestepping.md

@@ -0,0 +1,64 @@
+# [Time-stepping](@id timestepping)
+
+FourierFlows.jl includes several time-stepping algorithms.
+
+Most of the time-stepping algorithms are fully explicit schemes: [`ForwardEulerTimeStepper`](@ref), [`AB3TimeStepper`](@ref), [`RK4TimeStepper`](@ref), and [`LSRK54TimeStepper`](@ref)
+but also implemented is the [`ETDRK4TimeStepper`](@ref).
+
+## High-wavenumber filtering
+
+Most of the time steppers also come with their `Filtered` equivalents: [`FilteredForwardEulerTimeStepper`](@ref), [`FilteredAB3TimeStepper`](@ref), [`FilteredRK4TimeStepper`](@ref), [`FilteredLSRK54TimeStepper`](@ref), and [`FilteredETDRK4TimeStepper`](@ref).
+
+The filtered time steppers apply a high-wavenumber filter to the solution at the end of each step.
+The motivation behind filtering is to remove enstrophy accumulating at high wavenumbers and creating 
+noise at grid-scale level.
+
+The high-wavenumber filter used in the filtered timesteppers is:
+
+```math
+\mathrm{filter}(\boldsymbol{k}) = 
+     \begin{cases}
+       1 & \quad |\boldsymbol{k}| ≤ k_{\textrm{cutoff}} \, ,\\ 
+       \exp{ \left [- \alpha (|\boldsymbol{k}| - k_{\textrm{cutoff}})^p \right]} & \quad |\boldsymbol{k}| > k_{\textrm{cutoff}} \, .
+     \end{cases}
+```
+
+For fluid equations with quadratic non-linearities it makes sense to choose a cutoff wavenumber
+at 2/3 of the highest wavenumber resolved in our domain, ``k_{\textrm{cutoff}} = \tfrac{2}{3} k_{\textrm{max}}`` (see discussion in [Aliasing section](@ref aliasing)).
+
+Given the order ``p``, we calculate the coefficient ``\alpha`` so that the the filter value
+that corresponds to the highest allowed wavenumber in our domain is a small value, ``\delta``,
+taken to be close to machine precision.
+
+That is:
+
+```math
+\alpha = \frac{- \log\delta}{(k_{\textrm{max}} - k_{\textrm{cutoff}})^p} \ .
+```
+
+The above filter originates from Canuto et al. (1988). In geophysical turbulence applications
+it was used by LaCasce (1996) and later by Arbic & Flierl (2003).
+
+Using the default parameters provided by the filtered time steppers (see
+[`FourierFlows.makefilter`](@ref)), the filter has the following form:
+
+```@setup 1
+using CairoMakie
+CairoMakie.activate!(type = "svg")
+set_theme!(Theme(linewidth = 3, fontsize = 20))
+```
+
+```@example 1
+using FourierFlows, CairoMakie
+
+K = 0:0.01:1 # non-dimensional wavenumber k * dx / π
+
+filter = FourierFlows.makefilter(K)
+
+fig = Figure()
+ax = Axis(fig[1, 1], xlabel = "k dx / π", ylabel = "filter", aspect=2.5, xticks=0:0.2:1)
+
+lines!(ax, K, filter)
+
+current_figure() # hide
+```

src/domains.jl

@@ -401,7 +401,7 @@
     getaliasedwavenumbers(nk, nkr, aliased_fraction)
 
 Return the top `aliased_fraction` highest wavenumbers, both for and real FFTs, `kalias` and 
-`kralias` respectively. For example, `aliased_fraction=1/3` should return the indices of the 
+`kralias` respectively. For example, `aliased_fraction = 1/3` should return the indices of the 
 top-most 1/6-th (in absolute value) for both positive and negative wavenumbers (i.e., 1/3 total) 
 that should be set to zero after performing an FFT. 
 """
@@ -476,40 +476,63 @@
 end
 
 """
-    makefilter(K; order=4, innerK=0.65, outerK=1)
+    makefilter(K; order=4, innerK=2/3, outerK=1, tol=1e-15)
 
-Return a filter acting on the non-dimensional wavenumber `K` that decays exponentially
-for `K > innerK`, thus removing high-wavenumber content from a spectrum it is multiplied
-with. The decay rate is determined by order and `outerK` determines the outer wavenumber
-at which the filter is smaller than `Float64` machine precision.
+Return a filter acting on the non-dimensional wavenumber `K`. For a
+one-dimensional grid, the non-dimensional wavenumber `K` is
+
+```julia
+K = k * dx / π
+```
+
+and thus it takes values `K` ``∈ [-1, 1]``.
+
+For `K ≤ innerK` the filter is inactive, i.e., equal to 1. For `K > innerK`,
+the filter decays exponentially to remove high-wavenumber content from 
+the solution, i.e.,
+
+```julia
+filter(K) = exp(- decay * (K - innerK)^order)
+```
+
+For a given `order` and , the `decay` rate is determined so that the filter value at the
+outer wavenumber `outerK` is `tol`, where `tol` is a small number, close to machine precision.
+
+```julia
+decay = - log(tol) / (outerK - innerK)^order
+```
 """
-function makefilter(K::Array; order=4, innerK=0.65, outerK=1)
+function makefilter(K::Array; order=4, innerK=2/3, outerK=1, tol=1e-15)
   TK = typeof(K)
   K = Array(K)
-  decay = 15*log(10) / (outerK - innerK)^order # decay rate for filtering function
-  filter = @. exp(-decay*(K - innerK)^order)
+
+  decay = -log(tol) / (outerK - innerK)^order
+
+  filter = @. exp(- decay * (K - innerK)^order)
   filter[K .< innerK] .= 1
-  
+
   return TK(filter)
 end
 
+makefilter(K::AbstractRange; kwargs...) = makefilter(Array(K); kwargs...)
+
 function makefilter(g::OneDGrid; realvars=true, kwargs...)
-  K = realvars ? g.kr*g.dx/π : @.(abs(g.k*g.dx/π))
-  
+  K = realvars ? g.kr * g.dx / π : @.(abs(g.k * g.dx / π))
+
   return makefilter(K; kwargs...)
 end
 
 function makefilter(g::TwoDGrid; realvars=true, kwargs...)
-  K = realvars ?
-      @.(sqrt((g.kr*g.dx/π)^2 + (g.l*g.dy/π)^2)) : @.(sqrt((g.k*g.dx/π)^2 + (g.l*g.dy/π)^2))
-  
+  K = realvars ? @.(sqrt((g.kr * g.dx / π)^2 + (g.l * g.dy / π)^2)) :
+                 @.(sqrt((g.k  * g.dx / π)^2 + (g.l * g.dy / π)^2))
+
   return makefilter(K; kwargs...)
 end
 
 function makefilter(g::ThreeDGrid; realvars=true, kwargs...)
-  K = realvars ?
-      @.(sqrt((g.kr*g.dx/π)^2 + (g.l*g.dy/π)^2 + (g.m*g.dz/π)^2)) : @.(sqrt((g.k*g.dx/π)^2 + (g.l*g.dy/π)^2 + (g.m*g.dz/π)^2))
-  
+  K = realvars ? @.(sqrt((g.kr * g.dx / π)^2 + (g.l * g.dy / π)^2 + (g.m * g.dz / π)^2)) :
+                 @.(sqrt((g.k  * g.dx / π)^2 + (g.l * g.dy / π)^2 + (g.m * g.dz / π)^2))
+
   return makefilter(K; kwargs...)
 end
 

src/timesteppers.jl

@@ -50,7 +50,8 @@ isexplicit(stepper) = any(Symbol(stepper) .== fullyexplicitsteppers)
     TimeStepper(stepper, equation, dt=nothing, dev=CPU(); kw...)
 
 Instantiate the Time`stepper` for `equation` with timestep `dt` and
-on the `dev`ice. The `kw` are passed to the timestepper constructor.
+on the `dev`ice. The keyword arguments, `kw`, are passed to the
+timestepper constructor.
 """
 function TimeStepper(stepper, equation, dt=nothing, dev::Device=CPU(); kw...)
   fullsteppername = Symbol(stepper, :TimeStepper)