From 00a839597d61ae7091f5b0c9be0655110e360940 Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Tue, 10 Sep 2019 17:50:53 +0300 Subject: [PATCH 1/3] =?UTF-8?q?changes=20beta->=CE=B2=20and=20rho->=CF=81?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/multilayerqg.jl | 50 ++++++++++---------- test/test_multilayerqg.jl | 96 +++++++++++++++++++-------------------- 2 files changed, 73 insertions(+), 73 deletions(-) diff --git a/src/multilayerqg.jl b/src/multilayerqg.jl index 15eb29b1..4426cac9 100644 --- a/src/multilayerqg.jl +++ b/src/multilayerqg.jl @@ -42,11 +42,11 @@ function Problem(; # Physical parameters nlayers = 2, # number of fluid layers f0 = 1.0, # Coriolis parameter - beta = 0.0, # y-gradient of Coriolis parameter + β = 0.0, # y-gradient of Coriolis parameter g = 1.0, # gravitational constant U = zeros(ny, nlayers), # imposed zonal flow U(y) in each layer H = [0.2, 0.8], # rest fluid height of each layer - rho = [4.0, 5.0], # density of each layer + ρ = [4.0, 5.0], # density of each layer eta = zeros(nx, ny), # topographic PV # Bottom Drag and/or (hyper)-viscosity μ = 0.0, @@ -59,7 +59,7 @@ function Problem(; T = Float64) grid = TwoDGrid(nx, Lx, ny, Ly; T=T) - params = Params(nlayers, g, f0, beta, Array{T}(rho), Array{T}(H), Array{T}(U), Array{T}(eta), μ, ν, nν, grid, calcFq=calcFq) + params = Params(nlayers, g, f0, β, Array{T}(ρ), Array{T}(H), Array{T}(U), Array{T}(eta), μ, ν, nν, grid, calcFq=calcFq) vars = calcFq == nothingfunction ? Vars(grid, params) : ForcedVars(grid, params) eqn = linear ? LinearEquation(params, grid) : Equation(params, grid) @@ -72,9 +72,9 @@ struct Params{T} <: AbstractParams # prescribed params nlayers :: Int # Number of fluid layers g :: T # Gravitational constant - f0:: T # Constant planetary vorticity - beta :: T # Planetary vorticity y-gradient - rho :: Array{T,3} # Array with density of each fluid layer + f0 :: T # Constant planetary vorticity + β :: T # Planetary vorticity y-gradient + ρ :: Array{T,3} # Array with density of each fluid layer H :: Array{T,3} # Array with rest height of each fluid layer U :: Array{T,3} # Array with imposed constant zonal flow U(y) in each fluid layer eta :: Array{T,2} # Array containing topographic PV @@ -85,8 +85,8 @@ struct Params{T} <: AbstractParams # derived params g′ :: Array{T,1} # Array with the reduced gravity constants for each fluid interface - Qx :: Array{T,3} # Array containing x-gradient of PV due to U and eta in each fluid layer - Qy :: Array{T,3} # Array containing y-gradient of PV due to beta, U, and eta in each fluid layer + Qx :: Array{T,3} # Array containing x-gradient of PV due to eta in each fluid layer + Qy :: Array{T,3} # Array containing y-gradient of PV due to β, U, and eta in each fluid layer S :: Array{T,4} # Array containing coeffients for getting PV from streamfunction invS :: Array{T,4} # Array containing coeffients for inverting PV to streamfunction rfftplan :: FFTW.rFFTWPlan{T,-1,false,3} # rfft plan for FFTs @@ -96,8 +96,8 @@ struct SingleLayerParams{T} <: BarotropicParams # prescribed params g :: T # Gravitational constant f0 :: T # Constant planetary vorticity - beta :: T # Planetary vorticity y-gradient - rho :: T # Fluid layer density + β :: T # Planetary vorticity y-gradient + ρ :: T # Fluid layer density H :: T # Fluid layer rest height U :: Array{T,1} # Imposed constant zonal flow U(y) eta :: Array{T,2} # Array containing topographic PV @@ -107,12 +107,12 @@ struct SingleLayerParams{T} <: BarotropicParams calcFq! :: Function # Function that calculates the forcing on QGPV q # derived params - Qx :: Array{T,2} # Array containing zonal PV gradient due to beta, U, and eta in each fluid layer - Qy :: Array{T,2} # Array containing meridional PV gradient due to beta, U, and eta in each fluid layer + Qx :: Array{T,2} # Array containing x-gradient of PV due to eta in each fluid layer + Qy :: Array{T,2} # Array containing meridional PV gradient due to β, U, and eta in each fluid layer rfftplan :: FFTW.rFFTWPlan{T,-1,false,2} # rfft plan for FFTs end -function Params(nlayers, g, f0, beta, rho, H, U::Array{T,2}, eta, μ, ν, nν, grid::AbstractGrid{T}; calcFq=nothingfunction, effort=FFTW.MEASURE) where T +function Params(nlayers, g, f0, β, ρ, H, U::Array{T,2}, eta, μ, ν, nν, grid::AbstractGrid{T}; calcFq=nothingfunction, effort=FFTW.MEASURE) where T ny, nx = grid.ny , grid.nx nkr, nl = grid.nkr, grid.nl @@ -120,14 +120,14 @@ function Params(nlayers, g, f0, beta, rho, H, U::Array{T,2}, eta, μ, ν, nν, g U = reshape(U, (1, ny, nlayers)) - # g′ = g*(rho[2:nlayers]-rho[1:nlayers-1]) ./ rho[1:nlayers-1] # definition match PYQG - g′ = g*(rho[2:nlayers]-rho[1:nlayers-1]) ./ rho[2:nlayers] # correct definition + # g′ = g*(ρ[2:nlayers]-ρ[1:nlayers-1]) ./ ρ[1:nlayers-1] # definition match PYQG + g′ = g*(ρ[2:nlayers]-ρ[1:nlayers-1]) ./ ρ[2:nlayers] # correct definition Fm = @. f0^2 / ( g′*H[2:nlayers ] ) Fp = @. f0^2 / ( g′*H[1:nlayers-1] ) - rho = reshape(rho, (1, 1, nlayers)) - H = reshape(H , (1, 1, nlayers)) + ρ = reshape(ρ, (1, 1, nlayers)) + H = reshape(H, (1, 1, nlayers)) Uyy = repeat(irfft( -l.^2 .* rfft(U, [1, 2]), 1, [1, 2]), outer=(1, 1, 1)) Uyy = repeat(Uyy, outer=(nx, 1, 1)) @@ -140,11 +140,11 @@ function Params(nlayers, g, f0, beta, rho, H, U::Array{T,2}, eta, μ, ν, nν, g @views @. Qx[:, :, nlayers] += etax Qy = zeros(nx, ny, nlayers) - Qy[:, :, 1] = @. beta - Uyy[:, :, 1] - Fp[1]*( U[:, :, 2] - U[:, :, 1] ) + Qy[:, :, 1] = @. β - Uyy[:, :, 1] - Fp[1]*( U[:, :, 2] - U[:, :, 1] ) for j = 2:nlayers-1 - Qy[:, :, j] = @. beta - Uyy[:, :, j] - Fp[j]*( U[:, :, j+1] - U[:, :, j] ) - Fm[j-1]*( U[:, :, j-1] - U[:, :, j] ) + Qy[:, :, j] = @. β - Uyy[:, :, j] - Fp[j]*( U[:, :, j+1] - U[:, :, j] ) - Fm[j-1]*( U[:, :, j-1] - U[:, :, j] ) end - @views Qy[:, :, nlayers] = @. beta - Uyy[:, :, nlayers] - Fm[nlayers-1]*( U[:, :, nlayers-1] - U[:, :, nlayers] ) + @views Qy[:, :, nlayers] = @. β - Uyy[:, :, nlayers] - Fm[nlayers-1]*( U[:, :, nlayers-1] - U[:, :, nlayers] ) @views @. Qy[:, :, nlayers] += etay S = Array{T}(undef, (nkr, nl, nlayers, nlayers)) @@ -156,16 +156,16 @@ function Params(nlayers, g, f0, beta, rho, H, U::Array{T,2}, eta, μ, ν, nν, g rfftplanlayered = plan_rfft(Array{T,3}(undef, grid.nx, grid.ny, nlayers), [1, 2]; flags=effort) if nlayers == 1 - SingleLayerParams{T}(g, f0, beta, rho, H, U, eta, μ, ν, nν, calcFq, Qx, Qy, grid.rfftplan) + SingleLayerParams{T}(g, f0, β, ρ, H, U, eta, μ, ν, nν, calcFq, Qx, Qy, grid.rfftplan) else - Params{T}(nlayers, g, f0, beta, rho, H, U, eta, μ, ν, nν, calcFq, g′, Qx, Qy, S, invS, rfftplanlayered) + Params{T}(nlayers, g, f0, β, ρ, H, U, eta, μ, ν, nν, calcFq, g′, Qx, Qy, S, invS, rfftplanlayered) end end -function Params(nlayers, g, f0, beta, rho, H, U::Array{T,1}, eta, μ, ν, nν, grid::AbstractGrid{T}; calcFq=nothingfunction, effort=FFTW.MEASURE) where T +function Params(nlayers, g, f0, β, ρ, H, U::Array{T,1}, eta, μ, ν, nν, grid::AbstractGrid{T}; calcFq=nothingfunction, effort=FFTW.MEASURE) where T U = reshape(U, (1, nlayers)) U = repeat(U, outer=(grid.ny, 1)) - Params(nlayers, g, f0, beta, rho, H, U, eta, μ, ν, nν, grid; calcFq=calcFq, effort=effort) + Params(nlayers, g, f0, β, ρ, H, U, eta, μ, ν, nν, grid; calcFq=calcFq, effort=effort) end @@ -191,7 +191,7 @@ function Equation(p, gr::AbstractGrid{T}) where T end function Equation(p::SingleLayerParams, gr::AbstractGrid{T}) where T - L = @. -p.μ - p.ν*gr.Krsq^p.nν + im*p.beta*gr.kr*gr.invKrsq + L = @. -p.μ - p.ν*gr.Krsq^p.nν + im*p.β*gr.kr*gr.invKrsq L[1, 1] = 0 FourierFlows.Equation(L, calcN!, gr) end diff --git a/test/test_multilayerqg.jl b/test/test_multilayerqg.jl index e481e519..a2b2c996 100644 --- a/test/test_multilayerqg.jl +++ b/test/test_multilayerqg.jl @@ -42,15 +42,15 @@ Then given ψ1 and ψ2 it tests if pvfromstreamfunction gives the expected q1 and q2. Similarly, that streamfunctionfrompv gives ψ1 and ψ2 from q1 and q2. """ function test_pvtofromstreamfunction() - n, L = 128, 2π - gr = TwoDGrid(n, L) + n, L = 128, 2π + gr = TwoDGrid(n, L) - nlayers = 2 # these choice of parameters give the - f0, g = 1, 1 # desired PV-streamfunction relations - H = [0.2, 0.8] # q1 = Δψ1 + 25*(ψ2-ψ1), and - rho = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2). + nlayers = 2 # these choice of parameters give the + f0, 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). - prob = MultilayerQG.Problem(nlayers=nlayers, nx=n, Lx=L, f0=f0, g=g, H=H, rho=rho) + prob = MultilayerQG.Problem(nlayers=nlayers, nx=n, Lx=L, f0=f0, g=g, H=H, ρ=ρ) 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(gr) @@ -92,12 +92,12 @@ function test_mqg_nonlinearadvection(dt, stepper; n=128, L=2π, nlayers=2, μ=0. x, y = gridpoints(gr) k0, l0 = gr.k[2], gr.l[2] # fundamental wavenumbers - nlayers = 2 # these choice of parameters give the - f0, g = 1, 1 # desired PV-streamfunction relations - H = [0.2, 0.8] # q1 = Δψ1 + 25*(ψ2-ψ1), and - rho = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2). + nlayers = 2 # these choice of parameters give the + f0, 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). - beta = 0.35 + β = 0.35 U1, U2 = 0.1, 0.05 u1 = @. 0.5sech(gr.y/(Ly/15))^2 @@ -117,8 +117,8 @@ function test_mqg_nonlinearadvection(dt, stepper; n=128, L=2π, nlayers=2, μ=0. ψ1, ψ2, q1, q2, ψ1x, ψ2x, q1x, q2x, Δψ2, Δq1, Δq2 = constructtestfields(gr) - Ff1 = FourierFlows.jacobian(ψ1, q1, gr) + (beta .- uyy1 .- 25*(U2.+u2.-U1.-u1) ).*ψ1x + (U1.+u1).*q1x - ν*Δq1 - Ff2 = FourierFlows.jacobian(ψ2, q2 + η, gr) + (beta .- uyy2 .- 25/4*(U1.+u1.-U2.-u2) ).*ψ2x + (U2.+u2).*(q2x + ηx) + μ*Δψ2 - ν*Δq2 + Ff1 = FourierFlows.jacobian(ψ1, q1, gr) + (β .- uyy1 .- 25*(U2.+u2.-U1.-u1) ).*ψ1x + (U1.+u1).*q1x - ν*Δq1 + Ff2 = FourierFlows.jacobian(ψ2, q2 + η, gr) + (β .- uyy2 .- 25/4*(U1.+u1.-U2.-u2) ).*ψ2x + (U2.+u2).*(q2x + ηx) + μ*Δψ2 - ν*Δq2 Ff = zeros(gr.nx, gr.ny, nlayers) Ff[:, :, 1] .= Ff1 @@ -134,7 +134,7 @@ function test_mqg_nonlinearadvection(dt, stepper; n=128, L=2π, nlayers=2, μ=0. end prob = MultilayerQG.Problem(nlayers=nlayers, nx=nx, ny=ny, Lx=Lx, Ly=Ly, f0=f0, - g=g, H=H, rho=rho, U=U, eta=η, beta=beta, μ=μ, ν=ν, nν=nν, calcFq=calcFq!) + g=g, H=H, ρ=ρ, U=U, eta=η, β=β, μ=μ, ν=ν, nν=nν, calcFq=calcFq!) sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid qf = zeros(gr.nx, gr.ny, nlayers) @@ -175,12 +175,12 @@ function test_mqg_linearadvection(dt, stepper; n=128, L=2π, nlayers=2, μ=0.0, x, y = gridpoints(gr) k0, l0 = gr.k[2], gr.l[2] # fundamental wavenumbers - nlayers = 2 # these choice of parameters give the - f0, g = 1, 1 # desired PV-streamfunction relations - H = [0.2, 0.8] # q1 = Δψ1 + 25*(ψ2-ψ1), and - rho = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2). + nlayers = 2 # these choice of parameters give the + f0, 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). - beta = 0.35 + β = 0.35 U1, U2 = 0.1, 0.05 u1 = @. 0.5sech(gr.y/(Ly/15))^2 @@ -200,8 +200,8 @@ function test_mqg_linearadvection(dt, stepper; n=128, L=2π, nlayers=2, μ=0.0, ψ1, ψ2, q1, q2, ψ1x, ψ2x, q1x, q2x, Δψ2, Δq1, Δq2 = constructtestfields(gr) - Ff1 = (beta .- uyy1 .- 25*(U2.+u2.-U1.-u1) ).*ψ1x + (U1.+u1).*q1x - ν*Δq1 - Ff2 = FourierFlows.jacobian(ψ2, η, gr) + (beta .- uyy2 .- 25/4*(U1.+u1.-U2.-u2) ).*ψ2x + (U2.+u2).*(q2x + ηx) + μ*Δψ2 - ν*Δq2 + 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 Ff = zeros(gr.nx, gr.ny, nlayers) Ff[:, :, 1] .= Ff1 @@ -217,7 +217,7 @@ function test_mqg_linearadvection(dt, stepper; n=128, L=2π, nlayers=2, μ=0.0, end prob = MultilayerQG.Problem(nlayers=nlayers, nx=nx, ny=ny, Lx=Lx, Ly=Ly, f0=f0, - g=g, H=H, rho=rho, U=U, eta=η, beta=beta, μ=μ, ν=ν, nν=nν, calcFq=calcFq!, linear=true) + g=g, H=H, ρ=ρ, U=U, eta=η, β=β, μ=μ, ν=ν, nν=nν, calcFq=calcFq!, linear=true) sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid @@ -251,12 +251,12 @@ function test_mqg_energies(;dt=0.001, stepper="ForwardEuler", n=128, L=2π, nlay x, y = gridpoints(gr) k0, l0 = gr.k[2], gr.l[2] # fundamental wavenumbers - nlayers = 2 # these choice of parameters give the - f0, g = 1, 1 # desired PV-streamfunction relations - H = [0.2, 0.8] # q1 = Δψ1 + 25*(ψ2-ψ1), and - rho = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2). + nlayers = 2 # these choice of parameters give the + f0, 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). - prob = MultilayerQG.Problem(nlayers=nlayers, nx=nx, ny=ny, Lx=Lx, Ly=Ly, f0=f0, g=g, H=H, rho=rho) + prob = MultilayerQG.Problem(nlayers=nlayers, nx=nx, ny=ny, Lx=Lx, Ly=Ly, f0=f0, g=g, H=H, ρ=ρ) 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(gr) @@ -286,14 +286,14 @@ function test_mqg_fluxes(;dt=0.001, stepper="ForwardEuler", n=128, L=2π, nlayer x, y = gridpoints(gr) k0, l0 = gr.k[2], gr.l[2] # fundamental wavenumbers - nlayers = 2 # these choice of parameters give the - f0, g = 1, 1 # desired PV-streamfunction relations - H = [0.2, 0.8] # q1 = Δψ1 + 25*(ψ2-ψ1), and - rho = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2). - U = zeros(ny, nlayers) - U[:, 1] = @. sech(gr.y/0.2)^2 + nlayers = 2 # these choice of parameters give the + f0, 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(ny, nlayers) + U[:, 1] = @. sech(gr.y/0.2)^2 - prob = MultilayerQG.Problem(nlayers=nlayers, nx=nx, ny=ny, Lx=Lx, Ly=Ly, f0=f0, g=g, H=H, rho=rho, U=U) + prob = MultilayerQG.Problem(nlayers=nlayers, nx=nx, ny=ny, Lx=Lx, Ly=Ly, f0=f0, g=g, H=H, ρ=ρ, U=U) sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid ψ1 = @. cos(k0*x)*cos(l0*y) @@ -321,12 +321,12 @@ function test_mqg_setqsetpsi(;dt=0.001, stepper="ForwardEuler", n=64, L=2π, nla x, y = gridpoints(gr) k0, l0 = gr.k[2], gr.l[2] # fundamental wavenumbers - nlayers = 2 # these choice of parameters give the - f0, g = 1, 1 # desired PV-streamfunction relations - H = [0.2, 0.8] # q1 = Δψ1 + 25*(ψ2-ψ1), and - rho = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2). + nlayers = 2 # these choice of parameters give the + f0, 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). - prob = MultilayerQG.Problem(nlayers=nlayers, nx=nx, ny=ny, Lx=L, f0=f0, g=g, H=H, rho=rho) + prob = MultilayerQG.Problem(nlayers=nlayers, nx=nx, ny=ny, Lx=L, f0=f0, g=g, H=H, ρ=ρ) sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid f1 = @. 2cos(k0*x)*cos(l0*y) @@ -360,12 +360,12 @@ function test_mqg_paramsconstructor(;dt=0.001, stepper="ForwardEuler") L = 2π gr = TwoDGrid(nx, L, ny, L) - nlayers = 2 # these choice of parameters give the - f0, g = 1, 1 # desired PV-streamfunction relations - H = [0.2, 0.8] # q1 = Δψ1 + 25*(ψ2-ψ1), and - rho = [4.0, 5.0] # q2 = Δψ2 + 25/4*(ψ1-ψ2). - - U1, U2 = 0.1, 0.05 + nlayers = 2 # these choice of parameters give the + f0, 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 Uvectors = zeros(ny, nlayers) Uvectors[:, 1] .= U1 @@ -375,8 +375,8 @@ function test_mqg_paramsconstructor(;dt=0.001, stepper="ForwardEuler") Ufloats[1] = U1 Ufloats[2] = U2 - probUvectors = MultilayerQG.Problem(nlayers=nlayers, nx=nx, ny=ny, Lx=L, f0=f0, g=g, H=H, rho=rho, U=Uvectors) - probUfloats = MultilayerQG.Problem(nlayers=nlayers, nx=nx, ny=ny, Lx=L, f0=f0, g=g, H=H, rho=rho, U=Ufloats) + probUvectors = MultilayerQG.Problem(nlayers=nlayers, nx=nx, ny=ny, Lx=L, f0=f0, g=g, H=H, ρ=ρ, U=Uvectors) + probUfloats = MultilayerQG.Problem(nlayers=nlayers, nx=nx, ny=ny, Lx=L, f0=f0, g=g, H=H, ρ=ρ, U=Ufloats) isapprox(probUfloats.params.U, probUvectors.params.U, rtol=rtol_multilayerqg) end From 00b01422f923071b6ec21602e964d5de51df6e5e Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Wed, 11 Sep 2019 15:36:34 +0300 Subject: [PATCH 2/3] adds test for a 3-layer problem --- test/runtests.jl | 3 +- test/test_multilayerqg.jl | 103 ++++++++++++++++++++++++++++++++------ 2 files changed, 90 insertions(+), 16 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index d988dc03..e08a817b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -90,7 +90,8 @@ end @testset "MultilayerQG" begin include("test_multilayerqg.jl") - @test test_pvtofromstreamfunction() + @test test_pvtofromstreamfunction_2layer() + @test test_pvtofromstreamfunction_3layer() @test test_mqg_nonlinearadvection(0.001, "ForwardEuler") @test test_mqg_linearadvection(0.001, "ForwardEuler") @test test_mqg_energies() diff --git a/test/test_multilayerqg.jl b/test/test_multilayerqg.jl index a2b2c996..d44755bc 100644 --- a/test/test_multilayerqg.jl +++ b/test/test_multilayerqg.jl @@ -1,11 +1,11 @@ """ - constructtestfields() + constructtestfields_2layer(gr) Constructs flow fields for a 2-layer problem with parameters such that q1 = Δψ1 + 25*(ψ2-ψ1), q2 = Δψ2 + 25/4*(ψ1-ψ2). """ -function constructtestfields(gr) +function constructtestfields_2layer(gr) x, y = gridpoints(gr) k0, l0 = gr.k[2], gr.l[2] # fundamental wavenumbers @@ -32,6 +32,36 @@ function constructtestfields(gr) end +""" + constructtestfields_3layer(gr) + +Constructs flow fields for a 3-layer problem with parameters such that + q1 = Δψ1 + 20ψ2 - 20ψ1, + q2 = Δψ2 + 20ψ1 - 44ψ2 + 24ψ3, + q2 = Δψ2 + 12ψ2 - 12ψ3. +""" +function constructtestfields_3layer(gr) + x, y = gridpoints(gr) + k0, l0 = gr.k[2], gr.l[2] # fundamental wavenumbers + + # a set of streafunctions ψ1, ψ2, ψ3, ... + ψ1 = @. 1e-3 * ( 1/4*cos(2k0*x)*cos(5l0*y) + 1/3*cos(3k0*x)*cos(3l0*y) ) + ψ2 = @. 1e-3 * ( cos(3k0*x)*cos(4l0*y) + 1/2*cos(4k0*x)*cos(2l0*y) ) + ψ3 = @. 1e-3 * ( cos(1k0*x)*cos(3l0*y) + 1/2*cos(2k0*x)*cos(2l0*y) ) + + Δψ1 = @. -1e-3 * ( 1/4*((2k0)^2+(5l0)^2)*cos(2k0*x)*cos(5l0*y) + 1/3*((3k0)^2+(3l0)^2)*cos(3k0*x)*cos(3l0*y) ) + Δψ2 = @. -1e-3 * ( ((3k0)^2+(4l0)^2)*cos(3k0*x)*cos(4l0*y) + 1/2*((4k0)^2+(2l0)^2)*cos(4k0*x)*cos(2l0*y) ) + Δψ3 = @. -1e-3 * ( ((1k0)^2+(3l0)^2)*cos(1k0*x)*cos(3l0*y) + 1/2*((2k0)^2+(2l0)^2)*cos(2k0*x)*cos(2l0*y) ) + + # ... their corresponding PVs q1, q2, q3, + q1 = @. Δψ1 + 20ψ2 - 20ψ1 + q2 = @. Δψ2 + 20ψ1 - 44ψ2 + 24ψ3 + q3 = @. Δψ3 + 12ψ2 - 12ψ3 + + return ψ1, ψ2, ψ3, q1, q2, q3 +end + + """ test_pvtofromstreamfunction() @@ -41,19 +71,19 @@ constructs a 2-layer problem with parameters such that Then given ψ1 and ψ2 it tests if pvfromstreamfunction gives the expected q1 and q2. Similarly, that streamfunctionfrompv gives ψ1 and ψ2 from q1 and q2. """ -function test_pvtofromstreamfunction() +function test_pvtofromstreamfunction_2layer() n, L = 128, 2π gr = TwoDGrid(n, L) - nlayers = 2 # these choice of parameters give the - f0, g = 1, 1 # desired PV-streamfunction relations + nlayers = 2 # these choice of parameters give the + f0, 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). prob = MultilayerQG.Problem(nlayers=nlayers, nx=n, Lx=L, f0=f0, g=g, H=H, ρ=ρ) 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(gr) + ψ1, ψ2, q1, q2, ψ1x, ψ2x, q1x, q2x, Δψ2, Δq1, Δq2 = constructtestfields_2layer(gr) vs.psih[:, :, 1] .= rfft(ψ1) vs.psih[:, :, 2] .= rfft(ψ2) @@ -71,6 +101,49 @@ function test_pvtofromstreamfunction() end +""" + test_pvtofromstreamfunction() + +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. +""" +function test_pvtofromstreamfunction_3layer() + n, L = 128, 2π + gr = TwoDGrid(n, L) + + nlayers = 3 # these choice of parameters give the + f0, 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. + + prob = MultilayerQG.Problem(nlayers=nlayers, nx=n, Lx=L, f0=f0, g=g, H=H, ρ=ρ) + sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid + + ψ1, ψ2, ψ3, q1, q2, q3 = constructtestfields_3layer(gr) + + vs.psih[:, :, 1] .= rfft(ψ1) + vs.psih[:, :, 2] .= rfft(ψ2) + vs.psih[:, :, 3] .= rfft(ψ3) + + MultilayerQG.pvfromstreamfunction!(vs.qh, vs.psih, pr.S, gr) + MultilayerQG.invtransform!(vs.q, vs.qh, pr) + + vs.qh[:, :, 1] .= rfft(q1) + vs.qh[:, :, 2] .= rfft(q2) + vs.qh[:, :, 3] .= rfft(q3) + + MultilayerQG.streamfunctionfrompv!(vs.psih, vs.qh, pr.invS, gr) + MultilayerQG.invtransform!(vs.psi, vs.psih, pr) + + isapprox(q1, vs.q[:, :, 1], rtol=rtol_multilayerqg) && isapprox(q2, vs.q[:, :, 2], rtol=rtol_multilayerqg) && isapprox(q3, vs.q[:, :, 3], rtol=rtol_multilayerqg) && isapprox(ψ1, vs.psi[:, :, 1], rtol=rtol_multilayerqg) && isapprox(ψ2, vs.psi[:, :, 2], rtol=rtol_multilayerqg) && isapprox(ψ3, vs.psi[:, :, 3], rtol=rtol_multilayerqg) +end + + """ test_mqg_nonlinearadvection(dt, stepper; kwargs...) @@ -115,7 +188,7 @@ function test_mqg_nonlinearadvection(dt, stepper; n=128, L=2π, nlayers=2, μ=0. η = @. η0*exp( -(x+Lx/8)^2/(2σx^2) -(y-Ly/8)^2/(2σy^2) ) ηx = @. -(x+Lx/8)/(σx^2) * η - ψ1, ψ2, q1, q2, ψ1x, ψ2x, q1x, q2x, Δψ2, Δq1, Δq2 = constructtestfields(gr) + ψ1, ψ2, q1, q2, ψ1x, ψ2x, q1x, q2x, Δψ2, Δq1, Δq2 = constructtestfields_2layer(gr) Ff1 = FourierFlows.jacobian(ψ1, q1, gr) + (β .- uyy1 .- 25*(U2.+u2.-U1.-u1) ).*ψ1x + (U1.+u1).*q1x - ν*Δq1 Ff2 = FourierFlows.jacobian(ψ2, q2 + η, gr) + (β .- uyy2 .- 25/4*(U1.+u1.-U2.-u2) ).*ψ2x + (U2.+u2).*(q2x + ηx) + μ*Δψ2 - ν*Δq2 @@ -198,7 +271,7 @@ function test_mqg_linearadvection(dt, stepper; n=128, L=2π, nlayers=2, μ=0.0, η = @. η0*exp( -(x+Lx/8)^2/(2σx^2) -(y-Ly/8)^2/(2σy^2) ) ηx = @. -(x+Lx/8)/(σx^2) * η - ψ1, ψ2, q1, q2, ψ1x, ψ2x, q1x, q2x, Δψ2, Δq1, Δq2 = constructtestfields(gr) + ψ1, ψ2, q1, q2, ψ1x, ψ2x, q1x, q2x, Δψ2, Δq1, Δq2 = constructtestfields_2layer(gr) 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 @@ -243,7 +316,7 @@ end Tests the kinetic (KE) and potential (PE) energies function by constructing a 2-layer problem and initializing it with a flow field whose KE and PE are known. """ -function test_mqg_energies(;dt=0.001, stepper="ForwardEuler", n=128, L=2π, nlayers=2, μ=0.0, ν=0.0, nν=1) +function test_mqg_energies(; dt=0.001, stepper="ForwardEuler", n=128, L=2π, nlayers=2, μ=0.0, ν=0.0, nν=1) nx, ny = 64, 66 Lx, Ly = 2π, 2π gr = TwoDGrid(nx, Lx, ny, Ly) @@ -259,7 +332,7 @@ function test_mqg_energies(;dt=0.001, stepper="ForwardEuler", n=128, L=2π, nlay prob = MultilayerQG.Problem(nlayers=nlayers, nx=nx, ny=ny, Lx=Lx, Ly=Ly, f0=f0, g=g, H=H, ρ=ρ) 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(gr) + ψ1, ψ2, q1, q2, ψ1x, ψ2x, q1x, q2x, Δψ2, Δq1, Δq2 = constructtestfields_2layer(gr) qf = zeros(gr.nx, gr.ny, nlayers) qf[:, :, 1] .= q1 @@ -278,7 +351,7 @@ end Tests the lateral and vertical eddy fluxes by constructing a 2-layer problem and initializing it with a flow field whose fluxes are known. """ -function test_mqg_fluxes(;dt=0.001, stepper="ForwardEuler", n=128, L=2π, nlayers=2, μ=0.0, ν=0.0, nν=1) +function test_mqg_fluxes(; dt=0.001, stepper="ForwardEuler", n=128, L=2π, nlayers=2, μ=0.0, ν=0.0, nν=1) nx, ny = 128, 126 Lx, Ly = 2π, 2π gr = TwoDGrid(nx, Lx, ny, Ly) @@ -313,7 +386,7 @@ end Tests the set_q!() and set_psi!() functions that initialize sol with a flow with given `q` or `psi` respectively. """ -function test_mqg_setqsetpsi(;dt=0.001, stepper="ForwardEuler", n=64, L=2π, nlayers=2, μ=0.0, ν=0.0, nν=1) +function test_mqg_setqsetpsi(; 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(nx, L, ny, L) @@ -355,13 +428,13 @@ end Tests that `Params` constructor works with both mean flow `U` being a floats (i.e., constant `U` in each layer) or vectors (i.e., `U(y)` in each layer). """ -function test_mqg_paramsconstructor(;dt=0.001, stepper="ForwardEuler") +function test_mqg_paramsconstructor(; dt=0.001, stepper="ForwardEuler") nx, ny = 32, 34 L = 2π gr = TwoDGrid(nx, L, ny, L) - nlayers = 2 # these choice of parameters give the - f0, g = 1, 1 # desired PV-streamfunction relations + nlayers = 2 # these choice of parameters give the + f0, 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). From b683e7a325a4c9d00fffe4858a198a232dce6243 Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Wed, 11 Sep 2019 15:50:11 +0300 Subject: [PATCH 3/3] minor fix in docstrings [ci skip] --- test/test_multilayerqg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_multilayerqg.jl b/test/test_multilayerqg.jl index d44755bc..f14192a8 100644 --- a/test/test_multilayerqg.jl +++ b/test/test_multilayerqg.jl @@ -63,7 +63,7 @@ end """ - test_pvtofromstreamfunction() + test_pvtofromstreamfunction_2layer() Tests the pvfromstreamfunction function that gives qh from ψh. To do so, it constructs a 2-layer problem with parameters such that @@ -102,7 +102,7 @@ end """ - test_pvtofromstreamfunction() + test_pvtofromstreamfunction_3layer() Tests the pvfromstreamfunction function that gives qh from ψh. To do so, it constructs a 3-layer problem with parameters such that