Docs: simplify notation + update Docs/Forcing section + remove Ito comments from modules/tests

docs/make.jl

@@ -94,10 +94,11 @@ sitename = "GeophysicalFlows.jl",
               "modules/multilayerqg.md",
               "modules/surfaceqg.md"
             ],
-            "Forcing" => "forcing.md",
+            "Stochastic Forcing" => "stochastic_forcing.md",
             "DocStrings" => Any[
             "man/types.md",
-            "man/functions.md"]
+            "man/functions.md"
+            ]
            ]
 )
 

docs/src/assets/README.md

@@ -0,0 +1 @@
+### a placeholder directory for output generated by Docs

docs/src/modules/barotropicqgql.md

@@ -2,54 +2,75 @@
 
 ### Basic Equations
 
-This module solves the *quasi-linear* quasi-geostrophic barotropic vorticity equation on a beta-plane of variable fluid depth $H-h(x,y)$. 
-Quasi-linear refers to the dynamics that *neglect* the eddy--eddy interactions in the eddy evolution equation after an eddy--mean flow decomposition, e.g., 
+This module solves the *quasi-linear* quasi-geostrophic barotropic vorticity equation on a beta 
+plane of variable fluid depth ``H - h(x, y)``. Quasi-linear refers to the dynamics that *neglect* 
+the eddy--eddy interactions in the eddy evolution equation after an eddy--mean flow decomposition, e.g., 
 
-$$\phi(x, y, t) = \overline{\phi}(y, t) + \phi'(x,y,t) ,$$
+```math
+\phi(x, y, t) = \overline{\phi}(y, t) + \phi'(x, y, t) ,
+```
 
