From 16bd801fdbb4695147241eb947a37fcdf010012e Mon Sep 17 00:00:00 2001 From: mncrowe <87329513+mncrowe@users.noreply.github.com> Date: Mon, 17 Jun 2024 10:32:57 +0100 Subject: [PATCH 01/36] Update singlelayerqg.jl to include background flow --- src/singlelayerqg.jl | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 2ce0fc87..61715fe2 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -37,6 +37,7 @@ nothingfunction(args...) = nothing Ly = Lx, β = 0.0, deformation_radius = Inf, + U = 0.0, eta = nothing, ν = 0.0, nν = 1, @@ -62,6 +63,7 @@ Keyword arguments - `Ly`: Extent of the ``y``-domain. - `β`: Planetary vorticity ``y``-gradient. - `deformation_radius`: Rossby radius of deformation; set `Inf` for purely barotropic. + - `U`: Background flow in the ``x``-direction. - `eta`: Topographic potential vorticity. - `ν`: Small-scale (hyper)-viscosity coefficient. - `nν`: (Hyper)-viscosity order, `nν```≥ 1``. @@ -82,6 +84,7 @@ function Problem(dev::Device=CPU(); # Physical parameters β = 0.0, deformation_radius = Inf, + U = 0.0, eta = nothing, # Drag and (hyper-)viscosity ν = 0.0, @@ -103,7 +106,7 @@ function Problem(dev::Device=CPU(); # topographic PV eta === nothing && (eta = zeros(dev, T, (nx, ny))) - params = deformation_radius == Inf ? BarotropicQGParams(grid, T(β), eta, T(μ), T(ν), nν, calcF) : EquivalentBarotropicQGParams(grid, T(β), T(deformation_radius), eta, T(μ), T(ν), nν, calcF) + params = deformation_radius == Inf ? BarotropicQGParams(grid, T(β), T(U), eta, T(μ), T(ν), nν, calcF) : EquivalentBarotropicQGParams(grid, T(β), T(deformation_radius), T(U), eta, T(μ), T(ν), nν, calcF) vars = calcF == nothingfunction ? DecayingVars(grid) : (stochastic ? StochasticForcedVars(grid) : ForcedVars(grid)) @@ -131,6 +134,8 @@ struct Params{T, Aphys, Atrans, ℓ} <: SingleLayerQGParams β :: T "Rossby radius of deformation" deformation_radius :: ℓ + "Background flow in x direction" + U :: T "topographic potential vorticity" eta :: Aphys "Fourier transform of topographic potential vorticity" @@ -149,24 +154,24 @@ const BarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractAr const EquivalentBarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:AbstractFloat} """ - EquivalentBarotropicQGParams(grid, β, deformation_radius, eta, μ, ν, nν, calcF) + EquivalentBarotropicQGParams(grid, β, deformation_radius, U, eta, μ, ν, nν, calcF) Return the parameters for an Equivalent Barotropic QG problem (i.e., with finite Rossby radius of deformation). """ -function EquivalentBarotropicQGParams(grid::AbstractGrid{T, A}, β, deformation_radius, eta, μ, ν, nν::Int, calcF) where {T, A} +function EquivalentBarotropicQGParams(grid::AbstractGrid{T, A}, β, deformation_radius, U, eta, μ, ν, nν::Int, calcF) where {T, A} eta_on_grid = typeof(eta) <: AbstractArray ? A(eta) : FourierFlows.on_grid(eta, grid) etah_on_grid = rfft(eta_on_grid) - return Params(β, deformation_radius, eta_on_grid, etah_on_grid, μ, ν, nν, calcF) + return Params(β, deformation_radius, U, eta_on_grid, etah_on_grid, μ, ν, nν, calcF) end """ - BarotropicQGParams(grid, β, eta, μ, ν, nν, calcF) + BarotropicQGParams(grid, β, U|, eta, μ, ν, nν, calcF) Return the parameters for a Barotropic QG problem (i.e., with infinite Rossby radius of deformation). """ -BarotropicQGParams(grid::AbstractGrid, β, eta, μ, ν, nν::Int, calcF) = - EquivalentBarotropicQGParams(grid, β, nothing, eta, μ, ν, nν, calcF) +BarotropicQGParams(grid::AbstractGrid, β, U, eta, μ, ν, nν::Int, calcF) = + EquivalentBarotropicQGParams(grid, β, nothing, U, eta, μ, ν, nν, calcF) # --------- @@ -177,17 +182,17 @@ BarotropicQGParams(grid::AbstractGrid, β, eta, μ, ν, nν::Int, calcF) = Equation(params::BarotropicQGParams, grid) Return the equation for a barotropic QG problem with `params` and `grid`. Linear operator -``L`` includes bottom drag ``μ``, (hyper)-viscosity of order ``n_ν`` with coefficient ``ν`` -and the ``β`` term: +``L`` includes bottom drag ``μ``, (hyper)-viscosity of order ``n_ν`` with coefficient ``ν``, +the ``β`` term, and a constant background flow ``U``: ```math -L = - μ - ν |𝐤|^{2 n_ν} + i β k_x / |𝐤|² . +L = - μ - ν |𝐤|^{2 n_ν} + i β k_x / |𝐤|² - i U k_x . ``` The nonlinear term is computed via `calcN!` function. """ function Equation(params::BarotropicQGParams, grid::AbstractGrid) - L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr * grid.invKrsq + L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr * grid.invKrsq - im * params.U * grid.kr CUDA.@allowscalar L[1, 1] = 0 return FourierFlows.Equation(L, calcN!, grid) @@ -198,16 +203,16 @@ end Return the equation for an equivalent-barotropic QG problem with `params` and `grid`. Linear operator ``L`` includes bottom drag ``μ``, (hyper)-viscosity of order ``n_ν`` with -coefficient ``ν`` and the ``β`` term: +coefficient ``ν``, the ``β`` term and a constant background flow ``U``: ```math -L = -μ - ν |𝐤|^{2 n_ν} + i β k_x / (|𝐤|² + 1/ℓ²) . +L = -μ - ν |𝐤|^{2 n_ν} + i β k_x / (|𝐤|² + 1/ℓ²) - i U k_x . ``` The nonlinear term is computed via `calcN!` function. """ function Equation(params::EquivalentBarotropicQGParams, grid::AbstractGrid) - L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr / (grid.Krsq + 1 / params.deformation_radius^2) + L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr / (grid.Krsq + 1 / params.deformation_radius^2) - im * params.U * grid.kr CUDA.@allowscalar L[1, 1] = 0 return FourierFlows.Equation(L, calcN!, grid) From 30178e1b34792a4c60423cbb6edf22fc7476e56a Mon Sep 17 00:00:00 2001 From: mncrowe <87329513+mncrowe@users.noreply.github.com> Date: Mon, 17 Jun 2024 10:36:09 +0100 Subject: [PATCH 02/36] Fixed typo in previous change --- src/singlelayerqg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 61715fe2..48814125 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -166,7 +166,7 @@ function EquivalentBarotropicQGParams(grid::AbstractGrid{T, A}, β, deformation_ end """ - BarotropicQGParams(grid, β, U|, eta, μ, ν, nν, calcF) + BarotropicQGParams(grid, β, U, eta, μ, ν, nν, calcF) Return the parameters for a Barotropic QG problem (i.e., with infinite Rossby radius of deformation). """ From 4362a877e4d9be3b542f5992e40b2597332d4706 Mon Sep 17 00:00:00 2001 From: mncrowe <87329513+mncrowe@users.noreply.github.com> Date: Thu, 27 Jun 2024 12:23:31 +0100 Subject: [PATCH 03/36] added background flow test --- test/test_singlelayerqg.jl | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index 2cca28e3..fc6720dd 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -435,4 +435,28 @@ function test_1layerqg_enstrophy_drag(dev; deformation_radius=2.23) prob = SingleLayerQG.Problem(dev; nx=16, deformation_radius=deformation_radius) SingleLayerQG.enstrophy_drag(prob) return nothing -end \ No newline at end of file +end + +""" + test_1layerqg_background_flow(dev; n, L, dt) + +Evolves a lamb vortex with a background flow such that the vortex remains stationary. +""" +function test_1layerqg_background_flow(dev::Device=CPU(); n=256, L=10, dt=0.01) + + prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, dt, U = -1, stepper="FilteredRK4") + x, y = gridpoints(prob.grid) + + q₀ = lambdipole(1,1,prob.grid,center=(1e-10,0)) + SingleLayerQG.set_q!(prob, q₀) + + stepforward!(prob, Int(5/dt)) + + ~,i₀ = findmax(q₀) + ~,i₁ = findmax(prob.vars.q) + + x₀ = prob.grid.x[i₀[1]] # initial vortex position + x₁ = prob.grid.x[i₁[1]] # final vortex position + + return isapprox(x₀, x₁, atol=0.2) +end From 76f1643f268ec0a54e4fd603855db8bc0f713f8f Mon Sep 17 00:00:00 2001 From: mncrowe <87329513+mncrowe@users.noreply.github.com> Date: Thu, 27 Jun 2024 13:32:20 +0100 Subject: [PATCH 04/36] corrected typo --- test/test_singlelayerqg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index fc6720dd..84b8ca8d 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -452,8 +452,8 @@ function test_1layerqg_background_flow(dev::Device=CPU(); n=256, L=10, dt=0.01) stepforward!(prob, Int(5/dt)) - ~,i₀ = findmax(q₀) - ~,i₁ = findmax(prob.vars.q) + _,i₀ = findmax(q₀) + _,i₁ = findmax(prob.vars.q) x₀ = prob.grid.x[i₀[1]] # initial vortex position x₁ = prob.grid.x[i₁[1]] # final vortex position From 4781b8d492df351908689493d2e54ac235d8e8ef Mon Sep 17 00:00:00 2001 From: mncrowe <87329513+mncrowe@users.noreply.github.com> Date: Thu, 27 Jun 2024 23:07:11 +0100 Subject: [PATCH 05/36] Added some spaces to cheer up Greg Co-authored-by: Gregory L. Wagner --- test/test_singlelayerqg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index 84b8ca8d..8f767612 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -447,7 +447,7 @@ function test_1layerqg_background_flow(dev::Device=CPU(); n=256, L=10, dt=0.01) prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, dt, U = -1, stepper="FilteredRK4") x, y = gridpoints(prob.grid) - q₀ = lambdipole(1,1,prob.grid,center=(1e-10,0)) + q₀ = lambdipole(1, 1, prob.grid, center=(1e-10, 0)) SingleLayerQG.set_q!(prob, q₀) stepforward!(prob, Int(5/dt)) From eada586b9e7c58ec1dc7b898c7790ceac729b544 Mon Sep 17 00:00:00 2001 From: mncrowe <87329513+mncrowe@users.noreply.github.com> Date: Thu, 27 Jun 2024 23:18:18 +0100 Subject: [PATCH 06/36] more spaces --- test/test_singlelayerqg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index 8f767612..399f1d6b 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -452,8 +452,8 @@ function test_1layerqg_background_flow(dev::Device=CPU(); n=256, L=10, dt=0.01) stepforward!(prob, Int(5/dt)) - _,i₀ = findmax(q₀) - _,i₁ = findmax(prob.vars.q) + _, i₀ = findmax(q₀) + _, i₁ = findmax(prob.vars.q) x₀ = prob.grid.x[i₀[1]] # initial vortex position x₁ = prob.grid.x[i₁[1]] # final vortex position From cc1b78abfdb46ac43969fc4a3c0aa2fecdadc87e Mon Sep 17 00:00:00 2001 From: Matthew Crowe Date: Mon, 22 Jul 2024 17:32:31 +0100 Subject: [PATCH 07/36] added variable U functionality --- src/singlelayerqg.jl | 62 ++++++++++++++++++++++++++++++++++---------- 1 file changed, 48 insertions(+), 14 deletions(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 48814125..a8e61c17 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -106,7 +106,13 @@ function Problem(dev::Device=CPU(); # topographic PV eta === nothing && (eta = zeros(dev, T, (nx, ny))) - params = deformation_radius == Inf ? BarotropicQGParams(grid, T(β), T(U), eta, T(μ), T(ν), nν, calcF) : EquivalentBarotropicQGParams(grid, T(β), T(deformation_radius), T(U), eta, T(μ), T(ν), nν, calcF) + if U isa Number + U = T(U) + else + U = device_array(dev)(reshape(U,(1, grid.ny))) + end + + params = deformation_radius == Inf ? BarotropicQGParams(grid, T(β), U, eta, T(μ), T(ν), nν, calcF) : EquivalentBarotropicQGParams(grid, T(β), T(deformation_radius), U, eta, T(μ), T(ν), nν, calcF) vars = calcF == nothingfunction ? DecayingVars(grid) : (stochastic ? StochasticForcedVars(grid) : ForcedVars(grid)) @@ -135,7 +141,7 @@ struct Params{T, Aphys, Atrans, ℓ} <: SingleLayerQGParams "Rossby radius of deformation" deformation_radius :: ℓ "Background flow in x direction" - U :: T + U :: Union{T, Aphys} "topographic potential vorticity" eta :: Aphys "Fourier transform of topographic potential vorticity" @@ -172,7 +178,7 @@ Return the parameters for a Barotropic QG problem (i.e., with infinite Rossby ra """ BarotropicQGParams(grid::AbstractGrid, β, U, eta, μ, ν, nν::Int, calcF) = EquivalentBarotropicQGParams(grid, β, nothing, U, eta, μ, ν, nν, calcF) - + # --------- # Equations @@ -192,7 +198,14 @@ L = - μ - ν |𝐤|^{2 n_ν} + i β k_x / |𝐤|² - i U k_x . The nonlinear term is computed via `calcN!` function. """ function Equation(params::BarotropicQGParams, grid::AbstractGrid) - L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr * grid.invKrsq - im * params.U * grid.kr + + if params.U isa Number + U₀ = params.U + else + U₀ = 0.0 + end + + L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr * grid.invKrsq - im * U₀ * grid.kr CUDA.@allowscalar L[1, 1] = 0 return FourierFlows.Equation(L, calcN!, grid) @@ -212,7 +225,14 @@ L = -μ - ν |𝐤|^{2 n_ν} + i β k_x / (|𝐤|² + 1/ℓ²) - i U k_x . The nonlinear term is computed via `calcN!` function. """ function Equation(params::EquivalentBarotropicQGParams, grid::AbstractGrid) - L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr / (grid.Krsq + 1 / params.deformation_radius^2) - im * params.U * grid.kr + + if params.U isa Number + U₀ = params.U + else + U₀ = 0.0 + end + + L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr / (grid.Krsq + 1 / params.deformation_radius^2) - im * U₀ * grid.kr CUDA.@allowscalar L[1, 1] = 0 return FourierFlows.Equation(L, calcN!, grid) @@ -319,6 +339,7 @@ N = - \\widehat{𝖩(ψ, q + η)} = - i k_x \\widehat{u (q + η)} - i k_y \\wide ``` """ function calcN_advection!(N, sol, t, clock, vars, params, grid) + @. vars.qh = sol streamfunctionfrompv!(vars.ψh, vars.qh, params, grid) @. vars.uh = -im * grid.l * vars.ψh @@ -328,17 +349,30 @@ function calcN_advection!(N, sol, t, clock, vars, params, grid) ldiv!(vars.u, grid.rfftplan, vars.uh) ldiv!(vars.v, grid.rfftplan, vars.vh) - uq_plus_η = vars.u # use vars.u as scratch variable - @. uq_plus_η *= vars.q + params.eta # u * (q + η) - vq_plus_η = vars.v # use vars.v as scratch variable - @. vq_plus_η *= vars.q + params.eta # v * (q + η) + if params.U isa Number + + uq_plus_η = vars.u # use vars.u as scratch variable + @. uq_plus_η *= vars.q + params.eta # u * (q + η) + vq_plus_η = vars.v # use vars.v as scratch variable + @. vq_plus_η *= vars.q + params.eta # v * (q + η) + + else + + Q = params.eta .- real.(ifft(im * grid.l .* fft(U))) # PV background (η - ∂U/∂y) + + uq_plus_η = vars.u .+ U # use vars.u as scratch variable + @. uq_plus_η *= vars.q + Q # (u + U) * (q + η - ∂U/∂y) + vq_plus_η = vars.v # use vars.v as scratch variable + @. vq_plus_η *= vars.q + Q # v * (q + η - ∂U/∂y) + + end - uq_plus_ηh = vars.uh # use vars.uh as scratch variable - mul!(uq_plus_ηh, grid.rfftplan, uq_plus_η) # \hat{u * (q + η)} - vq_plus_ηh = vars.vh # use vars.vh as scratch variable - mul!(vq_plus_ηh, grid.rfftplan, vq_plus_η) # \hat{v * (q + η)} + uq_plus_ηh = vars.uh # use vars.uh as scratch variable + mul!(uq_plus_ηh, grid.rfftplan, uq_plus_η) # \hat{(u + U) * (q + η - ∂U/∂y)} + vq_plus_ηh = vars.vh # use vars.vh as scratch variable + mul!(vq_plus_ηh, grid.rfftplan, vq_plus_η) # \hat{v * (q + η - ∂U/∂y)} - @. N = -im * grid.kr * uq_plus_ηh - im * grid.l * vq_plus_ηh # - ∂[u*(q+η)]/∂x - ∂[v*(q+η)]/∂y + @. N = -im * grid.kr * uq_plus_ηh - im * grid.l * vq_plus_ηh # - ∂[(u+U)*(q+η-∂U/∂y)]/∂x - ∂[v*(q+η-∂U/∂y)]/∂y return nothing end From 3ce3abd6bf32b57686c2e53004bc48a5365c9b2f Mon Sep 17 00:00:00 2001 From: Matthew Crowe Date: Mon, 22 Jul 2024 21:51:22 +0100 Subject: [PATCH 08/36] reduced memory of Uy assignment and updated function descriptions --- src/singlelayerqg.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index a8e61c17..61df8cab 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -331,11 +331,11 @@ end """ calcN_advection!(N, sol, t, clock, vars, params, grid) -Calculate the Fourier transform of the advection term, ``- 𝖩(ψ, q+η)`` in conservative -form, i.e., ``- ∂_x[(∂_y ψ)(q+η)] - ∂_y[(∂_x ψ)(q+η)]`` and store it in `N`: +Calculate the Fourier transform of the advection term, ``- 𝖩(ψ+X, q+η-∂U/∂y)`` in conservative +form, i.e., ``- ∂[(u+U)*(q+η-∂U/∂y)]/∂x - ∂[v*(q+η-∂U/∂y)]/∂y`` and store it in `N`: ```math -N = - \\widehat{𝖩(ψ, q + η)} = - i k_x \\widehat{u (q + η)} - i k_y \\widehat{v (q + η)} . +N = - \\widehat{𝖩(ψ + X, q + η - ∂U/∂y)} = - i k_x \\widehat{(u+U) (q + η - ∂U/∂y)} - i k_y \\widehat{v (q + η - ∂U/∂y)} . ``` """ function calcN_advection!(N, sol, t, clock, vars, params, grid) @@ -358,12 +358,12 @@ function calcN_advection!(N, sol, t, clock, vars, params, grid) else - Q = params.eta .- real.(ifft(im * grid.l .* fft(U))) # PV background (η - ∂U/∂y) + Uy = real.(ifft(im * grid.l .* fft(params.U))) # PV background (η - ∂U/∂y) - uq_plus_η = vars.u .+ U # use vars.u as scratch variable - @. uq_plus_η *= vars.q + Q # (u + U) * (q + η - ∂U/∂y) + uq_plus_η = vars.u .+ params.U # use vars.u as scratch variable + @. uq_plus_η *= vars.q + params.eta .- Uy # (u + U) * (q + η - ∂U/∂y) vq_plus_η = vars.v # use vars.v as scratch variable - @. vq_plus_η *= vars.q + Q # v * (q + η - ∂U/∂y) + @. vq_plus_η *= vars.q + params.eta .- Uy # v * (q + η - ∂U/∂y) end From 9bd7576936d7ed7d3f988e8be6fbbcbf3e26538a Mon Sep 17 00:00:00 2001 From: Matthew Crowe Date: Mon, 22 Jul 2024 22:25:09 +0100 Subject: [PATCH 09/36] update function descriptions --- src/singlelayerqg.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 61df8cab..eb6047e1 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -63,7 +63,7 @@ Keyword arguments - `Ly`: Extent of the ``y``-domain. - `β`: Planetary vorticity ``y``-gradient. - `deformation_radius`: Rossby radius of deformation; set `Inf` for purely barotropic. - - `U`: Background flow in the ``x``-direction. + - `U`: Background flow in the ``x``-direction (``y``-dependent). - `eta`: Topographic potential vorticity. - `ν`: Small-scale (hyper)-viscosity coefficient. - `nν`: (Hyper)-viscosity order, `nν```≥ 1``. @@ -189,10 +189,10 @@ BarotropicQGParams(grid::AbstractGrid, β, U, eta, μ, ν, nν::Int, calcF) = Return the equation for a barotropic QG problem with `params` and `grid`. Linear operator ``L`` includes bottom drag ``μ``, (hyper)-viscosity of order ``n_ν`` with coefficient ``ν``, -the ``β`` term, and a constant background flow ``U``: +the ``β`` term, and a constant background flow ``U₀``: ```math -L = - μ - ν |𝐤|^{2 n_ν} + i β k_x / |𝐤|² - i U k_x . +L = - μ - ν |𝐤|^{2 n_ν} + i β k_x / |𝐤|² - i U₀ k_x . ``` The nonlinear term is computed via `calcN!` function. @@ -216,10 +216,10 @@ end Return the equation for an equivalent-barotropic QG problem with `params` and `grid`. Linear operator ``L`` includes bottom drag ``μ``, (hyper)-viscosity of order ``n_ν`` with -coefficient ``ν``, the ``β`` term and a constant background flow ``U``: +coefficient ``ν``, the ``β`` term and a constant background flow ``U₀``: ```math -L = -μ - ν |𝐤|^{2 n_ν} + i β k_x / (|𝐤|² + 1/ℓ²) - i U k_x . +L = -μ - ν |𝐤|^{2 n_ν} + i β k_x / (|𝐤|² + 1/ℓ²) - i U₀ k_x . ``` The nonlinear term is computed via `calcN!` function. @@ -337,6 +337,8 @@ form, i.e., ``- ∂[(u+U)*(q+η-∂U/∂y)]/∂x - ∂[v*(q+η-∂U/∂y)]/∂y` ```math N = - \\widehat{𝖩(ψ + X, q + η - ∂U/∂y)} = - i k_x \\widehat{(u+U) (q + η - ∂U/∂y)} - i k_y \\widehat{v (q + η - ∂U/∂y)} . ``` + +Note: here ∂X/∂y = U. """ function calcN_advection!(N, sol, t, clock, vars, params, grid) @@ -383,7 +385,7 @@ end Calculate the nonlinear term, that is the advection term and the forcing, ```math -N = - \\widehat{𝖩(ψ, q + η)} + F̂ . +N = - \\widehat{𝖩(ψ + X, q + η - ∂U/∂y)} + F̂ . ``` """ function calcN!(N, sol, t, clock, vars, params, grid) From f0baf3f36762128473438c174568fc27d1f99756 Mon Sep 17 00:00:00 2001 From: mncrowe <87329513+mncrowe@users.noreply.github.com> Date: Mon, 22 Jul 2024 22:26:05 +0100 Subject: [PATCH 10/36] T(U) -> convert(T, U) Co-authored-by: Gregory L. Wagner --- src/singlelayerqg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index eb6047e1..f508d95c 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -107,7 +107,7 @@ function Problem(dev::Device=CPU(); eta === nothing && (eta = zeros(dev, T, (nx, ny))) if U isa Number - U = T(U) + U = convert(T, U) else U = device_array(dev)(reshape(U,(1, grid.ny))) end From deea6265b21743cf3a009ec411f4ce0ff4c08e8f Mon Sep 17 00:00:00 2001 From: mncrowe <87329513+mncrowe@users.noreply.github.com> Date: Mon, 22 Jul 2024 22:29:16 +0100 Subject: [PATCH 11/36] clarify U array formulation Co-authored-by: Gregory L. Wagner --- src/singlelayerqg.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index f508d95c..80b3fc27 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -109,7 +109,8 @@ function Problem(dev::Device=CPU(); if U isa Number U = convert(T, U) else - U = device_array(dev)(reshape(U,(1, grid.ny))) + U = reshape(U, (1, grid.ny))) + U = device_array(dev, U) end params = deformation_radius == Inf ? BarotropicQGParams(grid, T(β), U, eta, T(μ), T(ν), nν, calcF) : EquivalentBarotropicQGParams(grid, T(β), T(deformation_radius), U, eta, T(μ), T(ν), nν, calcF) From f351da2d7d1f63988c07f2bb23f67db6e0fc4c2c Mon Sep 17 00:00:00 2001 From: Matthew Crowe Date: Mon, 22 Jul 2024 22:41:08 +0100 Subject: [PATCH 12/36] fix bug with previous push --- src/singlelayerqg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 80b3fc27..586e4a5b 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -109,8 +109,8 @@ function Problem(dev::Device=CPU(); if U isa Number U = convert(T, U) else - U = reshape(U, (1, grid.ny))) - U = device_array(dev, U) + U = reshape(U, (1, grid.ny)) + U = device_array(dev)(U) end params = deformation_radius == Inf ? BarotropicQGParams(grid, T(β), U, eta, T(μ), T(ν), nν, calcF) : EquivalentBarotropicQGParams(grid, T(β), T(deformation_radius), U, eta, T(μ), T(ν), nν, calcF) From 8c1859f42fac03ec9b90c9e94209b3df9447629f Mon Sep 17 00:00:00 2001 From: Matthew Crowe Date: Tue, 23 Jul 2024 14:20:37 +0100 Subject: [PATCH 13/36] add additional test function for U as an Array and include tests in --- test/runtests.jl | 2 ++ test/test_singlelayerqg.jl | 32 ++++++++++++++++++++++++++++++-- 2 files changed, 32 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 6d542b55..d32017f5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -73,6 +73,8 @@ for dev in devices @test test_1layerqg_stochasticforcing_energybudget(dev) @test test_1layerqg_deterministicforcing_enstrophybudget(dev) @test test_1layerqg_stochasticforcing_enstrophybudget(dev) + @test test_1layerqg_background_flow_Num(dev) + @test test_1layerqg_background_flow_Arr(dev) @test SingleLayerQG.nothingfunction() == nothing @test_throws ErrorException("not implemented for finite deformation radius") test_1layerqg_energy_dissipation(dev; deformation_radius=2.23) @test_throws ErrorException("not implemented for finite deformation radius") test_1layerqg_enstrophy_dissipation(dev; deformation_radius=2.23) diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index 399f1d6b..94dac558 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -438,11 +438,13 @@ function test_1layerqg_enstrophy_drag(dev; deformation_radius=2.23) end """ - test_1layerqg_background_flow(dev; n, L, dt) + test_1layerqg_background_flow_Num(dev; n, L, dt) Evolves a lamb vortex with a background flow such that the vortex remains stationary. + +This test uses U as a Number. """ -function test_1layerqg_background_flow(dev::Device=CPU(); n=256, L=10, dt=0.01) +function test_1layerqg_background_flow_Num(dev::Device=CPU(); n=256, L=10, dt=0.01) prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, dt, U = -1, stepper="FilteredRK4") x, y = gridpoints(prob.grid) @@ -460,3 +462,29 @@ function test_1layerqg_background_flow(dev::Device=CPU(); n=256, L=10, dt=0.01) return isapprox(x₀, x₁, atol=0.2) end + +""" + test_1layerqg_background_flow_Arr(dev; n, L, dt) + +Evolves a lamb vortex with a background flow such that the vortex remains stationary. + +This test uses U as an Array. +""" +function test_1layerqg_background_flow_Arr(dev::Device=CPU(); n=256, L=10, dt=0.01) + + prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, dt, U = -ones(n), stepper="FilteredRK4") + x, y = gridpoints(prob.grid) + + q₀ = lambdipole(1, 1, prob.grid, center=(1e-10, 0)) + SingleLayerQG.set_q!(prob, q₀) + + stepforward!(prob, Int(5/dt)) + + _, i₀ = findmax(q₀) + _, i₁ = findmax(prob.vars.q) + + x₀ = prob.grid.x[i₀[1]] # initial vortex position + x₁ = prob.grid.x[i₁[1]] # final vortex position + + return isapprox(x₀, x₁, atol=0.2) +end From dd93d2234bcda29d2ed5230b68ec43e596a22c1b Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 27 Jul 2024 14:59:03 +1000 Subject: [PATCH 14/36] simplify background U and use Rossby wave test for U that does not vary in y --- src/singlelayerqg.jl | 63 ++++++++++++++++++-------------------- test/runtests.jl | 24 +++++++-------- test/test_singlelayerqg.jl | 63 +++----------------------------------- 3 files changed, 46 insertions(+), 104 deletions(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 586e4a5b..8d7b2452 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -63,7 +63,7 @@ Keyword arguments - `Ly`: Extent of the ``y``-domain. - `β`: Planetary vorticity ``y``-gradient. - `deformation_radius`: Rossby radius of deformation; set `Inf` for purely barotropic. - - `U`: Background flow in the ``x``-direction (``y``-dependent). + - `U`: Imposed background constant zonal flow ``U(y)``. - `eta`: Topographic potential vorticity. - `ν`: Small-scale (hyper)-viscosity coefficient. - `nν`: (Hyper)-viscosity order, `nν```≥ 1``. @@ -107,10 +107,10 @@ function Problem(dev::Device=CPU(); eta === nothing && (eta = zeros(dev, T, (nx, ny))) if U isa Number - U = convert(T, U) + U = convert(T, U) else - U = reshape(U, (1, grid.ny)) - U = device_array(dev)(U) + U = reshape(U, (1, grid.ny)) + U = device_array(dev)(U) end params = deformation_radius == Inf ? BarotropicQGParams(grid, T(β), U, eta, T(μ), T(ν), nν, calcF) : EquivalentBarotropicQGParams(grid, T(β), T(deformation_radius), U, eta, T(μ), T(ν), nν, calcF) @@ -141,7 +141,7 @@ struct Params{T, Aphys, Atrans, ℓ} <: SingleLayerQGParams β :: T "Rossby radius of deformation" deformation_radius :: ℓ - "Background flow in x direction" + "Background flow in ``x`` direction" U :: Union{T, Aphys} "topographic potential vorticity" eta :: Aphys @@ -190,25 +190,24 @@ BarotropicQGParams(grid::AbstractGrid, β, U, eta, μ, ν, nν::Int, calcF) = Return the equation for a barotropic QG problem with `params` and `grid`. Linear operator ``L`` includes bottom drag ``μ``, (hyper)-viscosity of order ``n_ν`` with coefficient ``ν``, -the ``β`` term, and a constant background flow ``U₀``: +the ``β`` term, and a constant background flow ``U``: ```math -L = - μ - ν |𝐤|^{2 n_ν} + i β k_x / |𝐤|² - i U₀ k_x . +L = - μ - ν |𝐤|^{2 n_ν} + i β k_x / |𝐤|² - i k_x U . ``` The nonlinear term is computed via `calcN!` function. """ function Equation(params::BarotropicQGParams, grid::AbstractGrid) + L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr * grid.invKrsq + if params.U isa Number - U₀ = params.U - else - U₀ = 0.0 + @. L -= im * params.U * grid.kr end - L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr * grid.invKrsq - im * U₀ * grid.kr CUDA.@allowscalar L[1, 1] = 0 - + return FourierFlows.Equation(L, calcN!, grid) end @@ -217,25 +216,23 @@ end Return the equation for an equivalent-barotropic QG problem with `params` and `grid`. Linear operator ``L`` includes bottom drag ``μ``, (hyper)-viscosity of order ``n_ν`` with -coefficient ``ν``, the ``β`` term and a constant background flow ``U₀``: +coefficient ``ν``, the ``β`` term and a constant background flow ``U``: ```math -L = -μ - ν |𝐤|^{2 n_ν} + i β k_x / (|𝐤|² + 1/ℓ²) - i U₀ k_x . +L = -μ - ν |𝐤|^{2 n_ν} + i β k_x / (|𝐤|² + 1/ℓ²) - i k_x U . ``` The nonlinear term is computed via `calcN!` function. """ function Equation(params::EquivalentBarotropicQGParams, grid::AbstractGrid) + L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr / (grid.Krsq + 1 / params.deformation_radius^2) if params.U isa Number - U₀ = params.U - else - U₀ = 0.0 + @. L -= im * params.U * grid.kr end - L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr / (grid.Krsq + 1 / params.deformation_radius^2) - im * U₀ * grid.kr CUDA.@allowscalar L[1, 1] = 0 - + return FourierFlows.Equation(L, calcN!, grid) end @@ -354,28 +351,28 @@ function calcN_advection!(N, sol, t, clock, vars, params, grid) if params.U isa Number - uq_plus_η = vars.u # use vars.u as scratch variable - @. uq_plus_η *= vars.q + params.eta # u * (q + η) - vq_plus_η = vars.v # use vars.v as scratch variable - @. vq_plus_η *= vars.q + params.eta # v * (q + η) + uq_plus_η = vars.u # use vars.u as scratch variable + @. uq_plus_η *= vars.q + params.eta # u * (q + η) + vq_plus_η = vars.v # use vars.v as scratch variable + @. vq_plus_η *= vars.q + params.eta # v * (q + η) else - Uy = real.(ifft(im * grid.l .* fft(params.U))) # PV background (η - ∂U/∂y) + Uy = real.(ifft(im * grid.l .* fft(params.U))) # PV background (η - ∂U/∂y) - uq_plus_η = vars.u .+ params.U # use vars.u as scratch variable - @. uq_plus_η *= vars.q + params.eta .- Uy # (u + U) * (q + η - ∂U/∂y) - vq_plus_η = vars.v # use vars.v as scratch variable - @. vq_plus_η *= vars.q + params.eta .- Uy # v * (q + η - ∂U/∂y) + uq_plus_η = vars.u .+ params.U # use vars.u as scratch variable + @. uq_plus_η *= vars.q + params.eta .- Uy # (u + U) * (q + η - ∂U/∂y) + vq_plus_η = vars.v # use vars.v as scratch variable + @. vq_plus_η *= vars.q + params.eta .- Uy # v * (q + η - ∂U/∂y) end - uq_plus_ηh = vars.uh # use vars.uh as scratch variable - mul!(uq_plus_ηh, grid.rfftplan, uq_plus_η) # \hat{(u + U) * (q + η - ∂U/∂y)} - vq_plus_ηh = vars.vh # use vars.vh as scratch variable - mul!(vq_plus_ηh, grid.rfftplan, vq_plus_η) # \hat{v * (q + η - ∂U/∂y)} + uq_plus_ηh = vars.uh # use vars.uh as scratch variable + mul!(uq_plus_ηh, grid.rfftplan, uq_plus_η) # \hat{(u + U) * (q + η - ∂U/∂y)} + vq_plus_ηh = vars.vh # use vars.vh as scratch variable + mul!(vq_plus_ηh, grid.rfftplan, vq_plus_η) # \hat{v * (q + η - ∂U/∂y)} - @. N = -im * grid.kr * uq_plus_ηh - im * grid.l * vq_plus_ηh # - ∂[(u+U)*(q+η-∂U/∂y)]/∂x - ∂[v*(q+η-∂U/∂y)]/∂y + @. N = -im * grid.kr * uq_plus_ηh - im * grid.l * vq_plus_ηh # - ∂[(u+U)*(q+η-∂U/∂y)]/∂x - ∂[v*(q+η-∂U/∂y)]/∂y return nothing end diff --git a/test/runtests.jl b/test/runtests.jl index d32017f5..e5dd9100 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -54,16 +54,16 @@ for dev in devices @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) - @test test_1layerqg_rossbywave("RK4", 1e-2, 20, dev, deformation_radius=deformation_radius) - @test test_1layerqg_rossbywave("FilteredRK4", 1e-2, 20, dev, deformation_radius=deformation_radius) - @test test_1layerqg_rossbywave("AB3", 1e-3, 200, dev, deformation_radius=deformation_radius) - @test test_1layerqg_rossbywave("FilteredAB3", 1e-3, 200, dev, deformation_radius=deformation_radius) - @test test_1layerqg_rossbywave("ForwardEuler", 1e-4, 2000, dev, deformation_radius=deformation_radius) - @test test_1layerqg_rossbywave("FilteredForwardEuler", 1e-4, 2000, dev, deformation_radius=deformation_radius) - @test test_1layerqg_problemtype(dev, Float32, deformation_radius=deformation_radius) + for deformation_radius in (Inf, 1.23), U₀ in (0, 0.3) + for (timestepper, dt, nsteps) in zip(("ETDRK4", "FilteredETDRK4", "RK4", "FilteredRK4", "AB3", "FilteredAB3", "ForwardEuler", "FilteredForwardEuler",), + (1e-2, 1e-2, 1e-2, 1e-2, 1e-3, 1e-3, 1e-4, 1e-4,), + (20, 20, 20, 20, 200, 200, 2000, 2000,)) + + nx = 64 + @test test_1layerqg_rossbywave(timestepper, dt, nsteps, dev, nx; deformation_radius, U₀) + @test test_1layerqg_rossbywave(timestepper, dt, nsteps, dev, nx; deformation_radius, U₀=U₀*zeros(nx)) + @test test_1layerqg_problemtype(dev, Float32; deformation_radius) + end end @test test_1layerqg_advection(0.0005, "ForwardEuler", dev) @test test_streamfunctionfrompv(dev; deformation_radius=1.23) @@ -73,9 +73,7 @@ for dev in devices @test test_1layerqg_stochasticforcing_energybudget(dev) @test test_1layerqg_deterministicforcing_enstrophybudget(dev) @test test_1layerqg_stochasticforcing_enstrophybudget(dev) - @test test_1layerqg_background_flow_Num(dev) - @test test_1layerqg_background_flow_Arr(dev) - @test SingleLayerQG.nothingfunction() == nothing + @test SingleLayerQG.nothingfunction() === nothing @test_throws ErrorException("not implemented for finite deformation radius") test_1layerqg_energy_dissipation(dev; deformation_radius=2.23) @test_throws ErrorException("not implemented for finite deformation radius") test_1layerqg_enstrophy_dissipation(dev; deformation_radius=2.23) @test_throws ErrorException("not implemented for finite deformation radius") test_1layerqg_energy_work(dev; deformation_radius=2.23) diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index 94dac558..4ef30fa7 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -3,8 +3,7 @@ Evolves a Rossby wave and compares with the analytic solution. """ -function test_1layerqg_rossbywave(stepper, dt, nsteps, dev::Device=CPU(); deformation_radius=Inf) - nx = 64 +function test_1layerqg_rossbywave(stepper, dt, nsteps, dev::Device=CPU(), nx=64; deformation_radius=Inf, U₀ = 0) Lx = 2π β = 2.0 μ = 0.0 @@ -19,15 +18,15 @@ function test_1layerqg_rossbywave(stepper, dt, nsteps, dev::Device=CPU(); deform eta(x, y) = 0 * x end - prob = SingleLayerQG.Problem(dev; nx=nx, Lx=Lx, eta=eta, deformation_radius=deformation_radius, β=β, μ=μ, ν=ν, stepper=stepper, dt=dt) + prob = SingleLayerQG.Problem(dev; nx=nx, Lx=Lx, eta=eta, deformation_radius=deformation_radius, β=β, U=U₀, μ=μ, ν=ν, stepper=stepper, dt=dt) sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid x, y = gridpoints(grid) # the Rossby wave initial condition ampl = 1e-2 - kwave = 3 * 2π/grid.Lx - lwave = 2 * 2π/grid.Ly + kwave = 3 * 2π / grid.Lx + lwave = 2 * 2π / grid.Ly ω = -params.β * kwave / (kwave^2 + lwave^2 + 1 / deformation_radius^2) q₀ = @. ampl * cos(kwave * x) * cos(lwave * y) q₀h = rfft(q₀) @@ -38,7 +37,7 @@ function test_1layerqg_rossbywave(stepper, dt, nsteps, dev::Device=CPU(); deform dealias!(sol, grid) SingleLayerQG.updatevars!(prob) - q_theory = @. ampl * cos(kwave * (x - ω / kwave * clock.t)) * cos(lwave * y) + q_theory = @. ampl * cos(kwave * (x - (ω / kwave + U₀) * clock.t)) * cos(lwave * y) return isapprox(q_theory, vars.q, rtol=grid.nx * grid.ny * nsteps * 1e-12) end @@ -436,55 +435,3 @@ function test_1layerqg_enstrophy_drag(dev; deformation_radius=2.23) SingleLayerQG.enstrophy_drag(prob) return nothing end - -""" - test_1layerqg_background_flow_Num(dev; n, L, dt) - -Evolves a lamb vortex with a background flow such that the vortex remains stationary. - -This test uses U as a Number. -""" -function test_1layerqg_background_flow_Num(dev::Device=CPU(); n=256, L=10, dt=0.01) - - prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, dt, U = -1, stepper="FilteredRK4") - x, y = gridpoints(prob.grid) - - q₀ = lambdipole(1, 1, prob.grid, center=(1e-10, 0)) - SingleLayerQG.set_q!(prob, q₀) - - stepforward!(prob, Int(5/dt)) - - _, i₀ = findmax(q₀) - _, i₁ = findmax(prob.vars.q) - - x₀ = prob.grid.x[i₀[1]] # initial vortex position - x₁ = prob.grid.x[i₁[1]] # final vortex position - - return isapprox(x₀, x₁, atol=0.2) -end - -""" - test_1layerqg_background_flow_Arr(dev; n, L, dt) - -Evolves a lamb vortex with a background flow such that the vortex remains stationary. - -This test uses U as an Array. -""" -function test_1layerqg_background_flow_Arr(dev::Device=CPU(); n=256, L=10, dt=0.01) - - prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, dt, U = -ones(n), stepper="FilteredRK4") - x, y = gridpoints(prob.grid) - - q₀ = lambdipole(1, 1, prob.grid, center=(1e-10, 0)) - SingleLayerQG.set_q!(prob, q₀) - - stepforward!(prob, Int(5/dt)) - - _, i₀ = findmax(q₀) - _, i₁ = findmax(prob.vars.q) - - x₀ = prob.grid.x[i₀[1]] # initial vortex position - x₁ = prob.grid.x[i₁[1]] # final vortex position - - return isapprox(x₀, x₁, atol=0.2) -end From 122b713631058882baaebe1b256d17d5fd0f3af7 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 27 Jul 2024 14:59:15 +1000 Subject: [PATCH 15/36] homogenize notation/docs --- src/multilayerqg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/multilayerqg.jl b/src/multilayerqg.jl index 4e35d679..05c41b0e 100644 --- a/src/multilayerqg.jl +++ b/src/multilayerqg.jl @@ -66,7 +66,7 @@ Keyword arguments - `Ly`: Extent of the ``y``-domain. - `f₀`: Constant planetary vorticity. - `β`: Planetary vorticity ``y``-gradient. - - `U`: The imposed constant zonal flow ``U(y)`` in each fluid layer. + - `U`: Imposed background constant zonal flow ``U(y)`` in each fluid layer. - `H`: Rest height of each fluid layer. - `b`: Boussinesq buoyancy of each fluid layer. - `eta`: Periodic component of the topographic potential vorticity. From ab0fe027f7a16febb24d7e23028a8cb1fa037436 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 27 Jul 2024 14:59:23 +1000 Subject: [PATCH 16/36] bump patch release --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a49c334c..452978c4 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 contributors"] -version = "0.16.1" +version = "0.16.2" [deps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" From 85d4393176f30f1fcbfa3fa778364a27e3cb7cd3 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 27 Jul 2024 15:17:04 +1000 Subject: [PATCH 17/36] fix test for gpu --- test/runtests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index e5dd9100..1ceffbfc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -48,7 +48,7 @@ for dev in devices @test test_twodnavierstokes_stochasticforcing_enstrophybudget(dev) @test test_twodnavierstokes_energyenstrophypalinstrophy(dev) @test test_twodnavierstokes_problemtype(dev, Float32) - @test TwoDNavierStokes.nothingfunction() == nothing + @test TwoDNavierStokes.nothingfunction() === nothing end @testset "SingleLayerQG" begin @@ -61,7 +61,7 @@ for dev in devices nx = 64 @test test_1layerqg_rossbywave(timestepper, dt, nsteps, dev, nx; deformation_radius, U₀) - @test test_1layerqg_rossbywave(timestepper, dt, nsteps, dev, nx; deformation_radius, U₀=U₀*zeros(nx)) + @test test_1layerqg_rossbywave(timestepper, dt, nsteps, dev, nx; deformation_radius, U₀=U₀*zeros(dev, Float64, (nx,))) @test test_1layerqg_problemtype(dev, Float32; deformation_radius) end end From 619f3988e96a925c25605b5f52f0f3646b219249 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 27 Jul 2024 15:20:02 +1000 Subject: [PATCH 18/36] == -> === --- test/runtests.jl | 4 ++-- test/test_singlelayerqg.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 1ceffbfc..70059ba4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -98,7 +98,7 @@ for dev in devices @test test_bqgql_advection(0.0005, "ForwardEuler", dev) @test test_bqgql_energyenstrophy(dev) @test test_bqgql_problemtype(dev, Float32) - @test BarotropicQGQL.nothingfunction() == nothing + @test BarotropicQGQL.nothingfunction() === nothing end @testset "SurfaceQG" begin @@ -112,7 +112,7 @@ for dev in devices @test test_sqg_problemtype(dev, Float32) @test test_sqg_paramsconstructor(dev) @test test_sqg_noforcing(dev) - @test SurfaceQG.nothingfunction() == nothing + @test SurfaceQG.nothingfunction() === nothing end @testset "MultiLayerQG" begin diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index 4ef30fa7..f9554b17 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -328,7 +328,7 @@ function test_1layerqg_energyenstrophy_BarotropicQG(dev::Device=CPU()) enstrophyq₀ = SingleLayerQG.enstrophy(prob) return isapprox(energyq₀, energy_calc, rtol=rtol_singlelayerqg) && isapprox(enstrophyq₀, enstrophy_calc, rtol=rtol_singlelayerqg) && SingleLayerQG.potential_energy(prob)==0 && - SingleLayerQG.addforcing!(prob.timestepper.N, sol, clock.t, clock, vars, params, grid) == nothing + SingleLayerQG.addforcing!(prob.timestepper.N, sol, clock.t, clock, vars, params, grid) === nothing end """ @@ -361,7 +361,7 @@ function test_1layerqg_energies_EquivalentBarotropicQG(dev; deformation_radius=1 energyq₀ = SingleLayerQG.energy(prob) return isapprox(kinetic_energyq₀, kinetic_energy_calc, rtol=rtol_singlelayerqg) && isapprox(potential_energyq₀, potential_energy_calc, rtol=rtol_singlelayerqg) && isapprox(energyq₀, energy_calc, rtol=rtol_singlelayerqg) && - SingleLayerQG.addforcing!(prob.timestepper.N, sol, clock.t, clock, vars, params, grid) == nothing + SingleLayerQG.addforcing!(prob.timestepper.N, sol, clock.t, clock, vars, params, grid) === nothing end """ From ace06ebbd233c7f331a5fe5618e6d96f42a914dc Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 27 Jul 2024 16:10:53 +1000 Subject: [PATCH 19/36] use different calcN_advection! methods for U const and y-varying --- src/singlelayerqg.jl | 66 ++++++++++++++++++++++++++++---------------- 1 file changed, 42 insertions(+), 24 deletions(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 8d7b2452..3868def3 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -136,13 +136,13 @@ The parameters for the `SingleLayerQG` problem. $(TYPEDFIELDS) """ -struct Params{T, Aphys, Atrans, ℓ} <: SingleLayerQGParams +struct Params{T, Aphys, Atrans, ℓ, TU <: Union{T, Aphys}} <: SingleLayerQGParams "planetary vorticity ``y``-gradient" β :: T "Rossby radius of deformation" deformation_radius :: ℓ "Background flow in ``x`` direction" - U :: Union{T, Aphys} + U :: TU "topographic potential vorticity" eta :: Aphys "Fourier transform of topographic potential vorticity" @@ -160,6 +160,9 @@ end const BarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, Nothing} const EquivalentBarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:AbstractFloat} +const SingleLayerQGconstantUParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:Any, <:AbstractFloat} +const SingleLayerQGvaryingUParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:Any, <:AbstractArray} + """ EquivalentBarotropicQGParams(grid, β, deformation_radius, U, eta, μ, ν, nν, calcF) @@ -329,16 +332,16 @@ end """ calcN_advection!(N, sol, t, clock, vars, params, grid) -Calculate the Fourier transform of the advection term, ``- 𝖩(ψ+X, q+η-∂U/∂y)`` in conservative -form, i.e., ``- ∂[(u+U)*(q+η-∂U/∂y)]/∂x - ∂[v*(q+η-∂U/∂y)]/∂y`` and store it in `N`: +Calculate the Fourier transform of the advection term, ``- 𝖩(ψ, q+η) - U ∂(q+η)/∂x - v ∂U/∂y`` in conservative +form, i.e., ``- ∂[(u + U)(q + η - ∂U/∂y)]/∂x - ∂[v(q + η - ∂U/∂y)]/∂y`` and store it in `N`: ```math -N = - \\widehat{𝖩(ψ + X, q + η - ∂U/∂y)} = - i k_x \\widehat{(u+U) (q + η - ∂U/∂y)} - i k_y \\widehat{v (q + η - ∂U/∂y)} . +N = - \\widehat{𝖩(ψ + Ψ, q + η - ∂U/∂y)} = - i k_x \\widehat{(u+U) (q + η - ∂U/∂y)} - i k_y \\widehat{v (q + η - ∂U/∂y)} . ``` -Note: here ∂X/∂y = U. +Note: here ``- ∂Ψ/∂y = U``. """ -function calcN_advection!(N, sol, t, clock, vars, params, grid) +function calcN_advection!(N, sol, t, clock, vars, params::SingleLayerQGconstantUParams, grid) @. vars.qh = sol streamfunctionfrompv!(vars.ψh, vars.qh, params, grid) @@ -349,30 +352,45 @@ function calcN_advection!(N, sol, t, clock, vars, params, grid) ldiv!(vars.u, grid.rfftplan, vars.uh) ldiv!(vars.v, grid.rfftplan, vars.vh) - if params.U isa Number + uq_plus_η = vars.u # use vars.u as scratch variable + @. uq_plus_η *= vars.q + params.eta # u * (q + η) + vq_plus_η = vars.v # use vars.v as scratch variable + @. vq_plus_η *= vars.q + params.eta # v * (q + η) - uq_plus_η = vars.u # use vars.u as scratch variable - @. uq_plus_η *= vars.q + params.eta # u * (q + η) - vq_plus_η = vars.v # use vars.v as scratch variable - @. vq_plus_η *= vars.q + params.eta # v * (q + η) + uq_plus_ηh = vars.uh # use vars.uh as scratch variable + mul!(uq_plus_ηh, grid.rfftplan, uq_plus_η) # \hat{u * (q + η)} + vq_plus_ηh = vars.vh # use vars.vh as scratch variable + mul!(vq_plus_ηh, grid.rfftplan, vq_plus_η) # \hat{v * (q + η)} - else + @. N = -im * grid.kr * uq_plus_ηh - im * grid.l * vq_plus_ηh # - ∂[u*(q+η)]/∂x - ∂[v*(q+η)]/∂y - Uy = real.(ifft(im * grid.l .* fft(params.U))) # PV background (η - ∂U/∂y) + return nothing +end - uq_plus_η = vars.u .+ params.U # use vars.u as scratch variable - @. uq_plus_η *= vars.q + params.eta .- Uy # (u + U) * (q + η - ∂U/∂y) - vq_plus_η = vars.v # use vars.v as scratch variable - @. vq_plus_η *= vars.q + params.eta .- Uy # v * (q + η - ∂U/∂y) +function calcN_advection!(N, sol, t, clock, vars, params::SingleLayerQGvaryingUParams, grid) - end + @. vars.qh = sol + streamfunctionfrompv!(vars.ψh, vars.qh, params, grid) + @. vars.uh = -im * grid.l * vars.ψh + @. vars.vh = im * grid.kr * vars.ψh + + ldiv!(vars.q, grid.rfftplan, vars.qh) + ldiv!(vars.u, grid.rfftplan, vars.uh) + ldiv!(vars.v, grid.rfftplan, vars.vh) + + Uy = real.(ifft(im * grid.l .* fft(params.U))) # PV background (η - ∂U/∂y) + + uq_plus_η = vars.u .+ params.U # use vars.u as scratch variable + @. uq_plus_η *= vars.q + params.eta .- Uy # (u + U) * (q + η - ∂U/∂y) + vq_plus_η = vars.v # use vars.v as scratch variable + @. vq_plus_η *= vars.q + params.eta .- Uy # v * (q + η - ∂U/∂y) - uq_plus_ηh = vars.uh # use vars.uh as scratch variable - mul!(uq_plus_ηh, grid.rfftplan, uq_plus_η) # \hat{(u + U) * (q + η - ∂U/∂y)} - vq_plus_ηh = vars.vh # use vars.vh as scratch variable - mul!(vq_plus_ηh, grid.rfftplan, vq_plus_η) # \hat{v * (q + η - ∂U/∂y)} + uq_plus_ηh = vars.uh # use vars.uh as scratch variable + mul!(uq_plus_ηh, grid.rfftplan, uq_plus_η) # \hat{(u + U) * (q + η - ∂U/∂y)} + vq_plus_ηh = vars.vh # use vars.vh as scratch variable + mul!(vq_plus_ηh, grid.rfftplan, vq_plus_η) # \hat{v * (q + η - ∂U/∂y)} - @. N = -im * grid.kr * uq_plus_ηh - im * grid.l * vq_plus_ηh # - ∂[(u+U)*(q+η-∂U/∂y)]/∂x - ∂[v*(q+η-∂U/∂y)]/∂y + @. N = -im * grid.kr * uq_plus_ηh - im * grid.l * vq_plus_ηh # - ∂[(u+U)*(q+η-∂U/∂y)]/∂x - ∂[v*(q+η-∂U/∂y)]/∂y return nothing end From 2fcccc772bdca326ad49e9f99f570df7bbd7f51b Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 27 Jul 2024 17:26:33 +1000 Subject: [PATCH 20/36] minor tweaks --- src/singlelayerqg.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 3868def3..24d6ac69 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -130,17 +130,17 @@ end abstract type SingleLayerQGParams <: AbstractParams end """ - struct Params{T, Aphys, Atrans, ℓ} <: SingleLayerQGParams + struct Params{T, Aphys, Atrans, Tℓ, TU} <: SingleLayerQGParams The parameters for the `SingleLayerQG` problem. $(TYPEDFIELDS) """ -struct Params{T, Aphys, Atrans, ℓ, TU <: Union{T, Aphys}} <: SingleLayerQGParams +struct Params{T, Aphys, Atrans, Tℓ, TU <: Union{T, Aphys}} <: SingleLayerQGParams "planetary vorticity ``y``-gradient" β :: T "Rossby radius of deformation" - deformation_radius :: ℓ + deformation_radius :: Tℓ "Background flow in ``x`` direction" U :: TU "topographic potential vorticity" @@ -157,7 +157,7 @@ struct Params{T, Aphys, Atrans, ℓ, TU <: Union{T, Aphys}} <: SingleLayerQGPara calcF! :: Function end -const BarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, Nothing} +const BarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, Nothing} const EquivalentBarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:AbstractFloat} const SingleLayerQGconstantUParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:Any, <:AbstractFloat} @@ -202,7 +202,6 @@ L = - μ - ν |𝐤|^{2 n_ν} + i β k_x / |𝐤|² - i k_x U . The nonlinear term is computed via `calcN!` function. """ function Equation(params::BarotropicQGParams, grid::AbstractGrid) - L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr * grid.invKrsq if params.U isa Number From 20a361bd4843f3c5bde06d79a8e355a44ec88e58 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 09:20:18 +1000 Subject: [PATCH 21/36] update singlelayergq module in docs --- docs/src/modules/singlelayerqg.md | 30 ++++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/docs/src/modules/singlelayerqg.md b/docs/src/modules/singlelayerqg.md index 6f6e382c..1e63d9d5 100644 --- a/docs/src/modules/singlelayerqg.md +++ b/docs/src/modules/singlelayerqg.md @@ -20,15 +20,24 @@ radius follows is said to obey equivalent-barotropic dynamics. We denote the sum vorticity and the vortex stretching contributions to the QGPV with ``q \equiv \nabla^2 \psi - \psi / \ell^2``. Also, we denote the topographic PV with ``\eta \equiv f_0 h / H``. -The dynamical variable is ``q``. Thus, the equation solved by the module is: +The dynamical variable is ``q``. Including an imposed zonal flow ``U(y)``, the equation of motion is: ```math -\partial_t q + \mathsf{J}(\psi, q + \eta) + \beta \partial_x \psi = -\underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu} \right] q}_{\textrm{dissipation}} + F . +\partial_t q + \mathsf{J}(\psi, q) + (U - \partial_y\psi) \partial_x Q + U \partial_x q + (\partial_y Q)(\partial_x \psi) = \underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu} \right] q}_{\textrm{dissipation}} + F , ``` -where ``\mathsf{J}(a, b) = (\partial_x a)(\partial_y b)-(\partial_y a)(\partial_x b)`` is the -two-dimensional Jacobian. On the right hand side, ``F(x, y, t)`` is forcing, ``\mu`` is +with + +```math +\begin{aligned} +\partial_y Q &\equiv \beta - \partial_y^2 U + \partial_y \eta , \\ +\partial_x Q &\equiv \partial_x \eta , +\end{aligned} +``` + +the background PV gradient components, and with +``\mathsf{J}(a, b) = (\partial_x a)(\partial_y b) - (\partial_y a)(\partial_x b)`` +the two-dimensional Jacobian. On the right hand side, ``F(x, y, t)`` is forcing, ``\mu`` is linear drag, and ``\nu`` is hyperviscosity of order ``n_\nu``. Plain old viscosity corresponds to ``n_\nu = 1``. @@ -38,15 +47,16 @@ to ``n_\nu = 1``. The equation is time-stepped forward in Fourier space: ```math -\partial_t \widehat{q} = - \widehat{\mathsf{J}(\psi, q + \eta)} + \beta \frac{i k_x}{|𝐤|^2 + 1/\ell^2} \widehat{q} - \left(\mu + \nu |𝐤|^{2n_\nu} \right) \widehat{q} + \widehat{F} . +\partial_t \widehat{q} = - \widehat{\mathsf{J}(\psi, q)} - \widehat{U \partial_x Q} - \widehat{U \partial_x q} ++ \widehat{(\partial_y \psi) \partial_x Q} - \widehat{(\partial_x \psi)(\partial_y Q)} - \left(\mu + \nu |𝐤|^{2n_\nu} \right) \widehat{q} + \widehat{F} . ``` +In doing so the Jacobian is computed in the conservative form: ``\mathsf{J}(f,g) = +\partial_y [ (\partial_x f) g] - \partial_x[ (\partial_y f) g]``. + The state variable `sol` is the Fourier transform of the sum of relative vorticity and vortex stretching (when the latter is applicable), [`qh`](@ref GeophysicalFlows.SingleLayerQG.Vars). -The Jacobian is computed in the conservative form: ``\mathsf{J}(f, g) = -\partial_y [ (\partial_x f) g] - \partial_x[ (\partial_y f) g]``. - The linear operator is constructed in `Equation` ```@docs @@ -122,4 +132,4 @@ Other diagnostic include: [`energy_dissipation`](@ref GeophysicalFlows.SingleLay (barotropic quasi-geostrophic flow with ``\beta=0``) above topography. - [`examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl`](@ref singlelayerqg_decaying_barotropic_equivalentbarotropic_example): - Simulate two dimensional turbulence (``\beta=0``) with both infinite and finite Rossby radius of deformation and compares the evolution of the two. \ No newline at end of file + Simulate two dimensional turbulence (``\beta=0``) with both infinite and finite Rossby radius of deformation and compares the evolution of the two. From 346f9146b7a770f9f7a2ff5498e5ee546d9de8fd Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 09:20:36 +1000 Subject: [PATCH 22/36] homogenize notation/wording --- docs/src/modules/multilayerqg.md | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/src/modules/multilayerqg.md b/docs/src/modules/multilayerqg.md index 96d0488b..d9717c8c 100644 --- a/docs/src/modules/multilayerqg.md +++ b/docs/src/modules/multilayerqg.md @@ -71,10 +71,14 @@ with ```math \begin{aligned} \partial_y Q_j &\equiv \beta - \partial_y^2 U_j - (1-\delta_{j,1}) F_{j-1/2, j} (U_{j-1} - U_j) - (1 - \delta_{j,n}) F_{j+1/2, j} (U_{j+1} - U_j) + \delta_{j,n} \partial_y \eta , \\ -\partial_x Q_j &\equiv \delta_{j, n} \partial_x \eta . +\partial_x Q_j & \equiv \delta_{j, n} \partial_x \eta , \end{aligned} ``` +the background PV gradient components in each layer and with +``\mathsf{J}(a, b) = (\partial_x a)(\partial_y b)-(\partial_y a)(\partial_x b)`` is the +two-dimensional Jacobian. On the right hand side, ``\mu`` is linear bottom drag, and ``\nu`` is +hyperviscosity of order ``n_\nu``. Plain old viscosity corresponds to ``n_\nu = 1``. ### Implementation From e6e77077478c0783bac31009b3977a69d46d730c Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 10:11:12 +1000 Subject: [PATCH 23/36] add empty line --- src/multilayerqg.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/multilayerqg.jl b/src/multilayerqg.jl index 05c41b0e..450b7027 100644 --- a/src/multilayerqg.jl +++ b/src/multilayerqg.jl @@ -803,6 +803,7 @@ function calcN_linearadvection!(N, sol, vars, params, grid) @. vars.vh = im * grid.kr * vars.ψh invtransform!(vars.u, vars.uh, params) + @. vars.u += params.U # add the imposed zonal flow U uQx, uQxh = vars.q, vars.uh # use vars.q and vars.uh as scratch variables @. uQx = vars.u * params.Qx # (U+u)*∂Q/∂x From 73399f09805c28d93c0888859b1b45dcf437e380 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 10:23:58 +1000 Subject: [PATCH 24/36] code alignment --- src/multilayerqg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/multilayerqg.jl b/src/multilayerqg.jl index 450b7027..99668a8b 100644 --- a/src/multilayerqg.jl +++ b/src/multilayerqg.jl @@ -757,7 +757,7 @@ function calcN_advection!(N, sol, vars, params, grid) @. vars.u += params.U # add the imposed zonal flow U uQx, uQxh = vars.q, vars.uh # use vars.q and vars.uh as scratch variables - @. uQx = vars.u * params.Qx # (U+u)*∂Q/∂x + @. uQx = vars.u * params.Qx # (U+u)*∂Q/∂x fwdtransform!(uQxh, uQx, params) @. N = - uQxh # -\hat{(U+u)*∂Q/∂x} From 99eb5ff7d17651d1478a11075f37802709041eed Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 10:24:39 +1000 Subject: [PATCH 25/36] add Qx and Qy in params; refactor calcN_advection --- src/singlelayerqg.jl | 121 ++++++++++++++++++++++++++++--------------- 1 file changed, 80 insertions(+), 41 deletions(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 24d6ac69..f468fdf4 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -109,7 +109,7 @@ function Problem(dev::Device=CPU(); if U isa Number U = convert(T, U) else - U = reshape(U, (1, grid.ny)) + U = repeat(reshape(U, (1, ny)), outer=(nx, 1)) # convert to 2D U = device_array(dev)(U) end @@ -136,7 +136,7 @@ The parameters for the `SingleLayerQG` problem. $(TYPEDFIELDS) """ -struct Params{T, Aphys, Atrans, Tℓ, TU <: Union{T, Aphys}} <: SingleLayerQGParams +struct Params{T, Aphys, Atrans, Tℓ, TU <: Union{T, Aphys}, Aphys2D} <: SingleLayerQGParams "planetary vorticity ``y``-gradient" β :: T "Rossby radius of deformation" @@ -155,6 +155,11 @@ struct Params{T, Aphys, Atrans, Tℓ, TU <: Union{T, Aphys}} <: SingleLayerQGPar nν :: Int "function that calculates the Fourier transform of the forcing, ``F̂``" calcF! :: Function + # derived params + "array containing ``x``-gradient of PV due to eta" + Qx :: Aphys2D + "array containing ``y``-gradient of PV due to ``U`` and topographic PV" + Qy :: Aphys2D end const BarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, Nothing} @@ -169,10 +174,27 @@ const SingleLayerQGvaryingUParams = Params{<:AbstractFloat, <:AbstractArray, <: Return the parameters for an Equivalent Barotropic QG problem (i.e., with finite Rossby radius of deformation). """ function EquivalentBarotropicQGParams(grid::AbstractGrid{T, A}, β, deformation_radius, U, eta, μ, ν, nν::Int, calcF) where {T, A} + + if U isa Number + U = convert(T, U) + elseif length(U) == grid.ny + U = repeat(reshape(U, (1, grid.ny)), outer=(grid.nx, 1)) # convert to 2D + U = A(U) + end + eta_on_grid = typeof(eta) <: AbstractArray ? A(eta) : FourierFlows.on_grid(eta, grid) - etah_on_grid = rfft(eta_on_grid) + etah = rfft(eta_on_grid) + + Qx = irfft(im * grid.kr .* etah, grid.nx) # ∂η/∂x + Qy = irfft(im * grid.l .* etah, grid.nx) # ∂η/∂y - return Params(β, deformation_radius, U, eta_on_grid, etah_on_grid, μ, ν, nν, calcF) + if U isa AbstractArray + Uh = rfft(U) + Uyy = irfft(- grid.l.^2 .* Uh, grid.nx) + Qy .-= Uyy # -∂²U/∂y² + end + + return Params(β, deformation_radius, U, eta_on_grid, etah, μ, ν, nν, calcF, Qx, Qy) end """ @@ -190,18 +212,21 @@ BarotropicQGParams(grid::AbstractGrid, β, U, eta, μ, ν, nν::Int, calcF) = """ Equation(params::BarotropicQGParams, grid) + Equation(params::EquivalentBarotropicQGParams, grid) -Return the equation for a barotropic QG problem with `params` and `grid`. Linear operator -``L`` includes bottom drag ``μ``, (hyper)-viscosity of order ``n_ν`` with coefficient ``ν``, -the ``β`` term, and a constant background flow ``U``: +Return the equation for a `SingleLayerQG` problem with `params` and `grid`. +Linear operator ``L`` includes bottom drag ``μ``, (hyper)-viscosity of order ``n_ν`` with +coefficient ``ν``, and the ``β`` term. If there is a constant background flow ``U`` that +does not vary in ``y`` then the linear term ``L`` includes also the mean advection term +by ``U``, namely ``-i k_x U```. That is: ```math -L = - μ - ν |𝐤|^{2 n_ν} + i β k_x / |𝐤|² - i k_x U . +L = -μ - ν |𝐤|^{2 n_ν} + i β k_x / (|𝐤|² + 1/ℓ²) - i k_x U . ``` The nonlinear term is computed via `calcN!` function. """ -function Equation(params::BarotropicQGParams, grid::AbstractGrid) +function Equation(params::BarotropicQGParams, grid) L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr * grid.invKrsq if params.U isa Number @@ -213,20 +238,7 @@ function Equation(params::BarotropicQGParams, grid::AbstractGrid) return FourierFlows.Equation(L, calcN!, grid) end -""" - Equation(params::EquivalentBarotropicQGParams, grid) - -Return the equation for an equivalent-barotropic QG problem with `params` and `grid`. -Linear operator ``L`` includes bottom drag ``μ``, (hyper)-viscosity of order ``n_ν`` with -coefficient ``ν``, the ``β`` term and a constant background flow ``U``: - -```math -L = -μ - ν |𝐤|^{2 n_ν} + i β k_x / (|𝐤|² + 1/ℓ²) - i k_x U . -``` - -The nonlinear term is computed via `calcN!` function. -""" -function Equation(params::EquivalentBarotropicQGParams, grid::AbstractGrid) +function Equation(params::EquivalentBarotropicQGParams, grid) L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr / (grid.Krsq + 1 / params.deformation_radius^2) if params.U isa Number @@ -329,16 +341,16 @@ end # ------- """ - calcN_advection!(N, sol, t, clock, vars, params, grid) + calcN_advection!(N, sol, t, clock, vars, params::SingleLayerQGconstantUParams, grid) -Calculate the Fourier transform of the advection term, ``- 𝖩(ψ, q+η) - U ∂(q+η)/∂x - v ∂U/∂y`` in conservative -form, i.e., ``- ∂[(u + U)(q + η - ∂U/∂y)]/∂x - ∂[v(q + η - ∂U/∂y)]/∂y`` and store it in `N`: +Compute the advection term and stores it in `N`. The imposed zonal flow ``U`` is either +zero or constant, in which case is incorporated in the linear terms of the equation. +Thus, the nonlinear terms is ``- 𝖩(ψ, q + η)`` in conservative form, i.e., +``- ∂_x[(∂_y ψ)(q+η)] - ∂_y[(∂_x ψ)(q+η)]``: ```math -N = - \\widehat{𝖩(ψ + Ψ, q + η - ∂U/∂y)} = - i k_x \\widehat{(u+U) (q + η - ∂U/∂y)} - i k_y \\widehat{v (q + η - ∂U/∂y)} . +N = - \\widehat{𝖩(ψ, q + η)} = - i k_x \\widehat{u (q + η)} - i k_y \\widehat{v (q + η)} . ``` - -Note: here ``- ∂Ψ/∂y = U``. """ function calcN_advection!(N, sol, t, clock, vars, params::SingleLayerQGconstantUParams, grid) @@ -363,33 +375,59 @@ function calcN_advection!(N, sol, t, clock, vars, params::SingleLayerQGconstantU @. N = -im * grid.kr * uq_plus_ηh - im * grid.l * vq_plus_ηh # - ∂[u*(q+η)]/∂x - ∂[v*(q+η)]/∂y + @. N -= im * grid.kr * params.U * params.etah + return nothing end +""" + calcN_advection!(N, sol, t, clock, vars, params::SingleLayerQGvaryingUParams, grid) + +Compute the advection term and stores it in `N`. The imposed zonal flow ``U(y)`` varies +with ``y`` and therefore is not taken care by the linear term `L` but rather is +incorporated in the nonlinear term `N`. + +```math +N = - \\widehat{𝖩(ψ, q)} - \\widehat{U ∂_x Q} - \\widehat{U ∂_x q} + + \\widehat{(∂_y ψ)(∂_x Q)} - \\widehat{(∂_x ψ)(∂_y Q)} . +``` +""" function calcN_advection!(N, sol, t, clock, vars, params::SingleLayerQGvaryingUParams, grid) @. vars.qh = sol + streamfunctionfrompv!(vars.ψh, vars.qh, params, grid) + @. vars.uh = -im * grid.l * vars.ψh @. vars.vh = im * grid.kr * vars.ψh - ldiv!(vars.q, grid.rfftplan, vars.qh) ldiv!(vars.u, grid.rfftplan, vars.uh) + @. vars.u += params.U # add the imposed zonal flow U + + uQx, uQxh = vars.q, vars.uh # use vars.q and vars.uh as scratch variables + @. uQx = vars.u * params.Qx # (U+u)*∂Q/∂x + + mul!(uQxh, grid.rfftplan, uQx) + + @. N = - uQxh # -\hat{(U+u)*∂Q/∂x} + ldiv!(vars.v, grid.rfftplan, vars.vh) - Uy = real.(ifft(im * grid.l .* fft(params.U))) # PV background (η - ∂U/∂y) + vQy, vQyh = vars.q, vars.vh # use vars.q and vars.vh as scratch variables - uq_plus_η = vars.u .+ params.U # use vars.u as scratch variable - @. uq_plus_η *= vars.q + params.eta .- Uy # (u + U) * (q + η - ∂U/∂y) - vq_plus_η = vars.v # use vars.v as scratch variable - @. vq_plus_η *= vars.q + params.eta .- Uy # v * (q + η - ∂U/∂y) + @. vQy = vars.v * params.Qy # v*∂Q/∂y + mul!(vQyh, grid.rfftplan, vQy) + @. N -= vQyh # -\hat{v*∂Q/∂y} - uq_plus_ηh = vars.uh # use vars.uh as scratch variable - mul!(uq_plus_ηh, grid.rfftplan, uq_plus_η) # \hat{(u + U) * (q + η - ∂U/∂y)} - vq_plus_ηh = vars.vh # use vars.vh as scratch variable - mul!(vq_plus_ηh, grid.rfftplan, vq_plus_η) # \hat{v * (q + η - ∂U/∂y)} + ldiv!(vars.q, grid.rfftplan, vars.qh) + + @. vars.u = params.U + Uq , Uqh = vars.u , vars.uh # use vars.u and vars.uh as scratch variables + @. Uq *= vars.q # U*q + + mul!(Uqh, grid.rfftplan, Uq) - @. N = -im * grid.kr * uq_plus_ηh - im * grid.l * vq_plus_ηh # - ∂[(u+U)*(q+η-∂U/∂y)]/∂x - ∂[v*(q+η-∂U/∂y)]/∂y + @. N -= im * grid.kr * Uqh # -\hat{∂[U*q]/∂x} return nothing end @@ -400,7 +438,8 @@ end Calculate the nonlinear term, that is the advection term and the forcing, ```math -N = - \\widehat{𝖩(ψ + X, q + η - ∂U/∂y)} + F̂ . +N = - \\widehat{𝖩(ψ, q)} - \\widehat{U ∂_x Q} - \\widehat{U ∂_x q} + + \\widehat{(∂_y ψ)(∂_x Q)} - \\widehat{(∂_x ψ)(∂_y Q)} + F̂ . ``` """ function calcN!(N, sol, t, clock, vars, params, grid) From 6604120826b0bbbc98a7e4049a3e09a323af9535 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 10:32:21 +1000 Subject: [PATCH 26/36] clarify when U=const --- docs/src/modules/singlelayerqg.md | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/docs/src/modules/singlelayerqg.md b/docs/src/modules/singlelayerqg.md index 1e63d9d5..bc55d3dd 100644 --- a/docs/src/modules/singlelayerqg.md +++ b/docs/src/modules/singlelayerqg.md @@ -41,13 +41,20 @@ the two-dimensional Jacobian. On the right hand side, ``F(x, y, t)`` is forcing, linear drag, and ``\nu`` is hyperviscosity of order ``n_\nu``. Plain old viscosity corresponds to ``n_\nu = 1``. +In the case that the imposed background zonal flow is just a constant, the above simplifies to: + +```math +\partial_t q + \mathsf{J}(\psi, q + \eta) + U \partial_x (q + \eta) + β \partial_x \psi = \underbrace{-\left[\mu + \nu(-1)^{n_\nu} \nabla^{2n_\nu} \right] q}_{\textrm{dissipation}} + F , +``` + +and thus the advection of ``q`` can be incorporated in the linear term ``L``. ### Implementation The equation is time-stepped forward in Fourier space: ```math -\partial_t \widehat{q} = - \widehat{\mathsf{J}(\psi, q)} - \widehat{U \partial_x Q} - \widehat{U \partial_x q} +\partial_t \widehat{q} = - \widehat{\mathsf{J}(\psi, q)} - \widehat{U \partial_x Q} - \widehat{U \partial_x q} + \widehat{(\partial_y \psi) \partial_x Q} - \widehat{(\partial_x \psi)(\partial_y Q)} - \left(\mu + \nu |𝐤|^{2n_\nu} \right) \widehat{q} + \widehat{F} . ``` From f3873335a5618dcee7a36c129084c235a05ee95a Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 15:41:11 +1000 Subject: [PATCH 27/36] some fixes + cleanup --- src/singlelayerqg.jl | 49 +++++++++++++++----------------- test/runtests.jl | 8 +++--- test/test_singlelayerqg.jl | 58 +++++++++++++++++++++++--------------- 3 files changed, 62 insertions(+), 53 deletions(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index f468fdf4..995d87fa 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -101,19 +101,15 @@ function Problem(dev::Device=CPU(); # the grid grid = TwoDGrid(dev; nx, Lx, ny, Ly, aliased_fraction, T) - x, y = gridpoints(grid) # topographic PV - eta === nothing && (eta = zeros(dev, T, (nx, ny))) + eta isa Nothing && (eta = zeros(dev, T, (nx, ny))) - if U isa Number - U = convert(T, U) - else - U = repeat(reshape(U, (1, ny)), outer=(nx, 1)) # convert to 2D - U = device_array(dev)(U) - end + U = U isa Number ? convert(T, U) : U - params = deformation_radius == Inf ? BarotropicQGParams(grid, T(β), U, eta, T(μ), T(ν), nν, calcF) : EquivalentBarotropicQGParams(grid, T(β), T(deformation_radius), U, eta, T(μ), T(ν), nν, calcF) + params = (deformation_radius == Inf || + deformation_radius === nothing) ? BarotropicQGParams(grid, T(β), U, eta, T(μ), T(ν), nν, calcF) : + EquivalentBarotropicQGParams(grid, T(deformation_radius), T(β), U, eta, T(μ), T(ν), nν, calcF) vars = calcF == nothingfunction ? DecayingVars(grid) : (stochastic ? StochasticForcedVars(grid) : ForcedVars(grid)) @@ -136,7 +132,7 @@ The parameters for the `SingleLayerQG` problem. $(TYPEDFIELDS) """ -struct Params{T, Aphys, Atrans, Tℓ, TU <: Union{T, Aphys}, Aphys2D} <: SingleLayerQGParams +struct Params{T, Aphys, Atrans, Tℓ, TU <: Union{T, Aphys}} <: SingleLayerQGParams "planetary vorticity ``y``-gradient" β :: T "Rossby radius of deformation" @@ -157,32 +153,30 @@ struct Params{T, Aphys, Atrans, Tℓ, TU <: Union{T, Aphys}, Aphys2D} <: SingleL calcF! :: Function # derived params "array containing ``x``-gradient of PV due to eta" - Qx :: Aphys2D + Qx :: Aphys "array containing ``y``-gradient of PV due to ``U`` and topographic PV" - Qy :: Aphys2D + Qy :: Aphys end -const BarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, Nothing} -const EquivalentBarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:AbstractFloat} +const BarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:Nothing, <:Any} +const EquivalentBarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:AbstractFloat, <:Any} -const SingleLayerQGconstantUParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:Any, <:AbstractFloat} -const SingleLayerQGvaryingUParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:Any, <:AbstractArray} +const SingleLayerQGconstantUParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:Any, <:Number} +const SingleLayerQGvaryingUParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:Any, <:AbstractArray} """ - EquivalentBarotropicQGParams(grid, β, deformation_radius, U, eta, μ, ν, nν, calcF) + EquivalentBarotropicQGParams(grid, deformation_radius, β, U, eta, μ, ν, nν, calcF) Return the parameters for an Equivalent Barotropic QG problem (i.e., with finite Rossby radius of deformation). """ -function EquivalentBarotropicQGParams(grid::AbstractGrid{T, A}, β, deformation_radius, U, eta, μ, ν, nν::Int, calcF) where {T, A} +function EquivalentBarotropicQGParams(grid::AbstractGrid{T, A}, deformation_radius, β, U, eta, μ, ν, nν::Int, calcF) where {T, A} - if U isa Number - U = convert(T, U) - elseif length(U) == grid.ny + if U isa AbstractArray && length(U) == grid.ny U = repeat(reshape(U, (1, grid.ny)), outer=(grid.nx, 1)) # convert to 2D U = A(U) end - eta_on_grid = typeof(eta) <: AbstractArray ? A(eta) : FourierFlows.on_grid(eta, grid) + eta_on_grid = eta isa AbstractArray ? A(eta) : FourierFlows.on_grid(eta, grid) etah = rfft(eta_on_grid) Qx = irfft(im * grid.kr .* etah, grid.nx) # ∂η/∂x @@ -190,8 +184,8 @@ function EquivalentBarotropicQGParams(grid::AbstractGrid{T, A}, β, deformation_ if U isa AbstractArray Uh = rfft(U) - Uyy = irfft(- grid.l.^2 .* Uh, grid.nx) - Qy .-= Uyy # -∂²U/∂y² + Uyy = irfft(- grid.l.^2 .* Uh, grid.nx) # ∂²U/∂y² + Qy .-= Uyy # -∂²U/∂y² end return Params(β, deformation_radius, U, eta_on_grid, etah, μ, ν, nν, calcF, Qx, Qy) @@ -202,8 +196,11 @@ end Return the parameters for a Barotropic QG problem (i.e., with infinite Rossby radius of deformation). """ -BarotropicQGParams(grid::AbstractGrid, β, U, eta, μ, ν, nν::Int, calcF) = - EquivalentBarotropicQGParams(grid, β, nothing, U, eta, μ, ν, nν, calcF) +function BarotropicQGParams(grid, β, U, eta, μ, ν, nν::Int, calcF) + deformation_radius = nothing + + return EquivalentBarotropicQGParams(grid, deformation_radius, β, U, eta, μ, ν, nν, calcF) +end # --------- diff --git a/test/runtests.jl b/test/runtests.jl index 70059ba4..e3123178 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -60,15 +60,15 @@ for dev in devices (20, 20, 20, 20, 200, 200, 2000, 2000,)) nx = 64 - @test test_1layerqg_rossbywave(timestepper, dt, nsteps, dev, nx; deformation_radius, U₀) - @test test_1layerqg_rossbywave(timestepper, dt, nsteps, dev, nx; deformation_radius, U₀=U₀*zeros(dev, Float64, (nx,))) - @test test_1layerqg_problemtype(dev, Float32; deformation_radius) + @test test_1layerqg_rossbywave(timestepper, dt, nsteps, dev, nx; deformation_radius, U=U₀) + @test test_1layerqg_rossbywave(timestepper, dt, nsteps, dev, nx; deformation_radius, U=U₀*ones((nx,))) end + @test test_1layerqg_problemtype(dev, Float32; deformation_radius, U=U₀) end @test test_1layerqg_advection(0.0005, "ForwardEuler", dev) @test test_streamfunctionfrompv(dev; deformation_radius=1.23) - @test test_1layerqg_energies_EquivalentBarotropicQG(dev; deformation_radius=1.23) @test test_1layerqg_energyenstrophy_BarotropicQG(dev) + @test test_1layerqg_energies_EquivalentBarotropicQG(dev; deformation_radius=1.23) @test test_1layerqg_deterministicforcing_energybudget(dev) @test test_1layerqg_stochasticforcing_energybudget(dev) @test test_1layerqg_deterministicforcing_enstrophybudget(dev) diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index f9554b17..4b42713e 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -3,7 +3,8 @@ Evolves a Rossby wave and compares with the analytic solution. """ -function test_1layerqg_rossbywave(stepper, dt, nsteps, dev::Device=CPU(), nx=64; deformation_radius=Inf, U₀ = 0) +function test_1layerqg_rossbywave(stepper, dt, nsteps, dev::Device=CPU(), nx=64; deformation_radius=Inf, U = 0) + Lx = 2π β = 2.0 μ = 0.0 @@ -18,7 +19,7 @@ function test_1layerqg_rossbywave(stepper, dt, nsteps, dev::Device=CPU(), nx=64; eta(x, y) = 0 * x end - prob = SingleLayerQG.Problem(dev; nx=nx, Lx=Lx, eta=eta, deformation_radius=deformation_radius, β=β, U=U₀, μ=μ, ν=ν, stepper=stepper, dt=dt) + prob = SingleLayerQG.Problem(dev; nx=nx, Lx=Lx, eta=eta, deformation_radius=deformation_radius, β=β, U=U, μ=μ, ν=ν, stepper=stepper, dt=dt) sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid x, y = gridpoints(grid) @@ -37,7 +38,7 @@ function test_1layerqg_rossbywave(stepper, dt, nsteps, dev::Device=CPU(), nx=64; dealias!(sol, grid) SingleLayerQG.updatevars!(prob) - q_theory = @. ampl * cos(kwave * (x - (ω / kwave + U₀) * clock.t)) * cos(lwave * y) + q_theory = @CUDA.allowscalar @. ampl * cos(kwave * (x - (U[1] + ω / kwave) * clock.t)) * cos(lwave * y) return isapprox(q_theory, vars.q, rtol=grid.nx * grid.ny * nsteps * 1e-12) end @@ -308,50 +309,59 @@ Tests the energy and enstrophy function for a SingleLayerQG problem. function test_1layerqg_energyenstrophy_BarotropicQG(dev::Device=CPU()) nx, Lx = 64, 2π ny, Ly = 64, 3π - grid = TwoDGrid(dev; nx, Lx, ny, Ly) - k₀, l₀ = 2π/Lx, 2π/Ly # fundamental wavenumbers + + prob = SingleLayerQG.Problem(dev; nx=nx, Lx=Lx, ny=ny, Ly=Ly, stepper="ForwardEuler") + grid = prob.grid + + k₀, l₀ = 2π / grid.Lx, 2π / grid.Ly # fundamental wavenumbers x, y = gridpoints(grid) energy_calc = 29/9 enstrophy_calc = 10885/648 - η = @. cos(10k₀ * x) * cos(10l₀ * y) + η(x, y) = cos(10k₀ * x) * cos(10l₀ * y) ψ₀ = @. sin(2k₀ * x) * cos(2l₀ * y) + 2sin(k₀ * x) * cos(3l₀ * y) q₀ = @. - ((2k₀)^2 + (2l₀)^2) * sin(2k₀ * x) * cos(2l₀ * y) - (k₀^2 + (3l₀)^2) * 2sin(k₀ * x) * cos(3l₀*y) prob = SingleLayerQG.Problem(dev; nx=nx, Lx=Lx, ny=ny, Ly=Ly, eta=η, stepper="ForwardEuler") sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid + SingleLayerQG.set_q!(prob, q₀) SingleLayerQG.updatevars!(prob) energyq₀ = SingleLayerQG.energy(prob) enstrophyq₀ = SingleLayerQG.enstrophy(prob) - return isapprox(energyq₀, energy_calc, rtol=rtol_singlelayerqg) && isapprox(enstrophyq₀, enstrophy_calc, rtol=rtol_singlelayerqg) && SingleLayerQG.potential_energy(prob)==0 && - SingleLayerQG.addforcing!(prob.timestepper.N, sol, clock.t, clock, vars, params, grid) === nothing + return (isapprox(energyq₀, energy_calc, rtol=rtol_singlelayerqg) && + isapprox(enstrophyq₀, enstrophy_calc, rtol=rtol_singlelayerqg) && + SingleLayerQG.potential_energy(prob)==0 && + SingleLayerQG.addforcing!(prob.timestepper.N, sol, clock.t, clock, vars, params, grid) === nothing) end """ - test_1layerqg_energies_EquivalentBarotropicQG(dev) + test_1layerqg_energies_EquivalentBarotropicQG(dev; deformation_radius=1.23) Tests the kinetic and potential energy for an equivalent barotropic SingleLayerQG problem. """ function test_1layerqg_energies_EquivalentBarotropicQG(dev; deformation_radius=1.23) nx, Lx = 64, 2π ny, Ly = 64, 3π - grid = TwoDGrid(dev; nx, Lx, ny, Ly) - k₀, l₀ = 2π/Lx, 2π/Ly # fundamental wavenumbers + + prob = SingleLayerQG.Problem(dev; nx, Lx, ny, Ly, deformation_radius, stepper="ForwardEuler") + grid = prob.grid + + k₀, l₀ = 2π / grid.Lx, 2π / grid.Ly # fundamental wavenumbers x, y = gridpoints(grid) kinetic_energy_calc = 29/9 potential_energy_calc = 5/(8*deformation_radius^2) energy_calc = kinetic_energy_calc + potential_energy_calc - - η = @. cos(10k₀ * x) * cos(10l₀ * y) + + η(x, y) = cos(10k₀ * x) * cos(10l₀ * y) ψ₀ = @. sin(2k₀ * x) * cos(2l₀ * y) + 2sin(k₀ * x) * cos(3l₀ * y) q₀ = @. - ((2k₀)^2 + (2l₀)^2) * sin(2k₀ * x) * cos(2l₀ * y) - (k₀^2 + (3l₀)^2) * 2sin(k₀ * x) * cos(3l₀*y) - 1/deformation_radius^2 * ψ₀ - prob = SingleLayerQG.Problem(dev; nx=nx, Lx=Lx, ny=ny, Ly=Ly, eta=η, deformation_radius=deformation_radius, stepper="ForwardEuler") + prob = SingleLayerQG.Problem(dev; nx, Lx, ny, Ly, eta=η, deformation_radius, stepper="ForwardEuler") sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid SingleLayerQG.set_q!(prob, q₀) SingleLayerQG.updatevars!(prob) @@ -360,21 +370,24 @@ function test_1layerqg_energies_EquivalentBarotropicQG(dev; deformation_radius=1 potential_energyq₀ = SingleLayerQG.potential_energy(prob) energyq₀ = SingleLayerQG.energy(prob) - return isapprox(kinetic_energyq₀, kinetic_energy_calc, rtol=rtol_singlelayerqg) && isapprox(potential_energyq₀, potential_energy_calc, rtol=rtol_singlelayerqg) && isapprox(energyq₀, energy_calc, rtol=rtol_singlelayerqg) && - SingleLayerQG.addforcing!(prob.timestepper.N, sol, clock.t, clock, vars, params, grid) === nothing + return (isapprox(kinetic_energyq₀, kinetic_energy_calc, rtol=rtol_singlelayerqg) && + isapprox(potential_energyq₀, potential_energy_calc, rtol=rtol_singlelayerqg) && + isapprox(energyq₀, energy_calc, rtol=rtol_singlelayerqg) && + SingleLayerQG.addforcing!(prob.timestepper.N, sol, clock.t, clock, vars, params, grid) === nothing) end """ - test_1layerqg_problemtype(dev, T; deformation_radius=Inf) + test_1layerqg_problemtype(dev, T; deformation_radius=Inf, U=0) -Tests the SingleLayerQG problem constructor for different DataType `T`. +Test the SingleLayerQG problem constructor for different DataType `T`. """ -function test_1layerqg_problemtype(dev, T; deformation_radius=Inf) - prob = SingleLayerQG.Problem(dev; T=T, deformation_radius=deformation_radius) +function test_1layerqg_problemtype(dev, T; deformation_radius=Inf, U=0) + prob = SingleLayerQG.Problem(dev; T, deformation_radius, U) A = device_array(dev) - - return (typeof(prob.sol)<:A{Complex{T}, 2} && typeof(prob.grid.Lx)==T && eltype(prob.grid.x)==T && typeof(prob.vars.u)<:A{T, 2}) + + return (typeof(prob.sol)<:A{Complex{T}, 2} && typeof(prob.grid.Lx)==T && + eltype(prob.grid.x)==T && typeof(prob.vars.u)<:A{T, 2} && typeof(prob.params.U)==T) end function test_streamfunctionfrompv(dev; deformation_radius=1.23) @@ -394,7 +407,6 @@ function test_streamfunctionfrompv(dev; deformation_radius=1.23) SingleLayerQG.set_q!(prob_equivalentbarotropicQG, q_equivalentbarotropic) SingleLayerQG.streamfunctionfrompv!(prob_barotropicQG.vars.ψh, prob_barotropicQG.vars.qh, prob_barotropicQG.params, prob_barotropicQG.grid) - SingleLayerQG.streamfunctionfrompv!(prob_equivalentbarotropicQG.vars.ψh, prob_equivalentbarotropicQG.vars.qh, prob_equivalentbarotropicQG.params, prob_equivalentbarotropicQG.grid) return (prob_barotropicQG.vars.ψ ≈ ψ && prob_equivalentbarotropicQG.vars.ψ ≈ ψ) From 40cc4db6a722ded19cd1bb7d035477a02c8d27b5 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 16:24:16 +1000 Subject: [PATCH 28/36] add nonlinear advection test for finite deformation + topo --- test/runtests.jl | 3 ++- test/test_singlelayerqg.jl | 49 +++++++++++++++++++++++++++++++++++++- 2 files changed, 50 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index e3123178..1e728cef 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -65,7 +65,8 @@ for dev in devices end @test test_1layerqg_problemtype(dev, Float32; deformation_radius, U=U₀) end - @test test_1layerqg_advection(0.0005, "ForwardEuler", dev) + @test test_1layerqg_nonlinearadvection(0.0005, "ForwardEuler", dev) + @test test_1layerqg_nonlinearadvection_deformation(0.0005, "ForwardEuler", dev) @test test_streamfunctionfrompv(dev; deformation_radius=1.23) @test test_1layerqg_energyenstrophy_BarotropicQG(dev) @test test_1layerqg_energies_EquivalentBarotropicQG(dev; deformation_radius=1.23) diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index 4b42713e..4c5616c3 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -265,7 +265,7 @@ forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - ν∇²ζf. to the vorticity equation forced by this Ff is then ζf. (This solution may not be realized, at least at long times, if it is unstable.) """ -function test_1layerqg_advection(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1, μ=0.0) +function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1, μ=0.0) n, L = 128, 2π ν, nν = 1e-2, 1 μ = 0.0 @@ -301,6 +301,53 @@ function test_1layerqg_advection(dt, stepper, dev::Device=CPU(); n=128, L=2π, return isapprox(prob.vars.q, qf, rtol=rtol_singlelayerqg) end +""" + test_1layerqg_nonlinearadvection_deformation(dt, stepper, dev; kwargs...) + +Same as `test_1layerqg_advection` but with finite deformation radius. +""" +function test_1layerqg_nonlinearadvection_deformation(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1, μ=0.0) + n, L = 128, 2π + ν, nν = 1e-2, 1 + μ = 0.0 + tf = 1.0 + nt = round(Int, tf/dt) + + grid = TwoDGrid(dev; nx=n, Lx=L) + k₀, l₀ = 2π / grid.Lx, 2π / grid.Ly # fundamental wavenumbers + x, y = gridpoints(grid) + + deformation_radius = 1.23 + + η₀ = 0.4 + η(x, y) = η₀ * cos(10x) * cos(10y) + ψf = @. sin(2x) * cos(2y) + 2sin(x) * cos(3y) + qf = @. - (2^2 + 2^2) * sin(2x) * cos(2y) - (1^2 + 3^2) * 2sin(x) * cos(3y) - 1/deformation_radius^2 * ψf + + Ff = @. (- ν*(64sin(2x) * cos(2y) + 200sin(x) * cos(3y) + + (8sin(2x) * cos(2y) + 20sin(x) * cos(3y)) / deformation_radius^2) + + 4sin(x) * (sin(y) - sin(5y) + 2cos(2x) * (2sin(y) + sin(5y))) + - 20η₀ * (cos(10y) * sin(10x) * (sin(2x) * sin(2y) + 3sin(x) * sin(3y)) + + cos(10x) * sin(10y) * (cos(2x) * cos(2y) + cos(x) * cos(3y)))) + + Ffh = rfft(Ff) + + function calcF!(Fh, sol, t, clock, vars, params, grid) + Fh .= Ffh + return nothing + end + + prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, eta=η, ν=ν, nν=nν, μ=μ, dt=dt, deformation_radius=deformation_radius, stepper=stepper, calcF=calcF!) + + SingleLayerQG.set_q!(prob, qf) + + stepforward!(prob, round(Int, nt)) + + SingleLayerQG.updatevars!(prob) + + return isapprox(prob.vars.q, qf, rtol=rtol_singlelayerqg) +end + """ test_1layerqg_energyenstrophy_BarotropicQG(dev) From 63a7a8c8eff9bf21ba6d1904c9c9f52db01167bd Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 17:53:29 +1000 Subject: [PATCH 29/36] fixes nonlinear advection w U(y) + add tests --- src/singlelayerqg.jl | 19 ++++++++++--------- test/runtests.jl | 1 + test/test_singlelayerqg.jl | 28 +++++++++++++++++++++++++--- 3 files changed, 36 insertions(+), 12 deletions(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 995d87fa..684d51aa 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -1,5 +1,7 @@ module SingleLayerQG +using LinearAlgebra + export Problem, streamfunctionfrompv!, @@ -372,7 +374,7 @@ function calcN_advection!(N, sol, t, clock, vars, params::SingleLayerQGconstantU @. N = -im * grid.kr * uq_plus_ηh - im * grid.l * vq_plus_ηh # - ∂[u*(q+η)]/∂x - ∂[v*(q+η)]/∂y - @. N -= im * grid.kr * params.U * params.etah + @. N -= im * grid.kr * params.U * params.etah # - \hat{∂(Uη)/∂x} return nothing end @@ -403,28 +405,27 @@ function calcN_advection!(N, sol, t, clock, vars, params::SingleLayerQGvaryingUP uQx, uQxh = vars.q, vars.uh # use vars.q and vars.uh as scratch variables @. uQx = vars.u * params.Qx # (U+u)*∂Q/∂x - mul!(uQxh, grid.rfftplan, uQx) - @. N = - uQxh # -\hat{(U+u)*∂Q/∂x} ldiv!(vars.v, grid.rfftplan, vars.vh) vQy, vQyh = vars.q, vars.vh # use vars.q and vars.vh as scratch variables - @. vQy = vars.v * params.Qy # v*∂Q/∂y mul!(vQyh, grid.rfftplan, vQy) @. N -= vQyh # -\hat{v*∂Q/∂y} ldiv!(vars.q, grid.rfftplan, vars.qh) - @. vars.u = params.U - Uq , Uqh = vars.u , vars.uh # use vars.u and vars.uh as scratch variables - @. Uq *= vars.q # U*q + uq , vq = vars.u , vars.v # use vars.u and vars.v as scratch variables + uqh, vqh = vars.uh, vars.vh # use vars.uh and vars.vh as scratch variables + @. uq *= vars.q # (U+u)*q + @. vq *= vars.q # v*q - mul!(Uqh, grid.rfftplan, Uq) + mul!(uqh, grid.rfftplan, uq) + mul!(vqh, grid.rfftplan, vq) - @. N -= im * grid.kr * Uqh # -\hat{∂[U*q]/∂x} + @. N -= im * grid.kr * uqh + im * grid.l * vqh # -\hat{∂[(U+u)q]/∂x} - \hat{∂[vq]/∂y} return nothing end diff --git a/test/runtests.jl b/test/runtests.jl index 1e728cef..fd75a7f8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -66,6 +66,7 @@ for dev in devices @test test_1layerqg_problemtype(dev, Float32; deformation_radius, U=U₀) end @test test_1layerqg_nonlinearadvection(0.0005, "ForwardEuler", dev) + @test test_1layerqg_nonlinearadvection(0.0005, "ForwardEuler", dev, add_background_flow = true) @test test_1layerqg_nonlinearadvection_deformation(0.0005, "ForwardEuler", dev) @test test_streamfunctionfrompv(dev; deformation_radius=1.23) @test test_1layerqg_energyenstrophy_BarotropicQG(dev) diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index 4c5616c3..a7cb43a6 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -1,3 +1,5 @@ +using LinearAlgebra + """ test_1layerqg_rossbywave(stepper, dt, nsteps, dev; kwargs...) @@ -265,7 +267,7 @@ forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - ν∇²ζf. to the vorticity equation forced by this Ff is then ζf. (This solution may not be realized, at least at long times, if it is unstable.) """ -function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1, μ=0.0) +function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1, μ=0.0, add_background_flow = false) n, L = 128, 2π ν, nν = 1e-2, 1 μ = 0.0 @@ -283,6 +285,16 @@ function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); n=128, + 8 * ( cos(x) * cos(3y) * sin(2x) * sin(2y) - 3cos(2x) * cos(2y) * sin(x) * sin(3y) ) ) + if add_background_flow == true + U = @. 0.2 / cosh(2y)^2 + U = device_array(dev)(U) + Uh = rfft(U) + Uyy = irfft(- grid.l.^2 .* Uh, grid.nx) # ∂²U/∂y² + + @. Ff -= ( U * (16cos(2x) * cos(2y) + 20cos(x) * cos(3y)) + + Uyy * (2cos(2x) * cos(2y) + 2cos(x) * cos(3y))) + end + Ffh = rfft(Ff) function calcF!(Fh, sol, t, clock, vars, params, grid) @@ -290,14 +302,24 @@ function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); n=128, return nothing end - prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt, stepper=stepper, calcF=calcF!) + if add_background_flow == false + U₀ = 0 + elseif add_background_flow==true + U₀ = U[1, :] + end + + prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, U = U₀, ν=ν, nν=nν, μ=μ, dt=dt, stepper=stepper, calcF=calcF!) SingleLayerQG.set_q!(prob, qf) stepforward!(prob, round(Int, nt)) SingleLayerQG.updatevars!(prob) - + + println(norm(prob.vars.q)) + println(norm(qf)) + println(norm(prob.vars.q - qf)) + return isapprox(prob.vars.q, qf, rtol=rtol_singlelayerqg) end From 6430e59f5531ce4adf1176345e5dbfdffaf86344 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 18:33:51 +1000 Subject: [PATCH 30/36] =?UTF-8?q?add=20comment=20for=20=CE=B2=20term?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/singlelayerqg.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 684d51aa..6454bef0 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -190,6 +190,8 @@ function EquivalentBarotropicQGParams(grid::AbstractGrid{T, A}, deformation_radi Qy .-= Uyy # -∂²U/∂y² end + # Note: The β-term in Qy is included in the linear term L of Equation. + return Params(β, deformation_radius, U, eta_on_grid, etah, μ, ν, nν, calcF, Qx, Qy) end From 48f2f1b5ab72c6e980c26b40f803a34ad374b779 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 19:01:05 +1000 Subject: [PATCH 31/36] remove debug statements + fix test for julia 1.6 --- src/singlelayerqg.jl | 2 +- test/test_singlelayerqg.jl | 4 ---- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 6454bef0..e56f0837 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -134,7 +134,7 @@ The parameters for the `SingleLayerQG` problem. $(TYPEDFIELDS) """ -struct Params{T, Aphys, Atrans, Tℓ, TU <: Union{T, Aphys}} <: SingleLayerQGParams +struct Params{T, Aphys, Atrans, Tℓ, TU} <: SingleLayerQGParams "planetary vorticity ``y``-gradient" β :: T "Rossby radius of deformation" diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index a7cb43a6..3fc3d9cd 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -316,10 +316,6 @@ function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); n=128, SingleLayerQG.updatevars!(prob) - println(norm(prob.vars.q)) - println(norm(qf)) - println(norm(prob.vars.q - qf)) - return isapprox(prob.vars.q, qf, rtol=rtol_singlelayerqg) end From 42d90030b9b1abd3fafeebe615bad2a06cbe31bf Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 21:13:53 +1000 Subject: [PATCH 32/36] more tests for the nonlinear terms --- test/runtests.jl | 7 +++- test/test_singlelayerqg.jl | 83 ++++++++++++++++++++++++++------------ 2 files changed, 62 insertions(+), 28 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index fd75a7f8..95e193c1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -65,8 +65,11 @@ for dev in devices end @test test_1layerqg_problemtype(dev, Float32; deformation_radius, U=U₀) end - @test test_1layerqg_nonlinearadvection(0.0005, "ForwardEuler", dev) - @test test_1layerqg_nonlinearadvection(0.0005, "ForwardEuler", dev, add_background_flow = true) + @test test_1layerqg_nonlinearadvection(0.0005, "ForwardEuler", dev, add_background_flow = false, add_topography = false) + @test test_1layerqg_nonlinearadvection(0.0005, "ForwardEuler", dev, add_background_flow = false, add_topography = true) + @test test_1layerqg_nonlinearadvection(0.0005, "ForwardEuler", dev, add_background_flow = true, background_flow_vary_in_y = false) + @test test_1layerqg_nonlinearadvection(0.0005, "ForwardEuler", dev, add_background_flow = true, background_flow_vary_in_y = false, add_topography = true) + @test test_1layerqg_nonlinearadvection(0.0005, "ForwardEuler", dev, add_background_flow = true, background_flow_vary_in_y = true) @test test_1layerqg_nonlinearadvection_deformation(0.0005, "ForwardEuler", dev) @test test_streamfunctionfrompv(dev; deformation_radius=1.23) @test test_1layerqg_energyenstrophy_BarotropicQG(dev) diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index 3fc3d9cd..d70c34da 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -257,17 +257,26 @@ function test_1layerqg_deterministicforcing_enstrophybudget(dev::Device=CPU(); n end """ - test_1layerqg_nonlinearadvection(dt, stepper, dev; kwargs...) + test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); + add_topography = false, + add_background_flow = false, + background_flow_vary_in_y = true) -Tests the advection term in the SingleLayerQG module by timestepping a +Tests the advection term in the `SingleLayerQG` module by timestepping a test problem with timestep dt and timestepper identified by the string stepper. -The test problem is derived by picking a solution ζf (with associated -streamfunction ψf) for which the advection term J(ψf, ζf) is non-zero. Next, a -forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - ν∇²ζf. One solution -to the vorticity equation forced by this Ff is then ζf. (This solution may not +The test problem is derived by picking a solution qf (with associated +streamfunction ψf) for which the advection term J(qf, ζf) is non-zero. Next, a +forcing Ff is derived according to Ff = ∂qf/∂t + J(ψf, qf) - ν∇²qf. One solution +to the vorticity equation forced by this Ff is then qf. (This solution may not be realized, at least at long times, if it is unstable.) + +We can optionally add an imposed mean flow `U(y)` via the kwarg `add_background_flow = true` +or add topography via kwarg `add_topography = true`. """ -function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1, μ=0.0, add_background_flow = false) +function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); + add_topography = false, + add_background_flow = false, + background_flow_vary_in_y = true) n, L = 128, 2π ν, nν = 1e-2, 1 μ = 0.0 @@ -277,24 +286,42 @@ function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); n=128, grid = TwoDGrid(dev; nx=n, Lx=L) x, y = gridpoints(grid) + if add_topography == true + η₀ = 0.4 + elseif add_topography == false + η₀ = 0 + end + + η(x, y) = η₀ * cos(10x) * cos(10y) + ψf = @. sin(2x) * cos(2y) + 2sin(x) * cos(3y) qf = @. -8sin(2x) * cos(2y) - 20sin(x) * cos(3y) - Ff = @. -( - ν*( 64sin(2x) * cos(2y) + 200sin(x) * cos(3y) ) - + 8 * ( cos(x) * cos(3y) * sin(2x) * sin(2y) - 3cos(2x) * cos(2y) * sin(x) * sin(3y) ) - ) + # F = J(ψ, q + η) - ν ∇²q + U ∂(q+η)/∂x - v ∂²U/∂y² + # where q = ∇²ψ + Ff = @. (- ν * ( 64sin(2x) * cos(2y) + 200sin(x) * cos(3y) ) + - 8 * (cos(x) * cos(3y) * sin(2x) * sin(2y) - 3cos(2x) * cos(2y) * sin(x) * sin(3y)) + - 20η₀ * (sin(10x) * cos(10y) * (sin(2x) * sin(2y) + 3sin(x) * sin(3y)) + + cos(10x) * sin(10y) * (cos(2x) * cos(2y) + cos(x) * cos(3y)))) - if add_background_flow == true - U = @. 0.2 / cosh(2y)^2 + U_amplitude = 0.2 + + U = 0 + Uyy = 0 + + if add_background_flow == true && background_flow_vary_in_y == false + U = U_amplitude + elseif add_background_flow == true && background_flow_vary_in_y == true + U = @. U_amplitude / cosh(2y)^2 U = device_array(dev)(U) Uh = rfft(U) Uyy = irfft(- grid.l.^2 .* Uh, grid.nx) # ∂²U/∂y² - - @. Ff -= ( U * (16cos(2x) * cos(2y) + 20cos(x) * cos(3y)) - + Uyy * (2cos(2x) * cos(2y) + 2cos(x) * cos(3y))) end + @. Ff += (- U * (16cos(2x) * cos(2y) + 20cos(x) * cos(3y)) + - Uyy * (2cos(2x) * cos(2y) + 2cos(x) * cos(3y)) + - η₀ * U * 10sin(10x) * cos(10y)) + Ffh = rfft(Ff) function calcF!(Fh, sol, t, clock, vars, params, grid) @@ -304,11 +331,13 @@ function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); n=128, if add_background_flow == false U₀ = 0 - elseif add_background_flow==true + elseif add_background_flow == true && background_flow_vary_in_y == true U₀ = U[1, :] + elseif add_background_flow == true && background_flow_vary_in_y == false + U₀ = U_amplitude end - prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, U = U₀, ν=ν, nν=nν, μ=μ, dt=dt, stepper=stepper, calcF=calcF!) + prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, U = U₀, eta = η, ν=ν, nν=nν, μ=μ, dt=dt, stepper=stepper, calcF=calcF!) SingleLayerQG.set_q!(prob, qf) @@ -320,11 +349,13 @@ function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); n=128, end """ - test_1layerqg_nonlinearadvection_deformation(dt, stepper, dev; kwargs...) + test_1layerqg_nonlinearadvection_deformation(dt, stepper, dev::Device=CPU(); + deformation_radius=1.23) -Same as `test_1layerqg_advection` but with finite deformation radius. +Same as `test_1layerqg_nonlinearadvection` but with finite deformation radius. """ -function test_1layerqg_nonlinearadvection_deformation(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1, μ=0.0) +function test_1layerqg_nonlinearadvection_deformation(dt, stepper, dev::Device=CPU(); + deformation_radius=1.23) n, L = 128, 2π ν, nν = 1e-2, 1 μ = 0.0 @@ -335,16 +366,16 @@ function test_1layerqg_nonlinearadvection_deformation(dt, stepper, dev::Device=C k₀, l₀ = 2π / grid.Lx, 2π / grid.Ly # fundamental wavenumbers x, y = gridpoints(grid) - deformation_radius = 1.23 - η₀ = 0.4 η(x, y) = η₀ * cos(10x) * cos(10y) ψf = @. sin(2x) * cos(2y) + 2sin(x) * cos(3y) qf = @. - (2^2 + 2^2) * sin(2x) * cos(2y) - (1^2 + 3^2) * 2sin(x) * cos(3y) - 1/deformation_radius^2 * ψf - Ff = @. (- ν*(64sin(2x) * cos(2y) + 200sin(x) * cos(3y) - + (8sin(2x) * cos(2y) + 20sin(x) * cos(3y)) / deformation_radius^2) - + 4sin(x) * (sin(y) - sin(5y) + 2cos(2x) * (2sin(y) + sin(5y))) + # F = J(ψ, q + η) - ν ∇²q + # where q = ∇²ψ - ψ/ℓ² + Ff = @. (- ν * (64sin(2x) * cos(2y) + 200sin(x) * cos(3y) + + (8sin(2x) * cos(2y) + 20sin(x) * cos(3y)) / deformation_radius^2) + - 8 * (cos(x) * cos(3y) * sin(2x) * sin(2y) - 3cos(2x) * cos(2y) * sin(x) * sin(3y)) - 20η₀ * (cos(10y) * sin(10x) * (sin(2x) * sin(2y) + 3sin(x) * sin(3y)) + cos(10x) * sin(10y) * (cos(2x) * cos(2y) + cos(x) * cos(3y)))) From 22e3bde6601afda29d573e779e157ab039bffab3 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 21:35:46 +1000 Subject: [PATCH 33/36] fix test_1layerqg_nonlinearadvection on gpu --- test/test_singlelayerqg.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index d70c34da..8f6ffdfb 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -292,7 +292,8 @@ function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); η₀ = 0 end - η(x, y) = η₀ * cos(10x) * cos(10y) + η(x, y; η₀ = η₀) = η₀ * cos(10x) * cos(10y) + η_array = @. η(x, y; η₀) ψf = @. sin(2x) * cos(2y) + 2sin(x) * cos(3y) qf = @. -8sin(2x) * cos(2y) - 20sin(x) * cos(3y) @@ -337,7 +338,7 @@ function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); U₀ = U_amplitude end - prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, U = U₀, eta = η, ν=ν, nν=nν, μ=μ, dt=dt, stepper=stepper, calcF=calcF!) + prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, U = U₀, eta = η_array, ν=ν, nν=nν, μ=μ, dt=dt, stepper=stepper, calcF=calcF!) SingleLayerQG.set_q!(prob, qf) From c6cf36a27910e84833720b89ae344716200657fb Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 21:47:23 +1000 Subject: [PATCH 34/36] fix test_1layerqg_nonlinearadvection on gpu --- test/test_singlelayerqg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_singlelayerqg.jl b/test/test_singlelayerqg.jl index 8f6ffdfb..74289e03 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -292,8 +292,8 @@ function test_1layerqg_nonlinearadvection(dt, stepper, dev::Device=CPU(); η₀ = 0 end - η(x, y; η₀ = η₀) = η₀ * cos(10x) * cos(10y) - η_array = @. η(x, y; η₀) + η(x, y, η₀) = η₀ * cos(10x) * cos(10y) + η_array = @. η(x, y, η₀) ψf = @. sin(2x) * cos(2y) + 2sin(x) * cos(3y) qf = @. -8sin(2x) * cos(2y) - 20sin(x) * cos(3y) From 5c9eb93868fbbd63b634853280891580929b1907 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Jul 2024 22:14:34 +1000 Subject: [PATCH 35/36] improve docs --- docs/Project.toml | 1 + docs/src/modules/singlelayerqg.md | 2 +- src/singlelayerqg.jl | 8 +++----- 3 files changed, 5 insertions(+), 6 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index cf889bcc..dfbe88eb 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -3,6 +3,7 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" +GeophysicalFlows = "44ee3b1c-bc02-53fa-8355-8e347616e15e" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" diff --git a/docs/src/modules/singlelayerqg.md b/docs/src/modules/singlelayerqg.md index bc55d3dd..e8e09603 100644 --- a/docs/src/modules/singlelayerqg.md +++ b/docs/src/modules/singlelayerqg.md @@ -55,7 +55,7 @@ The equation is time-stepped forward in Fourier space: ```math \partial_t \widehat{q} = - \widehat{\mathsf{J}(\psi, q)} - \widehat{U \partial_x Q} - \widehat{U \partial_x q} -+ \widehat{(\partial_y \psi) \partial_x Q} - \widehat{(\partial_x \psi)(\partial_y Q)} - \left(\mu + \nu |𝐤|^{2n_\nu} \right) \widehat{q} + \widehat{F} . ++ \widehat{(\partial_y \psi) (\partial_x Q)} - \widehat{(\partial_x \psi)(\partial_y Q)} - \left(\mu + \nu |𝐤|^{2n_\nu} \right) \widehat{q} + \widehat{F} . ``` In doing so the Jacobian is computed in the conservative form: ``\mathsf{J}(f,g) = diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index e56f0837..baaf26e0 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -347,7 +347,7 @@ end Compute the advection term and stores it in `N`. The imposed zonal flow ``U`` is either zero or constant, in which case is incorporated in the linear terms of the equation. Thus, the nonlinear terms is ``- 𝖩(ψ, q + η)`` in conservative form, i.e., -``- ∂_x[(∂_y ψ)(q+η)] - ∂_y[(∂_x ψ)(q+η)]``: +``- ∂_x[(∂_y ψ)(q + η)] - ∂_y[(∂_x ψ)(q + η)]``: ```math N = - \\widehat{𝖩(ψ, q + η)} = - i k_x \\widehat{u (q + η)} - i k_y \\widehat{v (q + η)} . @@ -389,8 +389,7 @@ with ``y`` and therefore is not taken care by the linear term `L` but rather is incorporated in the nonlinear term `N`. ```math -N = - \\widehat{𝖩(ψ, q)} - \\widehat{U ∂_x Q} - \\widehat{U ∂_x q} - + \\widehat{(∂_y ψ)(∂_x Q)} - \\widehat{(∂_x ψ)(∂_y Q)} . +N = - \\widehat{𝖩(ψ, q + η)} - \\widehat{U ∂_x (q + η)} + \\widehat{(∂_x ψ)(∂_y² U)} . ``` """ function calcN_advection!(N, sol, t, clock, vars, params::SingleLayerQGvaryingUParams, grid) @@ -438,8 +437,7 @@ end Calculate the nonlinear term, that is the advection term and the forcing, ```math -N = - \\widehat{𝖩(ψ, q)} - \\widehat{U ∂_x Q} - \\widehat{U ∂_x q} - + \\widehat{(∂_y ψ)(∂_x Q)} - \\widehat{(∂_x ψ)(∂_y Q)} + F̂ . +N = - \\widehat{𝖩(ψ, q + η)} - \\widehat{U ∂_x (q + η)} + \\widehat{(∂_x ψ)(∂_y² U)} + F̂ . ``` """ function calcN!(N, sol, t, clock, vars, params, grid) From 1d2509fbbf6ee10be998f7c73290a7c8c29e983c Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 29 Jul 2024 13:56:47 +1000 Subject: [PATCH 36/36] drop GeophysicalFlows from docs/Project.toml --- docs/Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index dfbe88eb..cf889bcc 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -3,7 +3,6 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" -GeophysicalFlows = "44ee3b1c-bc02-53fa-8355-8e347616e15e" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"