diff --git a/Project.toml b/Project.toml index 3c607655..9ce4d9aa 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "GeophysicalFlows" uuid = "44ee3b1c-bc02-53fa-8355-8e347616e15e" license = "MIT" authors = ["Navid C. Constantinou ", "Gregory L. Wagner ", "and co-contributors"] -version = "0.15.0" +version = "0.15.1" [deps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" diff --git a/src/multilayerqg.jl b/src/multilayerqg.jl index 94a8889e..31924cd3 100644 --- a/src/multilayerqg.jl +++ b/src/multilayerqg.jl @@ -28,28 +28,29 @@ nothingfunction(args...) = nothing """ Problem(nlayers :: Int, - dev = CPU(); - nx = 128, - ny = nx, - Lx = 2π, - Ly = Lx, - f₀ = 1.0, - β = 0.0, - g = 1.0, - U = zeros(nlayers), - H = 1/nlayers * ones(nlayers), - ρ = Array{Float64}(1:nlayers), - eta = nothing, - μ = 0, - ν = 0, - nν = 1, - dt = 0.01, - stepper = "RK4", - calcFq = nothingfunction, - stochastic = false, - linear = false, - aliased_fraction = 1/3, - T = Float64) + dev = CPU(); + nx = 128, + ny = nx, + Lx = 2π, + Ly = Lx, + f₀ = 1.0, + β = 0.0, + g = 1.0, + U = zeros(nlayers), + H = 1/nlayers * ones(nlayers), + ρ = Array{Float64}(1:nlayers), + eta = nothing, + topographic_pv_gradient = (0, 0), + μ = 0, + ν = 0, + nν = 1, + dt = 0.01, + stepper = "RK4", + calcFq = nothingfunction, + stochastic = false, + linear = false, + aliased_fraction = 1/3, + T = Float64) Construct a multi-layer quasi-geostrophic problem with `nlayers` fluid layers on device `dev`. @@ -70,7 +71,8 @@ Keyword arguments - `U`: The imposed constant zonal flow ``U(y)`` in each fluid layer. - `H`: Rest height of each fluid layer. - `ρ`: Density of each fluid layer. - - `eta`: Topographic potential vorticity. + - `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. - `ν`: Small-scale (hyper)-viscosity coefficient. - `nν`: (Hyper)-viscosity order, `nν```≥ 1``. @@ -80,44 +82,45 @@ Keyword arguments - `stochastic`: `true` or `false` (default); boolean denoting whether `calcF` is temporally stochastic. - `linear`: `true` or `false` (default); boolean denoting whether the linearized equations of motions are used. - `aliased_fraction`: the fraction of high-wavenumbers that are zero-ed out by `dealias!()`. - - `T`: `Float32` or `Float64`; floating point type used for `problem` data. + - `T`: `Float32` or `Float64` (default); floating point type used for `problem` data. """ -function Problem(nlayers::Int, # number of fluid layers - dev = CPU(); +function Problem(nlayers::Int, # number of fluid layers + dev = CPU(); # Numerical parameters - nx = 128, - ny = nx, - Lx = 2π, - Ly = Lx, + nx = 128, + ny = nx, + 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, # topographic PV + 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 # Bottom Drag and/or (hyper)-viscosity - μ = 0, - ν = 0, - nν = 1, + μ = 0, + ν = 0, + nν = 1, # Timestepper and equation options - dt = 0.01, - stepper = "RK4", - calcFq = nothingfunction, - stochastic = false, - linear = false, + dt = 0.01, + stepper = "RK4", + calcFq = nothingfunction, + stochastic = false, + linear = false, # Float type and dealiasing - aliased_fraction = 1/3, - T = Float64) + aliased_fraction = 1/3, + T = Float64) if dev == GPU() && nlayers > 2 @warn """MultiLayerQG module is not optimized on the GPU yet for configurations with 3 fluid layers or more! - + See issues on Github at https://github.com/FourierFlows/GeophysicalFlows.jl/issues/112 and https://github.com/FourierFlows/GeophysicalFlows.jl/issues/267. - + To use MultiLayerQG with 3 fluid layers or more we suggest, for now, to restrict running on CPU.""" end @@ -134,7 +137,7 @@ function Problem(nlayers::Int, # number of fluid layers grid = TwoDGrid(dev; nx, Lx, ny, Ly, aliased_fraction, T) - params = Params(nlayers, g, f₀, β, ρ, H, U, eta, μ, ν, nν, grid; calcFq) + params = Params(nlayers, g, f₀, β, ρ, H, U, eta, topographic_pv_gradient, μ, ν, nν, grid; calcFq) vars = calcFq == nothingfunction ? DecayingVars(grid, params) : (stochastic ? StochasticForcedVars(grid, params) : ForcedVars(grid, params)) @@ -168,6 +171,8 @@ struct Params{T, Aphys3D, Aphys2D, Aphys1D, Atrans4D, Trfft} <: AbstractParams U :: Aphys3D "array containing the topographic PV" eta :: Aphys2D + "tuple containing the ``(x, y)`` components of topographic PV large-scale gradient" + topographic_pv_gradient :: Tuple{T, T} "linear bottom drag coefficient" μ :: T "small-scale (hyper)-viscosity coefficient" @@ -205,8 +210,10 @@ struct SingleLayerParams{T, Aphys3D, Aphys2D, Trfft} <: AbstractParams β :: T "array with imposed constant zonal flow ``U(y)``" U :: Aphys3D - "array containing topographic PV" + "array containing the periodic component of the topographic PV" eta :: Aphys2D + "tuple containing the ``(x, y)`` components of topographic PV large-scale gradient" + topographic_pv_gradient :: Tuple{T, T} "linear drag coefficient" μ :: T "small-scale (hyper)-viscosity coefficient" @@ -246,8 +253,10 @@ struct TwoLayerParams{T, Aphys3D, Aphys2D, Trfft} <: AbstractParams H :: Tuple "array with imposed constant zonal flow ``U(y)`` in each fluid layer" U :: Aphys3D - "array containing topographic PV" + "array containing periodic component of the topographic PV" eta :: Aphys2D + "tuple containing the ``(x, y)`` components of topographic PV large-scale gradient" + topographic_pv_gradient :: Tuple{T, T} "linear bottom drag coefficient" μ :: T "small-scale (hyper)-viscosity coefficient" @@ -297,7 +306,7 @@ function convert_U_to_U3D(dev, nlayers, grid, U::Number) return A(U_3D) end -function Params(nlayers, g, f₀, β, ρ, H, U, eta, μ, ν, nν, grid; calcFq=nothingfunction, effort=FFTW.MEASURE) +function Params(nlayers, g, f₀, β, ρ, H, U, eta, topographic_pv_gradient, μ, ν, nν, grid; calcFq=nothingfunction, effort=FFTW.MEASURE) dev = grid.device T = eltype(grid) A = device_array(dev) @@ -311,9 +320,15 @@ function Params(nlayers, g, f₀, β, ρ, H, U, eta, μ, ν, nν, grid; calcFq=n Uyy = real.(ifft(-l.^2 .* fft(U))) Uyy = CUDA.@allowscalar repeat(Uyy, outer=(nx, 1, 1)) + # Calculate periodic components of the topographic PV gradients. etah = rfft(A(eta)) - etax = irfft(im * kr .* etah, nx) - etay = irfft(im * l .* etah, nx) + etax = irfft(im * kr .* etah, nx) # ∂η/∂x + etay = irfft(im * l .* etah, nx) # ∂η/∂y + + # Add topographic PV large-scale gradient + topographic_pv_gradient = T.(topographic_pv_gradient) + @. etax += topographic_pv_gradient[1] + @. etay += topographic_pv_gradient[2] Qx = zeros(dev, T, (nx, ny, nlayers)) @views @. Qx[:, :, nlayers] += etax @@ -325,7 +340,7 @@ function Params(nlayers, g, f₀, β, ρ, H, U, eta, μ, ν, nν, grid; calcFq=n rfftplanlayered = plan_flows_rfft(A{T, 3}(undef, grid.nx, grid.ny, nlayers), [1, 2]; flags=effort) if nlayers==1 - return SingleLayerParams(T(β), U, eta, T(μ), T(ν), nν, calcFq, Qx, Qy, rfftplanlayered) + return SingleLayerParams(T(β), U, eta, topographic_pv_gradient, T(μ), T(ν), nν, calcFq, Qx, Qy, rfftplanlayered) else # if nlayers≥2 @@ -354,9 +369,9 @@ function Params(nlayers, g, f₀, β, ρ, H, U, eta, μ, ν, nν, grid; calcFq=n 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(β), A(ρ), (T(H[1]), T(H[2])), U, eta, T(μ), T(ν), nν, calcFq, T(g′[1]), Qx, Qy, rfftplanlayered) + return TwoLayerParams(T(g), T(f₀), T(β), A(ρ), (T(H[1]), T(H[2])), 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(β), A(ρ), A(H), U, eta, T(μ), T(ν), nν, calcFq, A(g′), Qx, Qy, S, S⁻¹, rfftplanlayered) + return Params(nlayers, T(g), T(f₀), T(β), A(ρ), A(H), U, eta, topographic_pv_gradient, T(μ), T(ν), nν, calcFq, A(g′), Qx, Qy, S, S⁻¹, rfftplanlayered) end end end @@ -615,7 +630,7 @@ case of a two fluid layer configuration. In this case we have, ``` where ``Δ = k² [k² + f₀² (H₁ + H₂) / (g′ H₁ H₂)]``. - + (Here, the PV-streamfunction relationship is hard-coded to avoid scalar operations on the GPU.) """ @@ -680,7 +695,7 @@ end """ calcN!(N, sol, t, clock, vars, params, grid) - + Compute the nonlinear term, that is the advection term, the bottom drag, and the forcing: ```math @@ -704,7 +719,7 @@ end """ calcNlinear!(N, sol, t, clock, vars, params, grid) - + Compute the nonlinear term of the linearized equations: ```math @@ -846,7 +861,7 @@ Update all problem variables using `sol`. """ function updatevars!(vars, params, grid, sol) dealias!(sol, grid) - + @. vars.qh = sol streamfunctionfrompv!(vars.ψh, vars.qh, params, grid) @. vars.uh = -im * grid.l * vars.ψh @@ -895,7 +910,7 @@ set_q!(prob, q) = set_q!(prob.sol, prob.params, prob.vars, prob.grid, q) set_ψ!(params, vars, grid, sol, ψ) set_ψ!(prob, ψ) -Set the solution `prob.sol` to the transform `qh` that corresponds to streamfunction `ψ` +Set the solution `prob.sol` to the transform `qh` that corresponds to streamfunction `ψ` and update variables. """ function set_ψ!(sol, params, vars, grid, ψ) @@ -936,7 +951,7 @@ The kinetic energy at the ``j``-th fluid layer is 𝖪𝖤_j = \\frac{H_j}{H} \\int \\frac1{2} |{\\bf ∇} ψ_j|^2 \\frac{𝖽x 𝖽y}{L_x L_y} = \\frac1{2} \\frac{H_j}{H} \\sum_{𝐤} |𝐤|² |ψ̂_j|², \\ j = 1, ..., n , ``` -while the potential energy that corresponds to the interface ``j+1/2`` (i.e., the interface +while the potential energy that corresponds to the interface ``j+1/2`` (i.e., the interface between the ``j``-th and ``(j+1)``-th fluid layer) is ```math @@ -949,7 +964,7 @@ function energies(vars, params, grid, sol) @. vars.qh = sol streamfunctionfrompv!(vars.ψh, vars.qh, params, grid) - + abs²∇𝐮h = vars.uh # use vars.uh as scratch variable @. abs²∇𝐮h = grid.Krsq * abs2(vars.ψh) @@ -1002,7 +1017,7 @@ energies(prob) = energies(prob.vars, prob.params, prob.grid, prob.sol) fluxes(prob) Return the lateral eddy fluxes within each fluid layer, lateralfluxes``_1,...,``lateralfluxes``_n`` -and also the vertical eddy fluxes at each fluid interface, +and also the vertical eddy fluxes at each fluid interface, verticalfluxes``_{3/2},...,``verticalfluxes``_{n-1/2}``, where ``n`` is the total number of layers in the fluid. (When ``n=1``, only the lateral fluxes are returned.) @@ -1017,7 +1032,7 @@ while the vertical eddy fluxes at the ``j+1/2``-th fluid interface (i.e., interf the ``j``-th and ``(j+1)``-th fluid layer) are ```math -\\textrm{verticalfluxes}_{j+1/2} = \\int \\frac{f₀²}{g'_{j+1/2} H} (U_j - U_{j+1}) \\, +\\textrm{verticalfluxes}_{j+1/2} = \\int \\frac{f₀²}{g'_{j+1/2} H} (U_j - U_{j+1}) \\, v_{j+1} ψ_{j} \\frac{𝖽x 𝖽y}{L_x L_y} , \\ j = 1, ..., n-1. ``` """ @@ -1047,7 +1062,7 @@ end function fluxes(vars, params::TwoLayerParams, grid, sol) nlayers = numberoflayers(params) - + lateralfluxes, verticalfluxes = zeros(nlayers), zeros(nlayers-1) updatevars!(vars, params, grid, sol) diff --git a/test/runtests.jl b/test/runtests.jl index 43c7fbbc..0664b009 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,25 +18,25 @@ using GeophysicalFlows: lambdipole, peakedisotropicspectrum # the devices on which tests will run devices = CUDA.functional() ? (CPU(), GPU()) : (CPU(),) -const rtol_lambdipole = 1e-2 # tolerance for lamb dipole tests +const rtol_lambdipole = 1e-2 # tolerance for lamb dipole tests const rtol_twodnavierstokes = 1e-13 # tolerance for twodnavierstokes forcing tests -const rtol_singlelayerqg = 1e-13 # tolerance for singlelayerqg forcing tests -const rtol_multilayerqg = 1e-13 # tolerance for multilayerqg forcing tests -const rtol_surfaceqg = 1e-13 # tolerance for surfaceqg forcing tests +const rtol_singlelayerqg = 1e-13 # tolerance for singlelayerqg forcing tests +const rtol_multilayerqg = 1e-13 # tolerance for multilayerqg forcing tests +const rtol_surfaceqg = 1e-13 # tolerance for surfaceqg forcing tests # Run tests testtime = @elapsed begin for dev in devices - + @info "testing on " * string(typeof(dev)) - + @testset "Utils" begin include("test_utils.jl") @test testpeakedisotropicspectrum(dev) @test_throws ErrorException("the domain is not square") testpeakedisotropicspectrum_rectangledomain() end - + @testset "TwoDNavierStokes" begin include("test_twodnavierstokes.jl") @@ -50,10 +50,10 @@ for dev in devices @test test_twodnavierstokes_problemtype(dev, Float32) @test TwoDNavierStokes.nothingfunction() == nothing end - + @testset "SingleLayerQG" begin include("test_singlelayerqg.jl") - + for deformation_radius in [Inf, 1.23] @test test_1layerqg_rossbywave("ETDRK4", 1e-2, 20, dev, deformation_radius=deformation_radius) @test test_1layerqg_rossbywave("FilteredETDRK4", 1e-2, 20, dev, deformation_radius=deformation_radius) @@ -81,7 +81,7 @@ for dev in devices @test_throws ErrorException("not implemented for finite deformation radius") test_1layerqg_energy_drag(dev; deformation_radius=2.23) @test_throws ErrorException("not implemented for finite deformation radius") test_1layerqg_enstrophy_drag(dev; deformation_radius=2.23) end - + @testset "BarotropicQGQL" begin include("test_barotropicqgql.jl") @@ -100,10 +100,10 @@ for dev in devices @test test_bqgql_problemtype(dev, Float32) @test BarotropicQGQL.nothingfunction() == nothing end - + @testset "SurfaceQG" begin include("test_surfaceqg.jl") - + @test test_sqg_kineticenergy_buoyancyvariance(dev) @test test_sqg_advection(0.0005, "ForwardEuler", dev) @test test_sqg_deterministicforcing_buoyancy_variance_budget(dev) @@ -114,10 +114,10 @@ for dev in devices @test test_sqg_noforcing(dev) @test SurfaceQG.nothingfunction() == nothing end - + @testset "MultiLayerQG" begin include("test_multilayerqg.jl") - + @test test_pvtofromstreamfunction_2layer(dev) @test test_pvtofromstreamfunction_3layer(dev) @test test_mqg_rossbywave("RK4", 1e-2, 20, dev) @@ -128,6 +128,7 @@ for dev in devices @test test_mqg_fluxes(dev) @test test_mqg_fluxessinglelayer(dev) @test test_mqg_setqsetψ(dev) + @test test_mqg_set_topographicPV_largescale_gradient(dev) @test test_mqg_paramsconstructor(dev) @test test_mqg_stochasticforcedproblemconstructor(dev) @test test_mqg_problemtype(dev, Float32) diff --git a/test/test_multilayerqg.jl b/test/test_multilayerqg.jl index db56f20a..9e8f717c 100644 --- a/test/test_multilayerqg.jl +++ b/test/test_multilayerqg.jl @@ -48,7 +48,7 @@ function constructtestfields_3layer(gr) ψ1 = @. 1e-3 * ( 1/4*cos(2k₀*x)*cos(5l₀*y) + 1/3*cos(3k₀*x)*cos(3l₀*y) ) ψ2 = @. 1e-3 * ( cos(3k₀*x)*cos(4l₀*y) + 1/2*cos(4k₀*x)*cos(2l₀*y) ) ψ3 = @. 1e-3 * ( cos(1k₀*x)*cos(3l₀*y) + 1/2*cos(2k₀*x)*cos(2l₀*y) ) - + Δψ1 = @. -1e-3 * ( 1/4*((2k₀)^2+(5l₀)^2)*cos(2k₀*x)*cos(5l₀*y) + 1/3*((3k₀)^2+(3l₀)^2)*cos(3k₀*x)*cos(3l₀*y) ) Δψ2 = @. -1e-3 * ( ((3k₀)^2+(4l₀)^2)*cos(3k₀*x)*cos(4l₀*y) + 1/2*((4k₀)^2+(2l₀)^2)*cos(4k₀*x)*cos(2l₀*y) ) Δψ3 = @. -1e-3 * ( ((1k₀)^2+(3l₀)^2)*cos(1k₀*x)*cos(3l₀*y) + 1/2*((2k₀)^2+(2l₀)^2)*cos(2k₀*x)*cos(2l₀*y) ) @@ -111,7 +111,7 @@ Tests the pvfromstreamfunction function that gives qh from ψh. To do so, it constructs a 3-layer problem with parameters such that q1 = Δψ1 + 20ψ2 - 20ψ1, q2 = Δψ2 + 20ψ1 - 44ψ2 + 24ψ3, q3 = Δψ3 + 12ψ2 - 12ψ3. Then given ψ1, ψ2, and ψ2 it tests if pvfromstreamfunction gives the expected -q1, q2, and q3. Similarly, that streamfunctionfrompv gives ψ1, ψ2, and ψ2 from +q1, q2, and q3. Similarly, that streamfunctionfrompv gives ψ1, ψ2, and ψ2 from q1, q2, and q3. """ function test_pvtofromstreamfunction_3layer(dev::Device=CPU()) @@ -164,7 +164,7 @@ that a solution to the problem forced by this Ff is then qf. """ function test_mqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); n=128, L=2π, nlayers=2, μ=0.0, ν=0.0, nν=1) - + A = device_array(dev) tf = 0.5 @@ -181,9 +181,9 @@ function test_mqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); f₀, g = 1, 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). - + β = 0.35 - + U1, U2 = 0.1, 0.05 u1 = @. 0.5sech(gr.y/(Ly/15))^2; u1 = A(reshape(u1, (1, gr.ny))) u2 = @. 0.02cos(3l₀*gr.y); u2 = A(reshape(u2, (1, gr.ny))) @@ -234,12 +234,12 @@ function test_mqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); @views ψf[:, :, 2] = ψ2 MultiLayerQG.set_q!(prob, qf) - + stepforward!(prob, nt) MultiLayerQG.updatevars!(prob) - + return isapprox(vs.q, qf, rtol=rtol_multilayerqg) && - isapprox(vs.ψ, ψf, rtol=rtol_multilayerqg) + isapprox(vs.ψ, ψf, rtol=rtol_multilayerqg) end """ @@ -255,9 +255,9 @@ is unstable.) """ function test_mqg_linearadvection(dt, stepper, dev::Device=CPU(); n=128, L=2π, nlayers=2, μ=0.0, ν=0.0, nν=1) - + A = device_array(dev) - + tf = 0.5 nt = round(Int, tf/dt) @@ -295,9 +295,9 @@ function test_mqg_linearadvection(dt, stepper, dev::Device=CPU(); Ff1 = (β .- uyy1 .- 25*(U2.+u2.-U1.-u1) ).*ψ1x + (U1.+u1).*q1x - ν*Δq1 Ff2 = FourierFlows.jacobian(ψ2, η, gr) + (β .- uyy2 .- 25/4*(U1.+u1.-U2.-u2) ).*ψ2x + (U2.+u2).*(q2x + ηx) + μ*Δψ2 - ν*Δq2 - + T = eltype(gr) - + Ff = zeros(dev, T, (gr.nx, gr.ny, nlayers)) @views Ff[:, :, 1] = Ff1 @views Ff[:, :, 2] = Ff2 @@ -313,7 +313,7 @@ function test_mqg_linearadvection(dt, stepper, dev::Device=CPU(); prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, g, H, ρ, U, eta=η, β, μ, ν, nν, calcFq=calcFq!, stepper, dt, linear=true) - + sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid qf = zeros(dev, T, (gr.nx, gr.ny, nlayers)) @@ -348,16 +348,16 @@ function test_mqg_energies(dev::Device=CPU(); x, y = gridpoints(gr) 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 - + prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, g, H, ρ) sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid ψ = zeros(dev, eltype(gr), (gr.nx, gr.ny, nlayers)) - + CUDA.@allowscalar @views @. ψ[:, :, 1] = cos(2k₀ * x) * cos(2l₀ * y) CUDA.@allowscalar @views @. ψ[:, :, 2] = 1/2 * cos(2k₀ * x) * cos(2l₀ * y) @@ -380,10 +380,10 @@ function test_mqg_energysinglelayer(dev::Device=CPU(); nx, Lx = 64, 2π ny, Ly = 64, 3π gr = TwoDGrid(dev; nx, Lx, ny, Ly) - + x, y = gridpoints(gr) k₀, l₀ = 2π/gr.Lx, 2π/gr.Ly # fundamental wavenumbers - + energy_calc = 29/9 ψ0 = @. sin(2k₀*x)*cos(2l₀*y) + 2sin(k₀*x)*cos(3l₀*y) @@ -439,12 +439,12 @@ end """ test_mqg_fluxessinglelayer(dt, stepper; kwargs...) -Tests the lateral eddy fluxes by constructing a 1-layer problem and initializing +Tests the lateral eddy fluxes by constructing a 1-layer problem and initializing it with a flow field whose fluxes are known. """ -function test_mqg_fluxessinglelayer(dev::Device=CPU(); dt=0.001, stepper="ForwardEuler", n=128, L=2π, μ=0.0, ν=0.0, nν=1) +function test_mqg_fluxessinglelayer(dev::Device=CPU(); dt=0.001, stepper="ForwardEuler", n=128, L=2π, μ=0.0, ν=0.0, nν=1) nlayers = 1 - + nx, ny = 128, 126 Lx, Ly = 2π, 2π gr = TwoDGrid(dev; nx, Lx, ny, Ly) @@ -486,11 +486,11 @@ function test_mqg_setqsetψ(dev::Device=CPU(); dt=0.001, stepper="ForwardEuler", ρ = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2). prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx=L, f₀, g, H, ρ) - + 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) f = zeros(dev, T, (gr.nx, gr.ny, nlayers)) @@ -512,6 +512,59 @@ function test_mqg_setqsetψ(dev::Device=CPU(); dt=0.001, stepper="ForwardEuler", isapprox(qtest, f, rtol=rtol_multilayerqg) end +""" + test_mqg_set_topographicPV_largescale_gradient(dt, stepper; kwargs...) + +Ensures that the error of setting a topographic PV large-scale gradient is small. +""" +function test_mqg_set_topographicPV_largescale_gradient(dev::Device=CPU(); dt=0.001, stepper="ForwardEuler", n=64, L=2π, nlayers=2, μ=0.0, ν=0.0, nν=1) + nx, ny = 32, 34 + L = 2π + gr = TwoDGrid(dev; nx, Lx=L, ny, Ly=L) + T = eltype(gr) + + x, y = gridpoints(gr) + 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 + H = [0.2, 0.8] # q1 = Δψ1 + 25 * (ψ2 - ψ1), and + ρ = [4.0, 5.0] # 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 + + # topographic amplitude + h₀ = 0.1 + + # topographic amplitude + topographic_slope_x, topographic_slope_y = 1e-2, 1e-1 + + # topography (bumps) + h = @. h₀ * cos(k₀*x) * cos(l₀*y) + + # topographic PV + eta = f₀ * h / H[2] + + prob = MultiLayerQG.Problem(nlayers, dev; + nx, ny, Lx=L, Ly=L, f₀, g, H, ρ, 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] + F1 = f₀^2 / (H[2] * g′) + Psi1y, Psi2y = -U[1], -U[2] + + Q2x_analytic = @. f₀ * (topographic_slope_x - h₀ * k₀ * sin(k₀*x) * cos(l₀*y)) / H[2] + Q2y_analytic = @. f₀ * (topographic_slope_y - h₀ * l₀ * cos(k₀*x) * sin(l₀*y)) / H[2] + F1 * (Psi1y - Psi2y) + + return prob.params.topographic_pv_gradient == f₀ / H[2] .* (topographic_slope_x, topographic_slope_y) && + isapprox(prob.params.Qx[:, :, 2], Q2x_analytic, rtol=rtol_multilayerqg) && + isapprox(prob.params.Qy[:, :, 2], Q2y_analytic, rtol=rtol_multilayerqg) +end + """ test_paramsconstructor(; kwargs...) @@ -527,11 +580,11 @@ function test_mqg_paramsconstructor(dev::Device=CPU(); dt=0.001, stepper="Forwar f₀, g = 1, 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). - + U1, U2 = 0.1, 0.05 T = eltype(gr) - + Uvectors = zeros(dev, T, (ny, nlayers)) Uvectors[:, 1] .= U1 Uvectors[:, 2] .= U2 @@ -549,9 +602,9 @@ end function test_mqg_problemtype(dev, T) prob1 = MultiLayerQG.Problem(1, dev; T) prob2 = MultiLayerQG.Problem(2, dev; T) - + A = device_array(dev) - + return typeof(prob1.sol)<:A{Complex{T}, 3} && typeof(prob1.grid.Lx)==T && typeof(prob1.vars.u)<:A{T, 3} && @@ -562,7 +615,7 @@ end """ test_mqg_rossbywave(; kwargs...) -Evolves a Rossby wave on a beta plane with an imposed zonal flow U and compares +Evolves a Rossby wave on a beta plane with an imposed zonal flow U and compares with the analytic solution. """ function test_mqg_rossbywave(stepper, dt, nsteps, dev::Device=CPU()) @@ -599,19 +652,19 @@ end function test_numberoflayers(dev::Device=CPU()) prob_nlayers1 = MultiLayerQG.Problem(1, dev) prob_nlayers2 = MultiLayerQG.Problem(2, dev) - + return MultiLayerQG.numberoflayers(prob_nlayers1)==1 && MultiLayerQG.numberoflayers(prob_nlayers2)==2 end function test_mqg_stochasticforcedproblemconstructor(dev::Device=CPU()) - + function calcFq!(Fqh, sol, t, clock, vars, params, grid) Fqh .= Ffh return nothing end - + prob = MultiLayerQG.Problem(2, dev; calcFq=calcFq!, stochastic=true) - + return typeof(prob.vars.prevsol)==typeof(prob.sol) end