-where overline above denotes a zonal mean, $\overline{\phi}(y, t) = \int \phi(x, y, t)\,\mathrm{d}x/L_x$, and prime denotes deviations from the zonal mean. This approximation is used in many process-model studies of zonation, e.g., 
+where overline above denotes a zonal mean, ``\overline{\phi}(y, t) = \int \phi(x, y, t) \, 𝖽x / L_x``, and prime denotes deviations from the zonal mean. This approximation is used in many process-model studies of zonation, e.g., 
 
 - Farrell, B. F. and Ioannou, P. J. (2003). [Structural stability of turbulent jets.](http://doi.org/10.1175/1520-0469(2003)060<2101:SSOTJ>2.0.CO;2) *J. Atmos. Sci.*, **60**, 2101-2118.
 - Srinivasan, K. and Young, W. R. (2012). [Zonostrophic instability.](http://doi.org/10.1175/JAS-D-11-0200.1) *J. Atmos. Sci.*, **69 (5)**, 1633-1656.
 - Constantinou, N. C., Farrell, B. F., and Ioannou, P. J. (2014). [Emergence and equilibration of jets in beta-plane turbulence: applications of Stochastic Structural Stability Theory.](http://doi.org/10.1175/JAS-D-13-076.1) *J. Atmos. Sci.*, **71 (5)**, 1818-1842.
 
-
-As in the [SingleLayerQG module](singlelayerqg.md), the flow is obtained through a streamfunction $\psi$ as $(u, v) = (-\partial_y\psi, \partial_x\psi)$. All flow fields can be obtained from the quasi-geostrophic potential vorticity (QGPV). Here the QGPV is
-
-$$\underbrace{f_0 + \beta y}_{\text{planetary PV}} + \underbrace{(\partial_x v
-	- \partial_y u)}_{\text{relative vorticity}} +
-	\underbrace{\frac{f_0 h}{H}}_{\text{topographic PV}}.$$
-
-The dynamical variable is the component of the vorticity of the flow normal to the plane of motion, $\zeta\equiv \partial_x v- \partial_y u = \nabla^2\psi$. Also, we denote the topographic PV with $\eta\equiv f_0 h/H$. After we apply the eddy-mean flow decomposition above, the QGPV dynamics are:
-
-$$\partial_t \overline{\zeta} + \mathsf{J}(\overline{\psi}, \underbrace{\overline{\zeta} + \overline{\eta}}_{\equiv \overline{q}}) + \overline{\mathsf{J}(\psi', \underbrace{\zeta' + \eta'}_{\equiv q'})} = \underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu}
-\right] \overline{\zeta} }_{\textrm{dissipation}} \ .$$
-
-$$\partial_t \zeta' + \mathsf{J}(\psi', \overline{q}) + \mathsf{J}(\overline{\psi}, q') + \underbrace{\mathsf{J}(\psi', q') - \overline{\mathsf{J}(\psi', q')}}_{\textrm{EENL}} + 
-\beta\partial_x\psi' = \underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu}
-\right] \zeta'}_{\textrm{dissipation}} + f\ .$$
-
-where $\mathsf{J}(a, b) = (\partial_x a)(\partial_y b)-(\partial_y a)(\partial_x b)$. On the right hand side, $f(x,y,t)$ is forcing (which is assumed to have zero mean, $\overline{f}=0$), $\mu$ is linear drag, and $\nu$ is hyperviscosity. Plain old viscosity corresponds to $n_{\nu}=1$. The sum of relative vorticity and topographic PV is denoted with $q\equiv\zeta+\eta$.
+As in the [SingleLayerQG module](singlelayerqg.md), the flow is obtained through a 
+streamfunction ``\psi`` as ``(u, v) = (-\partial_y \psi, \partial_x \psi)``. All flow fields 
+can be obtained from the quasi-geostrophic potential vorticity (QGPV). Here, the QGPV is
+
+```math
+\underbrace{f_0 + \beta y}_{\text{planetary PV}} + \underbrace{\partial_x v
+	- \partial_y u}_{\text{relative vorticity}} + \underbrace{\frac{f_0 h}{H}}_{\text{topographic PV}} .
+```
+
+The dynamical variable is the component of the vorticity of the flow normal to the plane of 
+motion, ``\zeta \equiv \partial_x v - \partial_y u = \nabla^2 \psi``. Also, we denote the 
+topographic PV with ``\eta \equiv f_0 h/H``. After we apply the eddy-mean flow decomposition 
+above, the QGPV dynamics are:
+
+```math
+\begin{aligned}
+\partial_t \overline{\zeta} & + \mathsf{J}(\overline{\psi}, \underbrace{\overline{\zeta} + \overline{\eta}}_{\equiv \overline{q}}) + \overline{\mathsf{J}(\psi', \underbrace{\zeta' + \eta'}_{\equiv q'})} = \underbrace{- \left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu}
+\right] \overline{\zeta} }_{\textrm{dissipation}} ,\\
+\partial_t \zeta' &+ \mathsf{J}(\psi', \overline{q}) + \mathsf{J}(\overline{\psi}, q') + \underbrace{\mathsf{J}(\psi', q') - \overline{\mathsf{J}(\psi', q')}}_{\textrm{EENL}} + 
+\beta \partial_x \psi' = \underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu}
+\right] \zeta'}_{\textrm{dissipation}} + F .
+\end{aligned}
+```
+
+where ``\mathsf{J}(a, b) = (\partial_x a)(\partial_y b) - (\partial_y a)(\partial_x b)``. On 
+the right hand side, ``F(x, y, t)`` is forcing (which is assumed to have zero zonal mean, 
+``\overline{F} = 0``), ``\mu`` is linear drag, and ``\nu`` is hyperviscosity. Plain old 
+viscosity corresponds to ``n_{\nu} = 1``. The sum of relative vorticity and topographic PV is 
+denoted with ``q \equiv \zeta + \eta``.
 
 *Quasi-linear* dynamics **neglect the term eddy-eddy nonlinearity (EENL) term** above.
 
 ### Implementation
 
 The equation is time-stepped forward in Fourier space:
 
