From 7eaad66153c87c57ae32c075f757a49f2e18fc2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pal=C3=B3czy?= Date: Wed, 6 Jul 2022 11:22:07 +0200 Subject: [PATCH 1/7] Fix PE calculation bug in MultiLayerQG.energies --- src/multilayerqg.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/multilayerqg.jl b/src/multilayerqg.jl index 3f42c310..940777a5 100644 --- a/src/multilayerqg.jl +++ b/src/multilayerqg.jl @@ -916,7 +916,7 @@ The kinetic energy at the ``j``-th fluid layer is 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 -𝖯𝖤_{j+1/2} = \\int \\frac1{2} \\frac{f₀^2}{g'_{j+1/2}} (ψ_j - ψ_{j+1})^2 \\frac{𝖽x 𝖽y}{L_x L_y} = \\frac1{2} \\frac{f₀^2}{g'_{j+1/2}} \\sum_{𝐤} |ψ_j - ψ_{j+1}|², \\ j = 1, ..., n-1 . +𝖯𝖤_{j+1/2} = \\int \\frac1{2} \\frac{f₀^2}{g'_{j+1/2} H} (ψ_j - ψ_{j+1})^2 \\frac{𝖽x 𝖽y}{L_x L_y} = \\frac1{2} \\frac{f₀^2}{g'_{j+1/2} H} \\sum_{𝐤} |ψ̂_j - ψ̂_{j+1}|², \\ j = 1, ..., n-1 . ``` """ function energies(vars, params, grid, sol) @@ -934,7 +934,7 @@ function energies(vars, params, grid, sol) end for j = 1:nlayers-1 - CUDA.@allowscalar PE[j] = 1 / (2 * grid.Lx * grid.Ly) * params.f₀^2 / params.g′[j] * parsevalsum(abs2.(vars.ψh[:, :, j+1] .- vars.ψh[:, :, j]), grid) + CUDA.@allowscalar PE[j] = 1 / (2 * grid.Lx * grid.Ly * sum(params.H)) * params.f₀^2 / params.g′[j] * parsevalsum(abs2.(vars.ψh[:, :, j+1] .- vars.ψh[:, :, j]), grid) end return KE, PE @@ -956,7 +956,7 @@ function energies(vars, params::TwoLayerParams, grid, sol) CUDA.@allowscalar KE[j] = @views 1 / (2 * grid.Lx * grid.Ly) * parsevalsum(abs²∇𝐮h[:, :, j], grid) * params.H[j] / sum(params.H) end - PE = @views 1 / (2 * grid.Lx * grid.Ly) * params.f₀^2 / params.g′ * parsevalsum(abs2.(ψ2h .- ψ1h), grid) + PE = @views 1 / (2 * grid.Lx * grid.Ly * sum(params.H)) * params.f₀^2 / params.g′ * parsevalsum(abs2.(ψ2h .- ψ1h), grid) return KE, PE end From 2d3062f4bd35255df2f8acb2ada92b6f004c3bbe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pal=C3=B3czy?= Date: Thu, 7 Jul 2022 09:39:09 +0200 Subject: [PATCH 2/7] Add test for MultiLayerQG.energies when total fluid depth is not unity --- test/runtests.jl | 1 + test/test_multilayerqg.jl | 43 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 44 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 43c7fbbc..9491a4e7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -124,6 +124,7 @@ for dev in devices @test test_mqg_nonlinearadvection(0.005, "ForwardEuler", dev) @test test_mqg_linearadvection(0.005, "ForwardEuler", dev) @test test_mqg_energies(dev) + @test test_mqg_energies_Hneq1(dev) @test test_mqg_energysinglelayer(dev) @test test_mqg_fluxes(dev) @test test_mqg_fluxessinglelayer(dev) diff --git a/test/test_multilayerqg.jl b/test/test_multilayerqg.jl index 849a0876..b49b8190 100644 --- a/test/test_multilayerqg.jl +++ b/test/test_multilayerqg.jl @@ -373,6 +373,49 @@ function test_mqg_energies(dev::Device=CPU(); MultiLayerQG.addforcing!(prob.timestepper.RHS₁, sol, cl.t, cl, vs, pr, gr)==nothing end +""" + test_mqg_energies_Hneq1(dev; kwargs...) + +Tests the kinetic (KE) and potential (PE) energies function (for a case +where the total fluid depth is not unity) by constructing a 2-layer problem +and initializing it with a flow field whose KE and PE are known. +""" +function test_mqg_energies_Hneq1(dev::Device=CPU(); + 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(dev, nx, Lx, ny, Ly) + T = eltype(gr) + + 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, makes no difference for KE1, KE2 but changes PE. + ρ = [1.0, 2.0] # Make g′ = 1/2 + + prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, g, H, ρ) + + ψf = zeros(dev, T, (gr.nx, gr.ny, nlayers)) + ψ1 = @. cos(2k₀*x)*cos(2l₀*y) + ψ2 = ψ1/2 + ψf[:, :, 1] = ψ1 + ψf[:, :, 2] = ψ2 + + MultiLayerQG.set_ψ!(prob, ψf) + + KE1_calc = 1/4 # = H1/H + KE2_calc = 3/4/4 # = H2/H/4 + PE_calc = 1/16/10 # = 1/16/sum(prob.params.H) + + KE, PE = MultiLayerQG.energies(prob) + + return isapprox(KE[1], KE1_calc, rtol=rtol_multilayerqg) && + isapprox(KE[2], KE2_calc, rtol=rtol_multilayerqg) && + isapprox(PE[1], PE_calc, rtol=rtol_multilayerqg) +end + function test_mqg_energysinglelayer(dev::Device=CPU(); dt=0.001, stepper="ForwardEuler", nlayers=1, μ=0.0, ν=0.0, nν=1) nx, Lx = 64, 2π From 3e22654ad2a713a7756704edbcecbcb49099cffe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20Pal=C3=B3czy?= Date: Thu, 7 Jul 2022 10:23:47 +0200 Subject: [PATCH 3/7] Minor clean-ups --- test/test_multilayerqg.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/test/test_multilayerqg.jl b/test/test_multilayerqg.jl index b49b8190..6b77ea96 100644 --- a/test/test_multilayerqg.jl +++ b/test/test_multilayerqg.jl @@ -385,7 +385,6 @@ function test_mqg_energies_Hneq1(dev::Device=CPU(); nx, ny = 64, 66 Lx, Ly = 2π, 2π gr = TwoDGrid(dev, nx, Lx, ny, Ly) - T = eltype(gr) x, y = gridpoints(gr) k₀, l₀ = 2π/gr.Lx, 2π/gr.Ly # fundamental wavenumbers @@ -397,17 +396,16 @@ function test_mqg_energies_Hneq1(dev::Device=CPU(); prob = MultiLayerQG.Problem(nlayers, dev; nx, ny, Lx, Ly, f₀, g, H, ρ) - ψf = zeros(dev, T, (gr.nx, gr.ny, nlayers)) + ψ = zeros(dev, eltype(gr), (gr.nx, gr.ny, nlayers)) ψ1 = @. cos(2k₀*x)*cos(2l₀*y) - ψ2 = ψ1/2 - ψf[:, :, 1] = ψ1 - ψf[:, :, 2] = ψ2 + CUDA.@allowscalar @views ψ[:, :, 1] .= ψ1 + CUDA.@allowscalar @views ψ[:, :, 2] .= ψ1/2 - MultiLayerQG.set_ψ!(prob, ψf) + MultiLayerQG.set_ψ!(prob, ψ) KE1_calc = 1/4 # = H1/H KE2_calc = 3/4/4 # = H2/H/4 - PE_calc = 1/16/10 # = 1/16/sum(prob.params.H) + PE_calc = 1/16/10 # = 1/16/H KE, PE = MultiLayerQG.energies(prob) From 840cf969c066759534c8e63e4d1e9e531f7597a9 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 8 Jul 2022 08:11:16 +1000 Subject: [PATCH 4/7] merge energy tests for multilayerqg --- test/runtests.jl | 1 - test/test_multilayerqg.jl | 58 ++++++--------------------------------- 2 files changed, 9 insertions(+), 50 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 9491a4e7..43c7fbbc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -124,7 +124,6 @@ for dev in devices @test test_mqg_nonlinearadvection(0.005, "ForwardEuler", dev) @test test_mqg_linearadvection(0.005, "ForwardEuler", dev) @test test_mqg_energies(dev) - @test test_mqg_energies_Hneq1(dev) @test test_mqg_energysinglelayer(dev) @test test_mqg_fluxes(dev) @test test_mqg_fluxessinglelayer(dev) diff --git a/test/test_multilayerqg.jl b/test/test_multilayerqg.jl index 6b77ea96..90356c2b 100644 --- a/test/test_multilayerqg.jl +++ b/test/test_multilayerqg.jl @@ -345,73 +345,33 @@ function test_mqg_energies(dev::Device=CPU(); Lx, Ly = 2π, 2π gr = TwoDGrid(dev, nx, Lx, ny, Ly) - 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). - - 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 - - ψ1, ψ2, q1, q2, ψ1x, ψ2x, q1x, q2x, Δψ2, Δq1, Δq2 = constructtestfields_2layer(gr) - - qf = zeros(dev, eltype(gr), (gr.nx, gr.ny, nlayers)) - CUDA.@allowscalar @views qf[:, :, 1] .= q1 - CUDA.@allowscalar @views qf[:, :, 2] .= q2 - - MultiLayerQG.set_q!(prob, qf) - - KE, PE = MultiLayerQG.energies(prob) - - return isapprox(KE[1], 61/640*1e-6, rtol=rtol_multilayerqg) && - isapprox(KE[2], 3*1e-6, rtol=rtol_multilayerqg) && - isapprox(PE[1], 1025/1152*1e-6, rtol=rtol_multilayerqg) && - MultiLayerQG.addforcing!(prob.timestepper.RHS₁, sol, cl.t, cl, vs, pr, gr)==nothing -end - -""" - test_mqg_energies_Hneq1(dev; kwargs...) - -Tests the kinetic (KE) and potential (PE) energies function (for a case -where the total fluid depth is not unity) by constructing a 2-layer problem -and initializing it with a flow field whose KE and PE are known. -""" -function test_mqg_energies_Hneq1(dev::Device=CPU(); - 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(dev, nx, Lx, ny, Ly) - 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, makes no difference for KE1, KE2 but changes PE. + 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, ρ) ψ = zeros(dev, eltype(gr), (gr.nx, gr.ny, nlayers)) - ψ1 = @. cos(2k₀*x)*cos(2l₀*y) - CUDA.@allowscalar @views ψ[:, :, 1] .= ψ1 - CUDA.@allowscalar @views ψ[:, :, 2] .= ψ1/2 + + CUDA.@allowscalar @views ψ[:, :, 1] .= cos(2k₀ * x) * cos(2l₀ * y) + CUDA.@allowscalar @views ψ[:, :, 2] .= 1/2 * cos(2k₀ * x) * cos(2l₀ * y) MultiLayerQG.set_ψ!(prob, ψ) - KE1_calc = 1/4 # = H1/H - KE2_calc = 3/4/4 # = H2/H/4 - PE_calc = 1/16/10 # = 1/16/H + KE1_calc = 1/4 # = H1/H + KE2_calc = 3/4/4 # = H2/H/4 + PE_calc = 1/16/10 # = 1/16/H KE, PE = MultiLayerQG.energies(prob) return isapprox(KE[1], KE1_calc, rtol=rtol_multilayerqg) && isapprox(KE[2], KE2_calc, rtol=rtol_multilayerqg) && - isapprox(PE[1], PE_calc, rtol=rtol_multilayerqg) + isapprox(PE[1], PE_calc, rtol=rtol_multilayerqg) && + MultiLayerQG.addforcing!(prob.timestepper.RHS₁, sol, cl.t, cl, vs, pr, gr)==nothing end function test_mqg_energysinglelayer(dev::Device=CPU(); From 7ef9de8ec436a93a61996d9569638a2359a52181 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 8 Jul 2022 08:12:00 +1000 Subject: [PATCH 5/7] bump patch releast --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6eb55c7b..48fe837a 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.14.0" +version = "0.14.1" [deps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" From d91872a8d98b228d403fd411c6bd9bc38c737780 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 8 Jul 2022 08:16:21 +1000 Subject: [PATCH 6/7] proper streamfunction initialization --- 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 90356c2b..992b76ae 100644 --- a/test/test_multilayerqg.jl +++ b/test/test_multilayerqg.jl @@ -357,8 +357,8 @@ function test_mqg_energies(dev::Device=CPU(); ψ = 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) + CUDA.@allowscalar @views @. ψ[:, :, 1] = cos(2k₀ * x) * cos(2l₀ * y) + CUDA.@allowscalar @views @. ψ[:, :, 2] = 1/2 * cos(2k₀ * x) * cos(2l₀ * y) MultiLayerQG.set_ψ!(prob, ψ) From 43bfa281d7e8a72ce55bbcdf13987dfb18d1b505 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 8 Jul 2022 08:18:08 +1000 Subject: [PATCH 7/7] add shorcuts in test_mqg_energies --- test/test_multilayerqg.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_multilayerqg.jl b/test/test_multilayerqg.jl index 992b76ae..5f973ddd 100644 --- a/test/test_multilayerqg.jl +++ b/test/test_multilayerqg.jl @@ -354,6 +354,7 @@ function test_mqg_energies(dev::Device=CPU(); ρ = [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))