Merge pull request #351 from FourierFlows/ncc/docstrings

Better phrasing and formatting for some docstrings

src/timesteppers.jl

@@ -1,13 +1,13 @@
 """
-    stepforward!(prob)
+    stepforward!(prob::Problem)
 
 Step forward `prob` one time step.
 """
 stepforward!(prob::Problem) =
   stepforward!(prob.sol, prob.clock, prob.timestepper, prob.eqn, prob.vars, prob.params, prob.grid)
 
 """
-    stepforward!(prob, nsteps::Int)
+    stepforward!(prob::Problem, nsteps::Int)
 
 Step forward `prob` for `nsteps`.
 """
@@ -87,7 +87,7 @@ end
     struct ForwardEulerTimeStepper{T} <: AbstractTimeStepper{T}
 
 A Forward Euler timestepper for time-stepping `∂u/∂t = RHS(u, t)` via:
-```
+```julia
 uⁿ⁺¹ = uⁿ + dt * RHS(uⁿ, tⁿ)
 ```
 """
@@ -154,21 +154,21 @@ end
     struct RK4TimeStepper{T} <: AbstractTimeStepper{T}
 
 A 4th-order Runge-Kutta timestepper for time-stepping `∂u/∂t = RHS(u, t)` via:
-```
+```julia
 uⁿ⁺¹ = uⁿ + dt/6 * (k₁ + 2 * k₂ + 2 * k₃ + k₄)
 ```
 where
-```
+```julia
 k₁ = RHS(uⁿ, tⁿ)
 k₂ = RHS(uⁿ + k₁ * dt/2, tⁿ + dt/2)
 k₃ = RHS(uⁿ + k₂ * dt/2, tⁿ + dt/2)
 k₄ = RHS(uⁿ + k₃ * dt, tⁿ + dt)
 ```
 
 !!! info "Usage"
-    If you are limited by memory then consider switching to [`LSRK54TimeStepper`](@ref).
-    The [`LSRK54TimeStepper`](@ref) timestepper has half the memory footprint compared
-    to the `RK4TimeStepper` with a 25~30% performance trade off.
+    If your simulation is limited by memory then consider switching to [`LSRK54TimeStepper`](@ref).
+    The [`LSRK54TimeStepper`](@ref) timestepper has about half the memory footprint compared
+    to the `RK4TimeStepper` with a about 25%-30% performance trade off.
 """
 struct RK4TimeStepper{T} <: AbstractTimeStepper{T}
   sol₁ :: T
@@ -286,7 +286,7 @@ end
 
 A 4th-order 5-stages 2-storage Runge-Kutta timestepper for time-stepping
 `∂u/∂t = RHS(u, t)` via:
-```
+```julia
 S² = 0
 
 for i = 1:5
@@ -298,16 +298,16 @@ uⁿ⁺¹ = uⁿ
 ```
 
 where `Aᵢ`, `Bᵢ`, and `Cᵢ` are the ``A``, ``B``, and ``C`` coefficients from
-the LSRK tableau table at the ``i``-th stage. For details, please refer to
+the LSRK table at the ``i``-th stage. For details, please refer to
 
 > Carpenter, M. H. and Kennedy, C. A. (1994). Fourth-order 2N-storage Runge–Kutta schemes, Technical Report NASA TM-109112, NASA Langley Research Center, VA.
 
 !!! info "Usage"
     The `LSRK54TimeStepper` is *slower* than the [`RK4TimeStepper`](@ref) but
-    with *less* memory footprint; half compared to [`RK4TimeStepper`](@ref).
+    has *less* memory footprint; half compared to [`RK4TimeStepper`](@ref).
     
-    If you are bound by performance then use [`RK4TimeStepper`](@ref); if your
-    simulation is bound by memory then consider using `LSRK54TimeStepper`.
+    If your simulation is bound by performance then use [`RK4TimeStepper`](@ref);
+    if your simulation is bound by memory then consider using `LSRK54TimeStepper`.
 """
 struct LSRK54TimeStepper{T,V} <: AbstractTimeStepper{T}
    S² :: T
@@ -320,7 +320,7 @@ end
 """
     LSRK54TimeStepper(equation::Equation, dev::Device=CPU())
 