-$$\partial_t \widehat{\zeta} = - \widehat{\mathsf{J}(\psi, q)}^{\textrm{QL}} +\beta\frac{\mathrm{i}k_x}{k^2}\widehat{\zeta} -\left(\mu
-+\nu k^{2n_\nu}\right) \widehat{\zeta}  + \widehat{f}\ .$$
+```math
+\partial_t \widehat{\zeta} = - \widehat{\mathsf{J}(\psi, q)}^{\textrm{QL}} + \beta \frac{i k_x}{|𝐤|^2} \widehat{\zeta} - \left ( \mu + \nu |𝐤|^{2n_\nu} \right ) \widehat{\zeta} + \widehat{F} .
+```
 
-In doing so the Jacobian is computed in the conservative form: $\mathsf{J}(f,g) =
-\partial_y [ (\partial_x f) g] -\partial_x[ (\partial_y f) g]$. The superscript QL in the Jacobian term above denotes that remove triad interactions that correspond to the EENL term.
+In doing so the Jacobian is computed in the conservative form: ``\mathsf{J}(f,g) =
+\partial_y [ (\partial_x f) g] -\partial_x[ (\partial_y f) g]``. The superscript QL on the 
+Jacobian term above denotes that triad interactions that correspond to the EENL term are 
+removed.
 
 Thus:
 
-$$\mathcal{L} = \beta\frac{\mathrm{i}k_x}{k^2} - \mu - \nu k^{2n_\nu}\ ,$$
-$$\mathcal{N}(\widehat{\zeta}) = - \mathrm{i}k_x \mathrm{FFT}(u q)^{\textrm{QL}}-
-	\mathrm{i}k_y \mathrm{FFT}(v q)^{\textrm{QL}}\ .$$
+```math
+\begin{aligned}
+L & = \beta \frac{i k_x}{|𝐤|^2} - \mu - \nu |𝐤|^{2n_\nu} ,\\
+N(\widehat{\zeta}) & = - i k_x \mathrm{FFT}(u q)^{\textrm{QL}} - i k_y \mathrm{FFT}(v q)^{\textrm{QL}} + \widehat{F}.
+\end{aligned}
+```
 
 
 ## Examples
 
-- `examples/barotropicqgql_betaforced.jl`: A script that simulates forced-dissipative quasi-linear quasi-geostrophic flow on a beta-plane demonstrating zonation. The forcing is temporally delta-correlated and its spatial structure is isotropic with power in a narrow annulus of total radius `kf` in wavenumber space. This example demonstrates that the anisotropic inverse energy cascade is not required for zonation.
+- `examples/barotropicqgql_betaforced.jl`: A script that simulates forced-dissipative quasi-linear quasi-geostrophic flow on a beta-plane demonstrating zonation. The forcing is temporally delta-correlated and its spatial structure is isotropic with power in a narrow annulus of total radius ``k_f`` in wavenumber space. This example demonstrates that the anisotropic inverse energy cascade is not required for zonation.

docs/src/modules/multilayerqg.md

@@ -10,10 +10,10 @@ is the number of fluid layers.
 The QGPV in each layer is
 
 ```math
-\mathrm{QGPV}_j = q_j  + \underbrace{f_0 + \beta y}_{\textrm{planetary PV}} + \delta_{j,n} \underbrace{\frac{f_0 h}{H_n}}_{\textrm{topographic PV}}, \quad j = 1, \dots, n \ .
+\mathrm{QGPV}_j = q_j + \underbrace{f_0 + \beta y}_{\textrm{planetary PV}} + \delta_{j, n} \underbrace{\frac{f_0 h}{H_n}}_{\textrm{topographic PV}}, \quad j = 1, \dots, n \ .
 ```
 
-where ``q_j`` incorporates the relative vorticity in each layer ``\nabla^2\psi_j`` and the 
+where ``q_j`` incorporates the relative vorticity in each layer ``\nabla^2 \psi_j`` and the 
 vortex stretching terms:
 
 ```math
@@ -33,9 +33,9 @@ In view of the relationships above, when we convert to Fourier space ``q``'s and
 related via the matrix equation
 
 ```math
