Merge pull request #343 from mpudig/density_to_buoyancy

Accepting buoyancy instead of density in `MultiLayerQG` module
FourierFlows · Dec 20, 2023 · 449f8be · 449f8be · navidcy · Dec 20, 2023
2 parents 6e9095c + cd22ddc
commit 449f8be
Show file tree

Hide file tree

Showing 6 changed files with 75 additions and 75 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,8 +1,8 @@
 name = "GeophysicalFlows"
 uuid = "44ee3b1c-bc02-53fa-8355-8e347616e15e"
 license = "MIT"
-authors = ["Navid C. Constantinou <navidcy@gmail.com>", "Gregory L. Wagner <wagner.greg@gmail.com>", "and co-contributors"]
-version = "0.15.4"
+authors = ["Navid C. Constantinou <navidcy@gmail.com>", "Gregory L. Wagner <wagner.greg@gmail.com>", "and contributors"]
+version = "0.16.0"
 
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -84,7 +84,7 @@ The bibtex entry for the paper is:
 
 ## Papers using `GeophysicalFlows.jl`
 
-1. Drivas, T. D. and Elgindi, T. M. (2023). Singularity formation in the incompressible Euler equation in finite and infinite time. _EMS Surv. Math. Sci._, **10(1)**, 1–100, doi:[10.4171/emss/66](https://doi.org/10.4171/emss/66).
+1. Drivas, T. D. and Elgindi, T. M. (2023). Singularity formation in the incompressible Euler equation in finite and infinite time. _EMS Surveys in Mathematical Sciences_, **10(1)**, 1–100, doi:[10.4171/emss/66](https://doi.org/10.4171/emss/66).
 
 1. Parfenyev, V. (2023) Statistical analysis of vortex condensate motion in two-dimensional turbulence. arXiv preprint arXiv:2311.03065, doi:[10.48550/arXiv.2311.03065](https://doi.org/10.48550/arXiv.2311.03065).
 

diff --git a/docs/src/modules/multilayerqg.md b/docs/src/modules/multilayerqg.md
@@ -2,7 +2,7 @@
 
 ### Basic Equations
 
-This module solves the layered quasi-geostrophic equations on a beta plane of variable fluid 
+This module solves the layered Boussinesq quasi-geostrophic equations on a beta plane of variable fluid 
 depth ``H - h(x, y)``. The flow in each layer is obtained through a streamfunction ``\psi_j`` as 
 ``(u_j, v_j) = (-\partial_y \psi_j, \partial_x \psi_j)``, ``j = 1, \dots, n``, where ``n`` 
 is the number of fluid layers.
@@ -28,9 +28,17 @@ with
 
 ```math
 F_{j+1/2, k} = \frac{f_0^2}{g'_{j+1/2} H_k} \quad \text{and} \quad
-g'_{j+1/2} = g \frac{\rho_{j+1} - \rho_j}{\rho_{j+1}} .
+g'_{j+1/2} = b_j - b_{j+1} ,
 ```
 
+where
+
+```math
+b_{j} = - g \frac{\delta \rho_j}{\rho_0}
+```
+
+is the Boussinesq buoyancy in each layer, with ``\rho = \rho_0 + \delta \rho`` the total density, ``\rho_0`` a constant reference density, and ``|\delta \rho| \ll \rho_0`` the perturbation density.
+
 In view of the relationships above, when we convert to Fourier space ``q``'s and ``\psi``'s are 
 related via the matrix equation
 

diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl
@@ -43,9 +43,9 @@ L = 2π                   # domain size
 β = 5                    # the y-gradient of planetary PV
 
 nlayers = 2              # number of layers
-f₀, g = 1, 1             # Coriolis parameter and gravitational constant
+f₀ = 1                   # Coriolis parameter
 H = [0.2, 0.8]           # the rest depths of each layer
-ρ = [4.0, 5.0]           # the density of each layer
+b = [-1.0, -1.2]         # Boussinesq buoyancy of each layer
 
 U = zeros(nlayers)       # the imposed mean zonal flow in each layer
 U[1] = 1.0
@@ -59,7 +59,7 @@ nothing #hide
 # at every time-step that removes enstrophy at high wavenumbers and, thereby,
 # stabilize the problem, despite that we use the default viscosity coefficient `ν=0`.
 
-prob = MultiLayerQG.Problem(nlayers, dev; nx=n, Lx=L, f₀, g, H, ρ, U, μ, β,
+prob = MultiLayerQG.Problem(nlayers, dev; nx=n, Lx=L, f₀, H, b, U, μ, β,
                             dt, stepper, aliased_fraction=0)
 nothing #hide
 

diff --git a/src/multilayerqg.jl b/src/multilayerqg.jl
@@ -35,10 +35,9 @@ nothingfunction(args...) = nothing
                          Ly = Lx,
                          f₀ = 1.0,
                           β = 0.0,
-                          g = 1.0,
                           U = zeros(nlayers),
                           H = 1/nlayers * ones(nlayers),
-                          ρ = Array{Float64}(1:nlayers),
+                          b = -(1 .+ 1/nlayers * Array{Float64}(0:nlayers-1)),
                         eta = nothing,
     topographic_pv_gradient = (0, 0),
                           μ = 0,
@@ -67,10 +66,9 @@ Keyword arguments
   - `Ly`: Extent of the ``y``-domain.
   - `f₀`: Constant planetary vorticity.
   - `β`: Planetary vorticity ``y``-gradient.
-  - `g`: Gravitational acceleration constant.
   - `U`: The imposed constant zonal flow ``U(y)`` in each fluid layer.
   - `H`: Rest height of each fluid layer.
-  - `ρ`: Density of each fluid layer.
+  - `b`: Boussinesq buoyancy of each fluid layer.
   - `eta`: Periodic component of the topographic potential vorticity.
   - `topographic_pv_gradient`: The ``(x, y)`` components of the topographic PV large-scale gradient.
   - `μ`: Linear bottom drag coefficient.
@@ -92,14 +90,13 @@ function Problem(nlayers::Int,                             # number of fluid lay
                            Lx = 2π,
                            Ly = Lx,
               # Physical parameters
-                           f₀ = 1.0,                       # Coriolis parameter
-                            β = 0.0,                       # y-gradient of Coriolis parameter
-                            g = 1.0,                       # gravitational constant
-                            U = zeros(nlayers),            # imposed zonal flow U(y) in each layer
-                            H = 1/nlayers * ones(nlayers), # rest fluid height of each layer
-                            ρ = Array{Float64}(1:nlayers), # density of each layer
-                          eta = nothing,                   # periodic component of the topographic PV
-      topographic_pv_gradient = (0, 0),                    # tuple with the ``(x, y)`` components of topographic PV large-scale gradient
+                           f₀ = 1.0,                                             # Coriolis parameter
+                            β = 0.0,                                             # y-gradient of Coriolis parameter
+                            U = zeros(nlayers),                                  # imposed zonal flow U(y) in each layer
+                            H = 1/nlayers * ones(nlayers),                       # rest fluid height of each layer
+                            b = -(1 .+ 1/nlayers * Array{Float64}(0:nlayers-1)), # Boussinesq buoyancy of each layer
+                          eta = nothing,                                         # periodic component of the topographic PV
+      topographic_pv_gradient = (0, 0),                                          # tuple with the ``(x, y)`` components of topographic PV large-scale gradient
               # Bottom Drag and/or (hyper)-viscosity
                             μ = 0,
                             ν = 0,
@@ -137,7 +134,7 @@ function Problem(nlayers::Int,                             # number of fluid lay
 
   grid = TwoDGrid(dev; nx, Lx, ny, Ly, aliased_fraction, T)
 
-  params = Params(nlayers, g, f₀, β, ρ, H, U, eta, topographic_pv_gradient, μ, ν, nν, grid; calcFq)
+  params = Params(nlayers, f₀, β, b, H, U, eta, topographic_pv_gradient, μ, ν, nν, grid; calcFq)
 
   vars = calcFq == nothingfunction ? DecayingVars(grid, params) : (stochastic ? StochasticForcedVars(grid, params) : ForcedVars(grid, params))
 
@@ -157,14 +154,12 @@ struct Params{T, Aphys3D, Aphys2D, Atrans4D, Trfft} <: AbstractParams
   # prescribed params
     "number of fluid layers"
    nlayers :: Int
-    "gravitational constant"
-         g :: T
     "constant planetary vorticity"
         f₀ :: T
     "planetary vorticity ``y``-gradient"
          β :: T
-    "array with density of each fluid layer"
-         ρ :: Tuple
+    "array with Boussinesq buoyancy of each fluid layer"
+         b :: Tuple
     "array with rest height of each fluid layer"
          H :: Tuple
     "array with imposed constant zonal flow ``U(y)`` in each fluid layer"
@@ -241,14 +236,12 @@ $(TYPEDFIELDS)
 """
 struct TwoLayerParams{T, Aphys3D, Aphys2D, Trfft} <: AbstractParams
   # prescribed params
-    "gravitational constant"
-         g :: T
     "constant planetary vorticity"
         f₀ :: T
     "planetary vorticity ``y``-gradient"
          β :: T
-    "array with density of each fluid layer"
-         ρ :: Tuple
+    "array with Boussinesq buoyancy of each fluid layer"
+         b :: Tuple
     "tuple with rest height of each fluid layer"
          H :: Tuple
    "array with imposed constant zonal flow ``U(y)`` in each fluid layer"
@@ -311,7 +304,7 @@ function convert_U_to_U3D(dev, nlayers, grid, U::Number)
   return A(U_3D)
 end
 
-function Params(nlayers::Int, g, f₀, β, ρ, H, U, eta, topographic_pv_gradient, μ, ν, nν, grid::TwoDGrid; calcFq=nothingfunction, effort=FFTW.MEASURE)
+function Params(nlayers::Int, f₀, β, b, H, U, eta, topographic_pv_gradient, μ, ν, nν, grid::TwoDGrid; calcFq=nothingfunction, effort=FFTW.MEASURE)
   dev = grid.device
   T = eltype(grid)
   A = device_array(dev)
@@ -349,10 +342,10 @@ function Params(nlayers::Int, g, f₀, β, ρ, H, U, eta, topographic_pv_gradien
 
   else # if nlayers≥2
 
-    ρ = reshape(T.(ρ), (1,  1, nlayers))
+    b = reshape(T.(b), (1,  1, nlayers))
     H = Tuple(T.(H))
 
-    g′ = T(g) * (ρ[2:nlayers] -   ρ[1:nlayers-1]) ./ ρ[2:nlayers] # reduced gravity at each interface
+    g′ = b[1:nlayers-1] - b[2:nlayers] # reduced gravity at each interface
 
     Fm = @. T(f₀^2 / (g′ * H[2:nlayers]))
     Fp = @. T(f₀^2 / (g′ * H[1:nlayers-1]))
@@ -374,9 +367,9 @@ function Params(nlayers::Int, g, f₀, β, ρ, H, U, eta, topographic_pv_gradien
     CUDA.@allowscalar @views Qy[:, :, nlayers] = @. Qy[:, :, nlayers] - Fm[nlayers-1] * (U[:, :, nlayers-1] - U[:, :, nlayers])
 
     if nlayers==2
-      return TwoLayerParams(T(g), T(f₀), T(β), Tuple(T.(ρ)), Tuple(T.(H)), U, eta, topographic_pv_gradient, T(μ), T(ν), nν, calcFq, T(g′[1]), Qx, Qy, rfftplanlayered)
+      return TwoLayerParams(T(f₀), T(β), Tuple(T.(b)), Tuple(T.(H)), U, eta, topographic_pv_gradient, T(μ), T(ν), nν, calcFq, T(g′[1]), Qx, Qy, rfftplanlayered)
     else # if nlayers>2
-      return Params(nlayers, T(g), T(f₀), T(β), Tuple(T.(ρ)), T.(H), U, eta, topographic_pv_gradient, T(μ), T(ν), nν, calcFq, Tuple(T.(g′)), Qx, Qy, S, S⁻¹, rfftplanlayered)
+      return Params(nlayers, T(f₀), T(β), Tuple(T.(b)), T.(H), U, eta, topographic_pv_gradient, T(μ), T(ν), nν, calcFq, Tuple(T.(g′)), Qx, Qy, S, S⁻¹, rfftplanlayered)
     end
   end
 end