-Construct a 4th-order 5-stages low storage Runge-Kutta timestepper for `equation` on device `dev`.
+Construct a 4th-order 5-stages low-storage Runge-Kutta timestepper for `equation` on device `dev`.
 """
 function LSRK54TimeStepper(equation::Equation, dev::Device=CPU())
   @devzeros typeof(dev) equation.T equation.dims S² RHS
@@ -379,7 +379,7 @@ end
 A 4th-order exponential-time-differencing Runge-Kutta timestepper for time-stepping
 `∂u/∂t = L * u + N(u)`. The scheme treats the linear term `L` exact while for the
 nonlinear terms `N(u)` it uses a 4th-order Runge-Kutta scheme. That is,
-```
+```julia
 uⁿ⁺¹ = exp(L * dt) * uⁿ + RK4(N(uⁿ))
 ```
 For more info refer to 
@@ -528,14 +528,14 @@ const ab3h3 = 5/12
     struct AB3TimeStepper{T} <: AbstractTimeStepper{T}
 
 A 3rd-order Adams-Bashforth timestepper for time-stepping `∂u/∂t = RHS(u, t)` via:
-```
+```julia
 uⁿ⁺¹ = uⁿ + dt/12 * (23 * RHS(uⁿ, tⁿ) - 16 * RHS(uⁿ⁻¹, tⁿ⁻¹) + 5 * RHS(uⁿ⁻², tⁿ⁻²))
 ```
 
 Adams-Bashforth is a multistep method, i.e., it not only requires information from the `n`-th time-step
 (`uⁿ`) but also from the previous two timesteps (`uⁿ⁻¹` and `uⁿ⁻²`). For the first two timesteps, it
 falls back to a forward Euler timestepping scheme:
-```
+```julia
 uⁿ⁺¹ = uⁿ + dt * RHS(uⁿ, tⁿ)
 ```
 """

src/utils.jl

@@ -97,82 +97,89 @@ end
 """
     parsevalsum2(uh, grid)
 
-Return `Σ |uh|²` on the `grid`, which is equal to the domain integral `u` squared.
-For example on a 2D grid, `parsevalsum2` returns
+Return the sum of `|uh|²` on the `grid`, which is equal to the domain integral of
+`u²`. For example on a 2D grid, `parsevalsum2` returns
 
 ```math
-\\sum_{𝐤} |û_{𝐤}|² L_x L_y = \\iint u² \\, 𝖽x 𝖽y \\,,
+\\sum_{𝐤} |û_{𝐤}|² L_x L_y = \\iint u² 𝖽x 𝖽y ,
 ```
 
 where ``û_{𝐤} =`` `uh` `` / (n_x e^{i 𝐤 ⋅ 𝐱₀})``. The elements of the vector ``𝐱₀`` are the
 left-most position in each direction, e.g., for a 2D grid `(grid.x[1], grid.y[1])`.
+
+When the input `uh` comes from a real-FFT transform, `parsevalsum2` takes care to
+count the contribution from certain ``k``-modes twice.
 """
 function parsevalsum2(uh, grid::TwoDGrid)
   if size(uh, 1) == grid.nkr  # uh is in conjugate symmetric form
-    U = sum(abs2, uh[1, :])                  # k=0 modes
-    U += sum(abs2, uh[grid.nkr, :])          # k=nx/2 modes
-    U += 2 * sum(abs2, uh[2:grid.nkr-1, :])  # sum twice for 0 < k < nx/2 modes
+    Σ = sum(abs2, uh[1, :])                  # k = 0 modes
+    Σ += sum(abs2, uh[grid.nkr, :])          # k = nx/2 modes
+    Σ += 2 * sum(abs2, uh[2:grid.nkr-1, :])  # sum twice for 0 < k < nx/2 modes
   else # count every mode once
-    U = sum(abs2, uh)
+    Σ = sum(abs2, uh)
   end
 
-  norm = grid.Lx * grid.Ly / (grid.nx^2 * grid.ny^2) # normalization for dft
+  Σ *= grid.Lx * grid.Ly / (grid.nx^2 * grid.ny^2) # normalization for dft
 
-  return norm * U
+  return Σ
 end
 
 function parsevalsum2(uh, grid::OneDGrid)
   if size(uh, 1) == grid.nkr  # uh is conjugate symmetric
-    U = sum(abs2, CUDA.@allowscalar uh[1])          # k=0 mode
-    U += sum(abs2, CUDA.@allowscalar uh[grid.nkr])  # k=nx/2 mode
-    U += @views 2 * sum(abs2, uh[2:grid.nkr-1])     # sum twice for 0 < k < nx/2 modes
+    Σ = sum(abs2, CUDA.@allowscalar uh[1])          # k = 0 mode
+    Σ += sum(abs2, CUDA.@allowscalar uh[grid.nkr])  # k = nx/2 mode
+    Σ += @views 2 * sum(abs2, uh[2:grid.nkr-1])     # sum twice for 0 < k < nx/2 modes
   else # count every mode once
-    U = sum(abs2, uh)
+    Σ = sum(abs2, uh)
   end
 
-  norm = grid.Lx / grid.nx^2 # normalization for dft
+  Σ *= grid.Lx / grid.nx^2 # normalization for dft
 