-\begin{pmatrix} \widehat{q}_{\boldsymbol{k}, 1}\\\vdots\\\widehat{q}_{\boldsymbol{k}, n} \end{pmatrix} =
-\underbrace{\left(-|\boldsymbol{k}|^2 \mathbb{1} + \mathbb{F} \right)}_{\equiv \mathbb{S}_{\boldsymbol{k}}}
-\begin{pmatrix} \widehat{\psi}_{\boldsymbol{k}, 1}\\\vdots\\\widehat{\psi}_{\boldsymbol{k}, n} \end{pmatrix}
+\begin{pmatrix} \widehat{q}_{𝐤, 1}\\\vdots\\\widehat{q}_{𝐤, n} \end{pmatrix} =
+\underbrace{\left(-|𝐤|^2 \mathbb{1} + \mathbb{F} \right)}_{\equiv \mathbb{S}_{𝐤}}
+\begin{pmatrix} \widehat{\psi}_{𝐤, 1}\\\vdots\\\widehat{\psi}_{𝐤, n} \end{pmatrix}
 ```
 
 where
@@ -53,16 +53,26 @@ where
 Including an imposed zonal flow ``U_j(y)`` in each layer, the equations of motion are:
 
 ```math
-\partial_t q_j + \mathsf{J}(\psi_j, q_j ) + (U_j - \partial_y\psi_j) \partial_x Q_j +  U_j \partial_x q_j  + (\partial_y Q_j)(\partial_x \psi_j) = -\delta_{j, n} \mu \nabla^2 \psi_n - \nu (-1)^{n_\nu} \nabla^{2n_\nu} q_j,
+\partial_t q_j + \mathsf{J}(\psi_j, q_j ) + (U_j - \partial_y\psi_j) \partial_x Q_j +  U_j \partial_x q_j  + (\partial_y Q_j)(\partial_x \psi_j) = -\delta_{j, n} \mu \nabla^2 \psi_n - \nu (-1)^{n_\nu} \nabla^{2 n_\nu} q_j ,
 ```
 
 with
 
 ```math
 \partial_y Q_j \equiv \beta - \partial_y^2 U_j - (1-\delta_{j,1}) F_{j-1/2, j} (U_{j-1} - U_j) - (1 - \delta_{j,n}) F_{j+1/2, j} (U_{j+1} - U_j) + \delta_{j,n} \partial_y \eta \ , \\
-\partial_x Q_j \equiv \delta_{j,n} \partial_x \eta \ .
+\partial_x Q_j \equiv \delta_{j, n} \partial_x \eta \ .
 ```
 
+### Helper functions
+
+```@docs
+GeophysicalFlows.MultiLayerQG.set_q!
+GeophysicalFlows.MultiLayerQG.set_ψ!
+GeophysicalFlows.MultiLayerQG.updatevars!
+```
+
+### Diagnostics
+
 The eddy kinetic energy in each layer and the eddy potential energy that corresponds to each 
 fluid interface is computed via `energies()`:
 
@@ -79,37 +89,50 @@ GeophysicalFlows.MultiLayerQG.fluxes
 
 ### Implementation
 
-Matrices ``\mathbb{S}_{\boldsymbol{k}}`` as well as ``\mathbb{S}^{-1}_{\boldsymbol{k}}`` are included 
-in `params` as `params.S` and `params.S⁻¹` respectively. Additionally, the background PV gradients 
-``\partial_x Q`` and ``\partial_y Q`` are also included in the `params` as `params.Qx` and `params.Qy`.
-
-You can get ``\widehat{\psi}_j`` from ``\widehat{q}_j`` with `streamfunctionfrompv!(psih, qh, params, grid)`, 
-while to get ``\widehat{q}_j`` from ``\widehat{\psi}_j`` you need to call `pvfromstreamfunction!(qh, psih, params, grid)`.
+Matrices ``\mathbb{S}_{𝐤}`` as well as ``\mathbb{S}^{-1}_{𝐤}`` are included in `params` as 
+`params.S` and `params.S⁻¹` respectively. Additionally, the background PV gradients 
+``\partial_x Q`` and ``\partial_y Q`` are also included in the `params` as `params.Qx` and 
+`params.Qy`.
 
+One can get the ``\widehat{\psi}_j`` from ``\widehat{q}_j`` via 
+`streamfunctionfrompv!(psih, qh, params, grid)`, while the inverse, i.e. obtain ``\widehat{q}_j`` from ``\widehat{\psi}_j``, is done via  `pvfromstreamfunction!(qh, psih, params, grid)`.
 
 The equations of motion are time-stepped forward in Fourier space:
 
 ```math
 \partial_t \widehat{q}_j = - \widehat{\mathsf{J}(\psi_j, q_j)}  - \widehat{U_j \partial_x Q_j} - \widehat{U_j \partial_x q_j}
