diff --git a/Project.toml b/Project.toml index 7175f83d..96b3344f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,8 +1,8 @@ name = "GeophysicalFlows" uuid = "44ee3b1c-bc02-53fa-8355-8e347616e15e" license = "MIT" -authors = ["Navid C. Constantinou ", "Gregory L. Wagner ", "and co-contributors"] -version = "0.15.4" +authors = ["Navid C. Constantinou ", "Gregory L. Wagner ", "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 index dc95d185..922fd2a2 100644 --- 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 index 9669fbbc..96d0488b 100644 --- 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 index 1caf0af4..f2ae13c1 100644 --- 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 index d25577bc..4e35d679 100644 --- 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 diff --git a/test/test_multilayerqg.jl b/test/test_multilayerqg.jl index 3a0df396..ab3e2892 100644 --- a/test/test_multilayerqg.jl +++ b/test/test_multilayerqg.jl @@ -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) @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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)) @@ -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 @@ -645,18 +645,18 @@ 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 T = eltype(gr) f1 = @. 2cos(k₀*x) * cos(l₀*y) - f2 = @. cos(k₀*x+π/10) * cos(2l₀*y) + f2 = @. cos(k₀*x + π/10) * cos(2l₀*y) f = zeros(dev, T, (gr.nx, gr.ny, nlayers)) view(f, :, :, 1) .= f1 @@ -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 @@ -712,13 +712,12 @@ 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[1] - b[2] F1 = f₀^2 / (H[2] * g′) Psi1y, Psi2y = -U[1], -U[2] @@ -742,9 +741,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 @@ -758,8 +757,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