-  return norm * U
+  return Σ
 end
 
 """
     parsevalsum(uh, grid)
 
-Return `real(Σ uh)` on the `grid`, i.e.,
+Return the real part of the sum of `uh` on the `grid`. For example on a 2D grid,
+`parsevalsum` returns
 
 ```math
-ℜ [ \\sum_{𝐤} û_{𝐤} L_x L_y ] \\,,
+ℜ [ \\sum_{𝐤} û_{𝐤} L_x L_y ] ,
 
 ```
 where ``û_{𝐤} =`` `uh` `` / (n_x e^{i 𝐤 ⋅ 𝐱₀})``. The elements of the vector ``𝐱₀`` are the
 left-most position in each direction, e.g., for a 2D grid `(grid.x[1], grid.y[1])`.
+
+When the input `uh` comes from a real-FFT transform, `parsevalsum` takes care to
+count the contribution from certain ``k``-modes twice.
 """
 function parsevalsum(uh, grid::TwoDGrid)
   if size(uh, 1) == grid.nkr  # uh is conjugate symmetric
-    U = sum(uh[1, :])                  # k = 0 modes
-    U += sum(uh[grid.nkr, :])          # k = nx/2 modes
-    U += 2 * sum(uh[2:grid.nkr-1, :])  # sum twice for 0 < k < nx/2 modes
+    Σ = sum(uh[1, :])                  # k = 0 modes
+    Σ += sum(uh[grid.nkr, :])          # k = nx/2 modes
+    Σ += 2 * sum(uh[2:grid.nkr-1, :])  # sum twice for 0 < k < nx/2 modes
   else # count every mode once
-    U = sum(uh)
+    Σ = sum(uh)
   end
 
-  norm = grid.Lx * grid.Ly / (grid.nx^2 * grid.ny^2) # normalization for dft
+  Σ *= grid.Lx * grid.Ly / (grid.nx^2 * grid.ny^2) # normalization for dft
 
-  return norm * real(U)
+  return real(Σ)
 end
 
 function parsevalsum(uh, grid::OneDGrid)
   if size(uh, 1) == grid.nkr        # uh is conjugate symmetric
-    U = uh[1]                       # k=0 mode
-    U += uh[grid.nkr]               # k=nx/2 mode
-    U += 2 * sum(uh[2:grid.nkr-1])  # sum twice for 0 < k < nx/2 modes
+    Σ = uh[1]                       # k = 0 mode
+    Σ += uh[grid.nkr]               # k = nx/2 mode
+    Σ += 2 * sum(uh[2:grid.nkr-1])  # sum twice for 0 < k < nx/2 modes
   else # count every mode once
-    U = sum(uh)
+    Σ = sum(uh)
   end
 
-  norm = grid.Lx / grid.nx^2 # normalization for dft
+  Σ *= grid.Lx / grid.nx^2 # normalization for dft
 
-  return norm * real(U)
+  return real(Σ)
 end
 
 """
@@ -219,14 +226,15 @@ compute the radial spectrum, we first interpolate ``f̂(k, l)`` onto a radial wa
 ``f̂ =`` `fh` `` / (n_x e^{i 𝐤 ⋅ 𝐱₀})``. The elements of the vector ``𝐱₀`` are the
 left-most position in each direction, e.g., for a 2D grid `(grid.x[1], grid.y[1])`.
 
-After interpolation, we integrate ``f̂``over angles ``θ`` to get `fρ`,
+After interpolation, we integrate ``f̂`` over angles ``θ`` to get `fρ`,
+
 ```math
-  f̂_ρ = \\int f̂(ρ, θ) ρ 𝖽ρ 𝖽θ \\, .
+  f̂_ρ = \\int f̂(ρ, θ) ρ 𝖽ρ 𝖽θ .
 ```
 
 The resolution `(n, m)` for the polar wavenumber grid is `n = refinement * maximum(grid.nk, grid.nl), 
-m = refinement * maximum(grid.nk, grid.nl)`, where `refinement = 2` by default. If `fh` is in conjugate 
-symmetric form then only the upper-half plane in ``θ`` is represented on the polar grid.
+m = refinement * maximum(grid.nk, grid.nl)`, where `refinement = 2` by default. If `fh` is in
+conjugate symmetric form then only the upper-half plane in ``θ`` is represented on the polar grid.
 """
 function radialspectrum(fh, grid::TwoDGrid; n=nothing, m=nothing, refinement=2)
 
@@ -280,15 +288,15 @@ end
 Return an array, of the type compatible with the `device` that the `grid` lives on,
 that contains the values of function `func` evaluated on the `grid`.
 """
-function on_grid(func, grid::OneDGrid{T}) where T
+function on_grid(func::Function, grid::OneDGrid{T}) where T
   f = zeros(grid.device, T, (grid.nx, ))
 
   @. f = func(grid.x)
 
   return f
 end
 
-function on_grid(func, grid::TwoDGrid{T}) where T
+function on_grid(func::Function, grid::TwoDGrid{T}) where T
   f = zeros(grid.device, T, (grid.nx, grid.ny))
 
   x = reshape(grid.x, (grid.nx, 1))