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" 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 diff --git a/docs/src/modules/singlelayerqg.md b/docs/src/modules/singlelayerqg.md index 6f6e382c..e8e09603 100644 --- a/docs/src/modules/singlelayerqg.md +++ b/docs/src/modules/singlelayerqg.md @@ -20,33 +20,50 @@ 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``. +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 + \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 +139,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. diff --git a/src/multilayerqg.jl b/src/multilayerqg.jl index 4e35d679..99668a8b 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. @@ -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} @@ -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 diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 2ce0fc87..baaf26e0 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -1,5 +1,7 @@ module SingleLayerQG +using LinearAlgebra + export Problem, streamfunctionfrompv!, @@ -37,6 +39,7 @@ nothingfunction(args...) = nothing Ly = Lx, β = 0.0, deformation_radius = Inf, + U = 0.0, eta = nothing, ν = 0.0, nν = 1, @@ -62,6 +65,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`: Imposed background constant zonal flow ``U(y)``. - `eta`: Topographic potential vorticity. - `ν`: Small-scale (hyper)-viscosity coefficient. - `nν`: (Hyper)-viscosity order, `nν```≥ 1``. @@ -82,6 +86,7 @@ function Problem(dev::Device=CPU(); # Physical parameters β = 0.0, deformation_radius = Inf, + U = 0.0, eta = nothing, # Drag and (hyper-)viscosity ν = 0.0, @@ -98,12 +103,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))) + + U = U isa Number ? convert(T, U) : U - 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 || + 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)) @@ -120,17 +128,19 @@ 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, ℓ} <: SingleLayerQGParams +struct Params{T, Aphys, Atrans, Tℓ, TU} <: 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" eta :: Aphys "Fourier transform of topographic potential vorticity" @@ -143,31 +153,59 @@ struct Params{T, Aphys, Atrans, ℓ} <: SingleLayerQGParams 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 :: Aphys + "array containing ``y``-gradient of PV due to ``U`` and topographic PV" + 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, <:Number} +const SingleLayerQGvaryingUParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:Any, <:AbstractArray} """ - 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} - eta_on_grid = typeof(eta) <: AbstractArray ? A(eta) : FourierFlows.on_grid(eta, grid) - etah_on_grid = rfft(eta_on_grid) +function EquivalentBarotropicQGParams(grid::AbstractGrid{T, A}, deformation_radius, β, U, eta, μ, ν, nν::Int, calcF) where {T, A} + + 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 = eta isa AbstractArray ? A(eta) : FourierFlows.on_grid(eta, grid) + etah = rfft(eta_on_grid) + + Qx = irfft(im * grid.kr .* etah, grid.nx) # ∂η/∂x + Qy = irfft(im * grid.l .* etah, grid.nx) # ∂η/∂y + + if U isa AbstractArray + Uh = rfft(U) + Uyy = irfft(- grid.l.^2 .* Uh, grid.nx) # ∂²U/∂y² + Qy .-= Uyy # -∂²U/∂y² + end + + # Note: The β-term in Qy is included in the linear term L of Equation. - return Params(β, deformation_radius, eta_on_grid, etah_on_grid, μ, ν, nν, calcF) + return Params(β, deformation_radius, U, eta_on_grid, etah, μ, ν, nν, calcF, Qx, Qy) 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) - +function BarotropicQGParams(grid, β, U, eta, μ, ν, nν::Int, calcF) + deformation_radius = nothing + + return EquivalentBarotropicQGParams(grid, deformation_radius, β, U, eta, μ, ν, nν, calcF) +end + # --------- # Equations @@ -175,41 +213,41 @@ BarotropicQGParams(grid::AbstractGrid, β, 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 ``ν`` -and the ``β`` term: +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 / |𝐤|² . +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 + @. L -= im * params.U * grid.kr + end + CUDA.@allowscalar L[1, 1] = 0 - + 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 ``ν`` and the ``β`` term: +function Equation(params::EquivalentBarotropicQGParams, grid) + L = @. - params.μ - params.ν * grid.Krsq^params.nν + im * params.β * grid.kr / (grid.Krsq + 1 / params.deformation_radius^2) -```math -L = -μ - ν |𝐤|^{2 n_ν} + i β k_x / (|𝐤|² + 1/ℓ²) . -``` + if params.U isa Number + @. L -= im * params.U * grid.kr + end -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) CUDA.@allowscalar L[1, 1] = 0 - + return FourierFlows.Equation(L, calcN!, grid) end @@ -304,16 +342,19 @@ 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+η)`` in conservative -form, i.e., ``- ∂_x[(∂_y ψ)(q+η)] - ∂_y[(∂_x ψ)(q+η)]`` 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 + η)} = - i k_x \\widehat{u (q + η)} - i k_y \\widehat{v (q + η)} . ``` """ -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) @. vars.uh = -im * grid.l * vars.ψh @@ -335,6 +376,58 @@ function calcN_advection!(N, sol, t, clock, vars, params, grid) @. 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 # - \hat{∂(Uη)/∂x} + + 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{(∂_x ψ)(∂_y² U)} . +``` +""" +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.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) + + 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) + + 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!(vqh, grid.rfftplan, vq) + + @. N -= im * grid.kr * uqh + im * grid.l * vqh # -\hat{∂[(U+u)q]/∂x} - \hat{∂[vq]/∂y} + return nothing end @@ -344,7 +437,7 @@ end Calculate the nonlinear term, that is the advection term and the forcing, ```math -N = - \\widehat{𝖩(ψ, q + η)} + F̂ . +N = - \\widehat{𝖩(ψ, q + η)} - \\widehat{U ∂_x (q + η)} + \\widehat{(∂_x ψ)(∂_y² U)} + F̂ . ``` """ function calcN!(N, sol, t, clock, vars, params, grid) diff --git a/test/runtests.jl b/test/runtests.jl index 6d542b55..95e193c1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -48,32 +48,37 @@ 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 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=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_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_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) @test test_1layerqg_stochasticforcing_enstrophybudget(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) @@ -98,7 +103,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 +117,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 2cca28e3..74289e03 100644 --- a/test/test_singlelayerqg.jl +++ b/test/test_singlelayerqg.jl @@ -1,10 +1,12 @@ +using LinearAlgebra + """ test_1layerqg_rossbywave(stepper, dt, nsteps, dev; kwargs...) 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 +21,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 +40,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 = @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 @@ -255,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_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(); + add_topography = false, + add_background_flow = false, + background_flow_vary_in_y = true) n, L = 128, 2π ν, nν = 1e-2, 1 μ = 0.0 @@ -275,13 +286,42 @@ function test_1layerqg_advection(dt, stepper, dev::Device=CPU(); n=128, L=2π, 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) + η_array = @. η(x, y, η₀) + ψ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)))) + + 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² + 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) @@ -290,14 +330,71 @@ function test_1layerqg_advection(dt, stepper, dev::Device=CPU(); n=128, L=2π, 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 && 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₀, eta = η_array, ν=ν, nν=nν, μ=μ, dt=dt, 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_nonlinearadvection_deformation(dt, stepper, dev::Device=CPU(); + deformation_radius=1.23) + +Same as `test_1layerqg_nonlinearadvection` but with finite deformation radius. +""" +function test_1layerqg_nonlinearadvection_deformation(dt, stepper, dev::Device=CPU(); + deformation_radius=1.23) + 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) + + η₀ = 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 + + # 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)))) + + 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 @@ -309,50 +406,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) @@ -361,21 +467,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) @@ -395,7 +504,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.ψ ≈ ψ) @@ -435,4 +543,4 @@ 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