Skip to content

Commit

Permalink
Changes to multilayerqg to accept Boussinesq buoyancy instead of dens…
Browse files Browse the repository at this point in the history 
…ity.
  • Loading branch information
Matt Pudig committed Nov 20, 2023
1 parent e7414a4 commit 6d01319
Show file tree
Hide file tree
Showing 3 changed files with 50 additions and 42 deletions.
12 changes: 10 additions & 2 deletions docs/src/modules/multilayerqg.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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+1} - b_j) .
```

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

Expand Down
2 changes: 1 addition & 1 deletion examples/multilayerqg_2layer.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ L = 2π # domain size
nlayers = 2 # number of layers
f₀ = 1 # Coriolis parameter
H = [0.2, 0.8] # the rest depths of each layer
b = [1.0, 1.5] # Boussinesq buoyancy 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
Expand Down
78 changes: 39 additions & 39 deletions test/test_multilayerqg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,11 +130,11 @@ function test_pvtofromstreamfunction_2layer(dev::Device=CPU())
gr = TwoDGrid(dev; nx=n, Lx=L)

nlayers = 2 # these choice of parameters give the
f₀, g = 1, 1 # desired PV-streamfunction relations
f₀ = 1 # desired PV-streamfunction relations
H = [0.2, 0.8] # q1 = Δψ1 + 25*(ψ2-ψ1), and
ρ = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2).
b = [-1.0, -1.2] # q2 = Δψ2 + 25/4*(ψ1-ψ2).

prob = MultiLayerQG.Problem(nlayers, dev; nx=n, Lx=L, f₀, g, H, ρ)
prob = MultiLayerQG.Problem(nlayers, dev; nx=n, Lx=L, f₀, H, b)
sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid

ψ1, ψ2, q1, q2, ψ1x, ψ2x, q1x, q2x, Δψ2, Δq1, Δq2 = constructtestfields_2layer(gr)
Expand Down Expand Up @@ -172,13 +172,13 @@ function test_pvtofromstreamfunction_3layer(dev::Device=CPU())
n, L = 128, 2π
gr = TwoDGrid(dev; nx=n, Lx=L)

nlayers = 3 # these choice of parameters give the
f₀, g = 1, 1 # desired PV-streamfunction relations
H = [0.25, 0.25, 0.5] # q1 = Δψ1 + 20ψ2 - 20ψ1,
ρ = [4.0, 5.0, 6.0] # q2 = Δψ2 + 20ψ1 - 44ψ2 + 24ψ3,
# q3 = Δψ3 + 12ψ2 - 12ψ3.
nlayers = 3 # these choice of parameters give the
f₀ = 1 # desired PV-streamfunction relations
H = [0.25, 0.25, 0.5] # q1 = Δψ1 + 20ψ2 - 20ψ1,
b = [-1.0, -1.2, -1.2-1/6] # q2 = Δψ2 + 20ψ1 - 44ψ2 + 24ψ3,
# q3 = Δψ3 + 12ψ2 - 12ψ3.

prob = MultiLayerQG.Problem(nlayers, dev; nx=n, Lx=L, f₀, g, H, ρ)
prob = MultiLayerQG.Problem(nlayers, dev; nx=n, Lx=L, f₀, H, b)
sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid

ψ1, ψ2, ψ3, q1, q2, q3, ψ1x, ψ2x, ψ3x, q1x, q2x, q3x, Δq1, Δq2, Δq3, Δψ3 = constructtestfields_3layer(gr)
Expand Down Expand Up @@ -237,9 +237,9 @@ function test_mqg_nonlinearadvection_2layers(dt, stepper, dev::Device=CPU())
k₀, l₀ = 2π/gr.Lx, 2π/gr.Ly # fundamental wavenumbers

nlayers = 2 # these choice of parameters give the
f₀, g = 1, 1 # desired PV-streamfunction relations
f₀ = 1 # desired PV-streamfunction relations
H = [0.2, 0.8] # q1 = Δψ1 + 25*(ψ2-ψ1), and
ρ = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2).
b = [-1.0, -1.2] # q2 = Δψ2 + 25/4*(ψ1-ψ2).

β = 0.35

Expand Down Expand Up @@ -282,7 +282,7 @@ function test_mqg_nonlinearadvection_2layers(dt, stepper, dev::Device=CPU())
return nothing
end

prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, g, H, ρ, U,
prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, H, b, U,
eta=η, β, μ, ν, nν, calcFq=calcFq!, stepper, dt)

sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid
Expand Down Expand Up @@ -323,11 +323,11 @@ function test_mqg_nonlinearadvection_3layers(dt, stepper, dev::Device=CPU())
x, y = gridpoints(gr)
k₀, l₀ = 2π/gr.Lx, 2π/gr.Ly # fundamental wavenumbers

nlayers = 3 # these choice of parameters give the
f₀, g = 1, 1 # desired PV-streamfunction relations
H = [0.25, 0.25, 0.5] # q1 = Δψ1 + 20ψ2 - 20ψ1,
ρ = [4.0, 5.0, 6.0] # q2 = Δψ2 + 20ψ1 - 44ψ2 + 24ψ3,
# q3 = Δψ3 + 12ψ2 - 12ψ3.
nlayers = 3 # these choice of parameters give the
f₀ = 1 # desired PV-streamfunction relations
H = [0.25, 0.25, 0.5] # q1 = Δψ1 + 20ψ2 - 20ψ1,
b = [-1.0, -1.2, -1.2-1/6] # q2 = Δψ2 + 20ψ1 - 44ψ2 + 24ψ3,
# q3 = Δψ3 + 12ψ2 - 12ψ3.

β = 0.35

Expand Down Expand Up @@ -375,7 +375,7 @@ function test_mqg_nonlinearadvection_3layers(dt, stepper, dev::Device=CPU())
return nothing
end

prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, g, H, ρ, U,
prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, H, b, U,
eta=η, β, μ, ν, nν, calcFq=calcFq!, stepper, dt)

sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid
Expand Down Expand Up @@ -426,9 +426,9 @@ function test_mqg_linearadvection(dt, stepper, dev::Device=CPU();
k₀, l₀ = 2π/gr.Lx, 2π/gr.Ly # fundamental wavenumbers

nlayers = 2 # these choice of parameters give the
f₀, g = 1, 1 # desired PV-streamfunction relations
f₀ = 1 # desired PV-streamfunction relations
H = [0.2, 0.8] # q1 = Δψ1 + 25*(ψ2-ψ1), and
ρ = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2).
b = [-1.0, -1.2] # q2 = Δψ2 + 25/4*(ψ1-ψ2).

β = 0.35

Expand Down Expand Up @@ -471,7 +471,7 @@ function test_mqg_linearadvection(dt, stepper, dev::Device=CPU();
return nothing
end

prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, g, H, ρ, U,
prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, H, b, U,
eta=η, β, μ, ν, nν, calcFq=calcFq!, stepper, dt, linear=true)

sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid
Expand Down Expand Up @@ -509,11 +509,11 @@ function test_mqg_energies(dev::Device=CPU();
k₀, l₀ = 2π/gr.Lx, 2π/gr.Ly # fundamental wavenumbers
nlayers = 2

f₀, g = 1, 1
H = [2.5, 7.5] # sum(params.H) = 10
ρ = [1.0, 2.0] # Make g′ = 1/2
f₀ = 1
H = [2.5, 7.5] # sum(params.H) = 10
b = [-1.0, -1.5] # Make g′ = 1/2

prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, g, H, ρ)
prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, H, b)
sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid

ψ = zeros(dev, eltype(gr), (gr.nx, gr.ny, nlayers))
Expand Down Expand Up @@ -573,13 +573,13 @@ function test_mqg_fluxes(dev::Device=CPU(); dt=0.001, stepper="ForwardEuler", n=
k₀, l₀ = 2π/gr.Lx, 2π/gr.Ly # fundamental wavenumbers

nlayers = 2 # these choice of parameters give the
f₀, g = 1, 1 # desired PV-streamfunction relations
f₀ = 1 # desired PV-streamfunction relations
H = [0.2, 0.8] # q1 = Δψ1 + 25*(ψ2-ψ1), and
ρ = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2).
b = [-1.0, -1.2] # q2 = Δψ2 + 25/4*(ψ1-ψ2).
U = zeros(ny, nlayers)
U[:, 1] = @. sech(gr.y / 0.2)^2

prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, g, H, ρ, U)
prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, H, b, U)

sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid

Expand Down Expand Up @@ -645,11 +645,11 @@ function test_mqg_setqsetψ(dev::Device=CPU(); dt=0.001, stepper="ForwardEuler",
k₀, l₀ = 2π/gr.Lx, 2π/gr.Ly # fundamental wavenumbers

nlayers = 2 # these choice of parameters give the
f₀, g = 1, 1 # desired PV-streamfunction relations
f₀ = 1 # desired PV-streamfunction relations
H = [0.2, 0.8] # q1 = Δψ1 + 25*(ψ2-ψ1), and
ρ = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2).
b = [-1.0, -1.2] # q2 = Δψ2 + 25/4*(ψ1-ψ2).

prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx=L, f₀, g, H, ρ)
prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx=L, f₀, H, b)

sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid

Expand Down Expand Up @@ -692,9 +692,9 @@ function test_mqg_set_topographicPV_largescale_gradient(dev::Device=CPU(); dt=0.
k₀, l₀ = 2π/gr.Lx, 2π/gr.Ly # fundamental wavenumbers

nlayers = 2 # these choice of parameters give the
f₀, g = 1, 1 # desired PV-streamfunction relations
f₀ = 1 # desired PV-streamfunction relations
H = [0.2, 0.8] # q1 = Δψ1 + 25 * (ψ2 - ψ1), and
ρ = [4.0, 5.0] # q2 = Δψ2 + 25/4 * (ψ1 - ψ2).
b = [-1.0, -1.2] # q2 = Δψ2 + 25/4 * (ψ1 - ψ2).
U = zeros(nlayers) # the imposed mean zonal flow in each layer
U[1] = 1.0
U[2] = 0.0
Expand All @@ -712,13 +712,13 @@ function test_mqg_set_topographicPV_largescale_gradient(dev::Device=CPU(); dt=0.
eta = f₀ * h / H[2]

prob = MultiLayerQG.Problem(nlayers, dev;
nx, ny, Lx=L, Ly=L, f₀, g, H, ρ, U, μ, β=0, T, eta,
nx, ny, Lx=L, Ly=L, f₀, H, b, U, μ, β=0, T, eta,
topographic_pv_gradient = f₀ / H[2] .* (topographic_slope_x, topographic_slope_y),
dt, stepper, aliased_fraction=0)

# Test to see if the internally-computed total bottom Qx and Qy are correct

g′ = g * (ρ[2] - ρ[1]) / ρ[2]
g′ = -(b[2] - b[1])
F1 = f₀^2 / (H[2] * g′)
Psi1y, Psi2y = -U[1], -U[2]

Expand All @@ -742,9 +742,9 @@ function test_mqg_paramsconstructor(dev::Device=CPU(); dt=0.001, stepper="Forwar
gr = TwoDGrid(dev; nx, Lx=L, ny, Ly=L)

nlayers = 2 # these choice of parameters give the
f₀, g = 1, 1 # desired PV-streamfunction relations
f₀ = 1 # desired PV-streamfunction relations
H = [0.2, 0.8] # q1 = Δψ1 + 25*(ψ2-ψ1), and
ρ = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2).
b = [-1.0, -1.2] # q2 = Δψ2 + 25/4*(ψ1-ψ2).

U1, U2 = 0.1, 0.05

Expand All @@ -758,8 +758,8 @@ function test_mqg_paramsconstructor(dev::Device=CPU(); dt=0.001, stepper="Forwar
CUDA.@allowscalar Ufloats[1] = U1
CUDA.@allowscalar Ufloats[2] = U2

probUvectors = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx=L, f₀, g, H, ρ, U=Uvectors)
probUfloats = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx=L, f₀, g, H, ρ, U=Ufloats)
probUvectors = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx=L, f₀, H, b, U=Uvectors)
probUfloats = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx=L, f₀, H, b, U=Ufloats)

return isapprox(probUfloats.params.U, probUvectors.params.U, rtol=rtol_multilayerqg)
end
Expand Down

0 comments on commit 6d01319

Please sign in to comment.