Add background zonal flow U or U(y) in the SingleLayerQG module

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 contributors"]
-version = "0.16.1"
+version = "0.16.2"
 
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"

docs/src/modules/multilayerqg.md

@@ -71,10 +71,14 @@ with
 ```math
 \begin{aligned}
 \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 ,
 \end{aligned}
 ```
 
+the background PV gradient components in each layer and with
+``\mathsf{J}(a, b) = (\partial_x a)(\partial_y b)-(\partial_y a)(\partial_x b)`` is the 
+two-dimensional Jacobian. On the right hand side, ``\mu`` is linear bottom drag, and ``\nu`` is
+hyperviscosity of order ``n_\nu``. Plain old viscosity corresponds to ``n_\nu = 1``.
 
 ### Implementation
 

docs/src/modules/singlelayerqg.md

@@ -20,33 +20,50 @@ radius follows is said to obey equivalent-barotropic dynamics. We denote the sum
 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:
+The dynamical variable is ``q``. Including an imposed zonal flow ``U(y)``, the equation of motion 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) + (U - \partial_y\psi) \partial_x Q +  U \partial_x q + (\partial_y Q)(\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)`` is the 
-two-dimensional Jacobian. On the right hand side, ``F(x, y, t)`` is forcing, ``\mu`` is 
+with
+
+```math
+\begin{aligned}
+\partial_y Q &\equiv \beta - \partial_y^2 U + \partial_y \eta , \\
+\partial_x Q &\equiv \partial_x \eta ,
+\end{aligned}
+```
+
+the background PV gradient components, and with
+``\mathsf{J}(a, b) = (\partial_x a)(\partial_y b) - (\partial_y a)(\partial_x b)``
+the two-dimensional Jacobian. 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``.
 
+In the case that the imposed background zonal flow is just a constant, the above simplifies to:
+
+```math
+\partial_t q + \mathsf{J}(\psi, q + \eta) + U \partial_x (q + \eta) + β \partial_x \psi = \underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu} \right] q}_{\textrm{dissipation}} + F ,
+```
+
+and thus the advection of ``q`` can be incorporated in the linear term ``L``.
 
 ### Implementation
 
 The equation is time-stepped forward in Fourier space:
 
 ```math
-\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} .
+\partial_t \widehat{q} = - \widehat{\mathsf{J}(\psi, q)} - \widehat{U \partial_x Q} - \widehat{U \partial_x q}
++ \widehat{(\partial_y \psi) (\partial_x Q)}  - \widehat{(\partial_x \psi)(\partial_y 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) =
+\partial_y [ (\partial_x f) g] - \partial_x[ (\partial_y f) g]``.
+
 The state variable `sol` is the Fourier transform of the sum of relative vorticity and vortex 
 stretching (when the latter is applicable), [`qh`](@ref GeophysicalFlows.SingleLayerQG.Vars).
 
-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 linear operator is constructed in `Equation`
 
 ```@docs
@@ -122,4 +139,4 @@ Other diagnostic include: [`energy_dissipation`](@ref GeophysicalFlows.SingleLay
   (barotropic quasi-geostrophic flow with ``\beta=0``) above topography.
 
 - [`examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl`](@ref singlelayerqg_decaying_barotropic_equivalentbarotropic_example):
-  Simulate two dimensional turbulence (``\beta=0``) with both infinite and finite Rossby radius of deformation and compares the evolution of the two.
+  Simulate two dimensional turbulence (``\beta=0``) with both infinite and finite Rossby radius of deformation and compares the evolution of the two.

src/multilayerqg.jl

@@ -66,7 +66,7 @@ Keyword arguments
   - `Ly`: Extent of the ``y``-domain.
   - `f₀`: Constant planetary vorticity.
   - `β`: Planetary vorticity ``y``-gradient.
-  - `U`: The imposed constant zonal flow ``U(y)`` in each fluid layer.
+  - `U`: Imposed background constant zonal flow ``U(y)`` in each fluid layer.
   - `H`: Rest height of each fluid layer.
   - `b`: Boussinesq buoyancy of each fluid layer.
   - `eta`: Periodic component of the topographic potential vorticity.
@@ -757,7 +757,7 @@ function calcN_advection!(N, sol, vars, params, grid)
   @. vars.u += params.U                    # add the imposed zonal flow U
 
   uQx, uQxh = vars.q, vars.uh              # use vars.q and vars.uh as scratch variables
-  @. uQx  = vars.u * params.Qx             # (U+u)*∂Q/∂x
+  @. uQx = vars.u * params.Qx              # (U+u)*∂Q/∂x
   fwdtransform!(uQxh, uQx, params)
   @. N = - uQxh                            # -\hat{(U+u)*∂Q/∂x}
 
@@ -803,6 +803,7 @@ function calcN_linearadvection!(N, sol, vars, params, grid)
   @. vars.vh =  im * grid.kr * vars.ψh
 
   invtransform!(vars.u, vars.uh, params)
+
   @. vars.u += params.U                    # add the imposed zonal flow U
   uQx, uQxh = vars.q, vars.uh              # use vars.q and vars.uh as scratch variables
   @. uQx  = vars.u * params.Qx             # (U+u)*∂Q/∂x