-+ \widehat{(\partial_y \psi_j) \partial_x Q_j}  - \widehat{(\partial_x\psi_j)(\partial_y Q_j)} + \delta_{j,n} \mu k^{2} \widehat{\psi}_n - \nu k^{2n_\nu} \widehat{q}_j \ .
++ \widehat{(\partial_y \psi_j) \partial_x Q_j}  - \widehat{(\partial_x \psi_j)(\partial_y Q_j)} + \delta_{j, n} \mu |𝐤|^{2} \widehat{\psi}_n - \nu |𝐤|^{2n_\nu} \widehat{q}_j \ .
 ```
 
 In doing so the Jacobian is computed in the conservative form: ``\mathsf{J}(f,g) =
-\partial_y [ (\partial_x f) g] -\partial_x[ (\partial_y f) g]``.
+\partial_y [ (\partial_x f) g] - \partial_x[ (\partial_y f) g]``.
 
 Equations are formulated using $q_j$ as the state variables, i.e., `sol = qh`.
 
-Thus:
+The linear operator is constructed in `Equation`
 
-```math
-\begin{aligned}
-\mathcal{L} & = - \nu k^{2n_\nu} \ , \\
-\mathcal{N}(\widehat{q}_j) & = - \widehat{\mathsf{J}(\psi_j, q_j)} - \widehat{U_j \partial_x Q_j} - \widehat{U_j \partial_x q_j}
- + \widehat{(\partial_y \psi_j)(\partial_x Q_j)} - \widehat{(\partial_x \psi_j)(\partial_y Q_j)} + \delta_{j,n} \mu k^{2} \widehat{\psi}_n\ .
-\end{aligned}
+```@docs
+GeophysicalFlows.MultiLayerQG.Equation
+GeophysicalFlows.MultiLayerQG.hyperviscosity
 ```
-
 
+The nonlinear terms is computed via
+
+```@docs
+GeophysicalFlows.MultiLayerQG.calcN!
+```
+which, in turn, calls 
+
+```@docs
+GeophysicalFlows.MultiLayerQG.calcN_advection!
+```
+and
+
+```@docs
+GeophysicalFlows.MultiLayerQG.addforcing!
+```
+
+
  ## Examples
 
  - `examples/multilayerqg_2layer.jl`: A script that simulates the growth and equilibration of baroclinic eddy turbulence in the Phillips 2-layer model.

docs/src/modules/singlelayerqg.md

@@ -3,8 +3,9 @@
 ### Basic Equations
 
 This module solves the barotropic or equivalent barotropic quasi-geostrophic vorticity equation 
