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

…#360) * Update singlelayerqg.jl to include background flow * Fixed typo in previous change * added background flow test * corrected typo * Added some spaces to cheer up Greg Co-authored-by: Gregory L. Wagner <[email protected]> * more spaces * added variable U functionality * reduced memory of Uy assignment and updated function descriptions * update function descriptions * T(U) -> convert(T, U) Co-authored-by: Gregory L. Wagner <[email protected]> * clarify U array formulation Co-authored-by: Gregory L. Wagner <[email protected]> * fix bug with previous push * add additional test function for U as an Array and include tests in * simplify background U and use Rossby wave test for U that does not vary in y * homogenize notation/docs * bump patch release * fix test for gpu * == -> === * use different calcN_advection! methods for U const and y-varying * minor tweaks * update singlelayergq module in docs * homogenize notation/wording * add empty line * code alignment * add Qx and Qy in params; refactor calcN_advection * clarify when U=const * some fixes + cleanup * add nonlinear advection test for finite deformation + topo * fixes nonlinear advection w U(y) + add tests * add comment for β term * remove debug statements + fix test for julia 1.6 * more tests for the nonlinear terms * fix test_1layerqg_nonlinearadvection on gpu * fix test_1layerqg_nonlinearadvection on gpu * improve docs * drop GeophysicalFlows from docs/Project.toml --------- Co-authored-by: Gregory L. Wagner <[email protected]> Co-authored-by: Navid C. Constantinou <[email protected]>
FourierFlows · Jul 29, 2024 · 14f581b · 14f581b · navidcy · Jul 29, 2024
1 parent 467e82f
commit 14f581b
Show file tree

Hide file tree

Showing 7 changed files with 340 additions and 112 deletions.
diff --git a/Project.toml b/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"

diff --git a/docs/src/modules/multilayerqg.md b/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
 

diff --git a/docs/src/modules/singlelayerqg.md b/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.
diff --git a/src/multilayerqg.jl b/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