-on a beta-plane of variable fluid depth ``H - h(x, y)``. The flow is obtained through a streamfunction ``\psi`` as ``(u, v) = (-\partial_y \psi, \partial_x \psi)``. All flow 
-fields can be obtained from the quasi-geostrophic potential vorticity (QGPV). Here the QGPV is
+on a beta-plane of variable fluid depth ``H - h(x, y)``. The flow is obtained through a 
+streamfunction ``\psi`` as ``(u, v) = (-\partial_y \psi, \partial_x \psi)``. All flow fields 
+can be obtained from the quasi-geostrophic potential vorticity (QGPV). Here the QGPV is
 
 ```math
 	\underbrace{f_0 + \beta y}_{\text{planetary PV}} + \underbrace{\partial_x v
@@ -16,19 +17,19 @@ fields can be obtained from the quasi-geostrophic potential vorticity (QGPV). He
 where ``\ell`` is the Rossby radius of deformation. Purely barotropic dynamics corresponds to 
 infinite Rossby radius of deformation (``\ell = \infty``), while a flow with a finite Rossby 
 radius follows is said to obey equivalent-barotropic dynamics. We denote the sum of the relative
-vorticity and the vortex stretching contributions to the QGPV with ``q = \nabla^2 \psi - \psi / \ell^2``.
+vorticity and the vortex stretching contributions to the QGPV with ``q \equiv \nabla^2 \psi - \psi / \ell^2``.
 Also, we denote the topographic PV with ``\eta \equiv f_0 h / H``.
 
 The dynamical variable is ``q``.  Thus, the equation solved by the module is:
 
 ```math
-\partial_t q + \mathsf{J}(\psi, q + \eta) +
-\beta \partial_x \psi = \underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu}
-\right] q}_{\textrm{dissipation}} + f \ .
+\partial_t q + \mathsf{J}(\psi, q + \eta) + \beta \partial_x \psi = 
+\underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu} \right] q}_{\textrm{dissipation}} + F \ .
 ```
 
 where ``\mathsf{J}(a, b) = (\partial_x a)(\partial_y b)-(\partial_y a)(\partial_x b)``. On 
-the right hand side, ``f(x, y, t)`` is forcing, ``\mu`` is linear drag, and ``\nu`` is hyperviscosity of order ``n_\nu``. Plain old viscosity corresponds to ``n_\nu = 1``.
+the right hand side, ``F(x, y, t)`` is forcing, ``\mu`` is linear drag, and ``\nu`` is 
+hyperviscosity of order ``n_\nu``. Plain old viscosity corresponds to ``n_\nu = 1``.
 
 The kinetic energy of the fluid is computed via:
 
@@ -53,7 +54,7 @@ GeophysicalFlows.SingleLayerQG.energy
 The equation is time-stepped forward in Fourier space:
 
 ```math
-\partial_t \widehat{q} = - \widehat{\mathsf{J}(\psi, q + \eta)} + \beta \frac{\mathrm{i} k_x}{k^2 + 1/\ell^2} \widehat{q} - \left(\mu + \nu k^{2n_\nu} \right) \widehat{q} + \widehat{f} \ .
+\partial_t \widehat{q} = - \widehat{\mathsf{J}(\psi, q + \eta)} + \beta \frac{i k_x}{|𝐤|^2 + 1/\ell^2} \widehat{q} - \left(\mu + \nu |𝐤|^{2n_\nu} \right) \widehat{q} + \widehat{F} \ .
 ```
 
 In doing so the Jacobian is computed in the conservative form: ``\mathsf{J}(f,g) =
@@ -63,8 +64,8 @@ Thus:
 
 ```math
 \begin{aligned}
-\mathcal{L} & = \beta \frac{\mathrm{i} k_x}{k^2 + 1/\ell^2} - \mu - \nu k^{2n_\nu} \ , \\
-\mathcal{N}(\widehat{q}) & = - \mathrm{i} k_x \mathrm{FFT}[u (q+\eta)] - \mathrm{i} k_y \mathrm{FFT}[v (q+\eta)] \ .
+L & = \beta \frac{i k_x}{|𝐤|^2 + 1/\ell^2} - \mu - \nu |𝐤|^{2n_\nu} \ , \\
+N(\widehat{q}) & = - i k_x \mathrm{FFT}[u (q + \eta)] - i k_y \mathrm{FFT}[v (q + \eta)] + \widehat{F}  \ .
 \end{aligned}
 ```