diff --git a/src/multilayerqg.jl b/src/multilayerqg.jl index 940777a5..8475d8b8 100644 --- a/src/multilayerqg.jl +++ b/src/multilayerqg.jl @@ -133,9 +133,9 @@ function Problem(nlayers::Int, # number of fluid layers end """ - Params{T, Aphys3D, Aphys2D, Aphys1D, Atrans4D, Trfft}(nlayers, g, f₀, β, ρ, H, U, eta, μ, ν, nν, calcFq!, g′, Qx, Qy, S, S⁻¹, rfftplan) + struct Params{T, Aphys3D, Aphys2D, Aphys1D, Atrans4D, Trfft} <: AbstractParams -A struct containing the parameters for the MultiLayerQG problem. Included are: +The parameters for the MultiLayerQG problem. $(TYPEDFIELDS) """ @@ -182,9 +182,9 @@ struct Params{T, Aphys3D, Aphys2D, Aphys1D, Atrans4D, Trfft} <: AbstractParams end """ - SingleLayerParams{T, Aphys3D, Aphys2D, Trfft}(β, U, eta, μ, ν, nν, calcFq!, Qx, Qy, rfftplan) + struct SingleLayerParams{T, Aphys3D, Aphys2D, Trfft} <: AbstractParams -A struct containing the parameters for the SingleLayerQG problem. Included are: +The parameters for the SingleLayerQG problem. $(TYPEDFIELDS) """ @@ -215,9 +215,9 @@ struct SingleLayerParams{T, Aphys3D, Aphys2D, Trfft} <: AbstractParams end """ - TwoLayerParams{T, Aphys3D, Aphys2D, Trfft}(g, f₀, β, ρ, H, U, eta, μ, ν, nν, calcFq!, g′, Qx, Qy, rfftplan) + TwoLayerParams{T, Aphys3D, Aphys2D, Trfft} <: AbstractParams -A struct containing the parameters for the TwoLayerQG problem. Included are: +The parameters for the TwoLayerQG problem. $(TYPEDFIELDS) """ @@ -378,7 +378,7 @@ end """ LinearEquation(dev, params, grid) -Return the `equation` for a multi-layer quasi-geostrophic problem with `params` and `grid`. +Return the equation for a multi-layer quasi-geostrophic problem with `params` and `grid`. The linear opeartor ``L`` includes only (hyper)-viscosity and is computed via `hyperviscosity(dev, params, grid)`. @@ -393,7 +393,7 @@ end """ Equation(dev, params, grid) -Return the `equation` for a multi-layer quasi-geostrophic problem with `params` and `grid`. +Return the equation for a multi-layer quasi-geostrophic problem with `params` and `grid`. The linear opeartor ``L`` includes only (hyper)-viscosity and is computed via `hyperviscosity(dev, params, grid)`. @@ -411,9 +411,9 @@ end # ---- """ - Vars{Aphys, Atrans, F, P}(q, ψ, u, v, qh, , ψh, uh, vh, Fh, prevsol) + struct Vars{Aphys, Atrans, F, P} <: AbstractVars -The variables for MultiLayer QG: +The variables for MultiLayer QG. $(FIELDS) """ @@ -447,7 +447,7 @@ const StochasticForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractAr """ DecayingVars(dev, grid, params) -Return the vars for unforced multi-layer QG problem with `grid` and `params`. +Return the variables for an unforced multi-layer QG problem with `grid` and `params`. """ function DecayingVars(dev::Dev, grid, params) where Dev T = eltype(grid) @@ -462,7 +462,7 @@ end """ ForcedVars(dev, grid, params) -Return the vars for forced multi-layer QG problem with `grid` and `params`. +Return the variables for a forced multi-layer QG problem with `grid` and `params`. """ function ForcedVars(dev::Dev, grid, params) where Dev T = eltype(grid) @@ -477,7 +477,7 @@ end """ StochasticForcedVars(dev, rid, params) -Return the vars for forced multi-layer QG problem with `grid` and `params`. +Return the variables for a forced multi-layer QG problem with `grid` and `params`. """ function StochasticForcedVars(dev::Dev, grid, params) where Dev T = eltype(grid) @@ -620,7 +620,7 @@ end """ calcS!(S, Fp, Fm, nlayers, grid) -Constructs the array ``𝕊``, which consists of `nlayer` x `nlayer` static arrays ``𝕊_𝐤`` that +Construct the array ``𝕊``, which consists of `nlayer` x `nlayer` static arrays ``𝕊_𝐤`` that relate the ``q̂_j``'s and ``ψ̂_j``'s for every wavenumber: ``q̂_𝐤 = 𝕊_𝐤 ψ̂_𝐤``. """ function calcS!(S, Fp, Fm, nlayers, grid) @@ -638,7 +638,7 @@ end """ calcS⁻¹!(S, Fp, Fm, nlayers, grid) -Constructs the array ``𝕊⁻¹``, which consists of `nlayer` x `nlayer` static arrays ``(𝕊_𝐤)⁻¹`` +Construct the array ``𝕊⁻¹``, which consists of `nlayer` x `nlayer` static arrays ``(𝕊_𝐤)⁻¹`` that relate the ``q̂_j``'s and ``ψ̂_j``'s for every wavenumber: ``ψ̂_𝐤 = (𝕊_𝐤)⁻¹ q̂_𝐤``. """ function calcS⁻¹!(S⁻¹, Fp, Fm, nlayers, grid) @@ -665,6 +665,7 @@ end calcN!(N, sol, t, clock, vars, params, grid) Compute the nonlinear term, that is the advection term, the bottom drag, and the forcing: + ```math N_j = - \\widehat{𝖩(ψ_j, q_j)} - \\widehat{U_j ∂_x Q_j} - \\widehat{U_j ∂_x q_j} + \\widehat{(∂_y ψ_j)(∂_x Q_j)} - \\widehat{(∂_x ψ_j)(∂_y Q_j)} + δ_{j, n} μ |𝐤|^2 ψ̂_n + F̂_j . @@ -688,6 +689,7 @@ end calcNlinear!(N, sol, t, clock, vars, params, grid) Compute the nonlinear term of the linearized equations: + ```math N_j = - \\widehat{U_j ∂_x Q_j} - \\widehat{U_j ∂_x q_j} + \\widehat{(∂_y ψ_j)(∂_x Q_j)} - \\widehat{(∂_x ψ_j)(∂_y Q_j)} + δ_{j, n} μ |𝐤|^2 ψ̂_n + F̂_j . @@ -707,6 +709,7 @@ end calcN_advection!(N, sol, vars, params, grid) Compute the advection term and stores it in `N`: + ```math N_j = - \\widehat{𝖩(ψ_j, q_j)} - \\widehat{U_j ∂_x Q_j} - \\widehat{U_j ∂_x q_j} + \\widehat{(∂_y ψ_j)(∂_x Q_j)} - \\widehat{(∂_x ψ_j)(∂_y Q_j)} . @@ -755,6 +758,7 @@ end calcN_linearadvection!(N, sol, vars, params, grid) Compute the advection term of the linearized equations and stores it in `N`: + ```math N_j = - \\widehat{U_j ∂_x Q_j} - \\widehat{U_j ∂_x q_j} + \\widehat{(∂_y ψ_j)(∂_x Q_j)} - \\widehat{(∂_x ψ_j)(∂_y Q_j)} . @@ -909,12 +913,15 @@ Return the kinetic energy of each fluid layer KE``_1, ...,`` KE``_{n}``, and the potential energy of each fluid interface PE``_{3/2}, ...,`` PE``_{n-1/2}``, where ``n`` is the number of layers in the fluid. (When ``n=1``, only the kinetic energy is returned.) -The kinetic energy at the ``j``-th fluid layer is +The kinetic energy at the ``j``-th fluid layer is + ```math 𝖪𝖤_j = \\frac{H_j}{H} \\int \\frac1{2} |{\\bf ∇} ψ_j|^2 \\frac{𝖽x 𝖽y}{L_x L_y} = \\frac1{2} \\frac{H_j}{H} \\sum_{𝐤} |𝐤|² |ψ̂_j|², \\ j = 1, ..., n , ``` + while the potential energy that corresponds to the interface ``j+1/2`` (i.e., the interface between the ``j``-th and ``(j+1)``-th fluid layer) is + ```math 𝖯𝖤_{j+1/2} = \\int \\frac1{2} \\frac{f₀^2}{g'_{j+1/2} H} (ψ_j - ψ_{j+1})^2 \\frac{𝖽x 𝖽y}{L_x L_y} = \\frac1{2} \\frac{f₀^2}{g'_{j+1/2} H} \\sum_{𝐤} |ψ̂_j - ψ̂_{j+1}|², \\ j = 1, ..., n-1 . ``` @@ -983,12 +990,15 @@ verticalfluxes``_{3/2},...,``verticalfluxes``_{n-1/2}``, where ``n`` is the tota (When ``n=1``, only the lateral fluxes are returned.) The lateral eddy fluxes within the ``j``-th fluid layer are + ```math \\textrm{lateralfluxes}_j = \\frac{H_j}{H} \\int U_j v_j ∂_y u_j \\frac{𝖽x 𝖽y}{L_x L_y} , \\ j = 1, ..., n , ``` + while the vertical eddy fluxes at the ``j+1/2``-th fluid interface (i.e., interface between the ``j``-th and ``(j+1)``-th fluid layer) are + ```math \\textrm{verticalfluxes}_{j+1/2} = \\int \\frac{f₀²}{g'_{j+1/2} H} (U_j - U_{j+1}) \\, v_{j+1} ψ_{j} \\frac{𝖽x 𝖽y}{L_x L_y} , \\ j = 1, ..., n-1. diff --git a/src/singlelayerqg.jl b/src/singlelayerqg.jl index 6c09aa84..d7fafa4a 100644 --- a/src/singlelayerqg.jl +++ b/src/singlelayerqg.jl @@ -116,9 +116,9 @@ end abstract type SingleLayerQGParams <: AbstractParams end """ - Params{T, Aphys, Atrans, ℓ}(β, deformation_radius, eta, etah, μ, ν, nν, calcF!) + struct Params{T, Aphys, Atrans, ℓ} <: SingleLayerQGParams -A struct containing the parameters for the SingleLayerQG problem. Included are: +The parameters for the `SingleLayerQG` problem. $(TYPEDFIELDS) """ @@ -145,9 +145,9 @@ const BarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractAr const EquivalentBarotropicQGParams = Params{<:AbstractFloat, <:AbstractArray, <:AbstractArray, <:AbstractFloat} """ - EquivalentBarotropicQGParams(grid::TwoDGrid, β, deformation_radius, eta, μ, ν, nν::Int, calcF + EquivalentBarotropicQGParams(grid, β, deformation_radius, eta, μ, ν, nν, calcF) -Constructor for EquivalentBarotropicQGParams (finite Rossby radius of deformation). +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_grid = typeof(eta) <: AbstractArray ? eta : FourierFlows.on_grid(eta, grid) @@ -156,9 +156,9 @@ function EquivalentBarotropicQGParams(grid::AbstractGrid{T, A}, β, deformation_ end """ - BarotropicQGParams(grid::TwoDGrid, β, eta, μ, ν, nν::Int, calcF + BarotropicQGParams(grid, β, eta, μ, ν, nν, calcF) -Constructor for BarotropicQGParams (infinite Rossby radius of deformation). +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) @@ -171,14 +171,15 @@ 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 +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: ```math L = - μ - ν |𝐤|^{2 n_ν} + i β k_x / |𝐤|² . ``` -Nonlinear term is computed via `calcN!` function. + +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 @@ -190,14 +191,15 @@ end """ Equation(params::EquivalentBarotropicQGParams, grid) -Return the `equation` for an equivalent-barotropic QG problem with `params` and `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: ```math L = -μ - ν |𝐤|^{2 n_ν} + i β k_x / (|𝐤|² + 1/ℓ²) . ``` -Nonlinear term is computed via `calcN!` function. + +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) @@ -214,7 +216,7 @@ end abstract type SingleLayerQGVars <: AbstractVars end """ - Vars{Aphys, Atrans, F, P}(q, ψ, u, v, qh, , ψh, uh, vh, Fh, prevsol) + struct Vars{Aphys, Atrans, F, P} <: SingleLayerQGVars The variables for SingleLayer QG: @@ -253,10 +255,10 @@ Return the `vars` for unforced single-layer QG problem on device `dev` and with """ function DecayingVars(dev::Dev, grid::AbstractGrid) where Dev T = eltype(grid) - + @devzeros Dev T (grid.nx, grid.ny) q u v ψ @devzeros Dev Complex{T} (grid.nkr, grid.nl) qh uh vh ψh - + Vars(q, ψ, u, v, qh, ψh, uh, vh, nothing, nothing) end @@ -267,10 +269,10 @@ Return the `vars` for forced single-layer QG problem on device dev and with `gri """ function ForcedVars(dev::Dev, grid::AbstractGrid) where Dev T = eltype(grid) - + @devzeros Dev T (grid.nx, grid.ny) q u v ψ @devzeros Dev Complex{T} (grid.nkr, grid.nl) qh uh vh ψh Fh - + return Vars(q, ψ, u, v, qh, ψh, uh, vh, Fh, nothing) end @@ -283,7 +285,7 @@ function StochasticForcedVars(dev::Dev, grid::AbstractGrid) where Dev T = eltype(grid) @devzeros Dev T (grid.nx, grid.ny) q u v ψ @devzeros Dev Complex{T} (grid.nkr, grid.nl) qh uh vh ψh Fh prevsol - + return Vars(q, ψ, u, v, qh, ψh, uh, vh, Fh, prevsol) end @@ -313,16 +315,17 @@ function calcN_advection!(N, sol, t, clock, vars, params, grid) 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+η) + @. 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+η) + @. 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+η)} + 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+η)} + mul!(vq_plus_ηh, grid.rfftplan, vq_plus_η) # \hat{v * (q + η)} + + @. 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*(q+η)]/∂x -∂[v*(q+η)]/∂y return nothing end @@ -430,14 +433,16 @@ end """ kinetic_energy(prob) - kinetic_energy(sol, grid, vars, params) -Return the domain-averaged kinetic energy of the fluid. Since ``u² + v² = |{\\bf ∇} ψ|²``, the -domain-averaged kinetic energy is +Return the problem's (`prob`) domain-averaged kinetic energy of the fluid. Since +``u² + v² = |{\\bf ∇} ψ|²``, the domain-averaged kinetic energy is + ```math \\int \\frac1{2} |{\\bf ∇} ψ|² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} |𝐤|² |ψ̂|² . ``` """ +@inline kinetic_energy(prob) = kinetic_energy(prob.sol, prob.vars, prob.params, prob.grid) + function kinetic_energy(sol, vars, params, grid) streamfunctionfrompv!(vars.ψh, sol, params, grid) @. vars.uh = sqrt.(grid.Krsq) * vars.ψh # vars.uh is a dummy variable @@ -445,33 +450,31 @@ function kinetic_energy(sol, vars, params, grid) return parsevalsum2(vars.uh , grid) / (2 * grid.Lx * grid.Ly) end -kinetic_energy(prob) = kinetic_energy(prob.sol, prob.vars, prob.params, prob.grid) - """ potential_energy(prob) - potential_energy(sol, grid, vars, params) -Return the domain-averaged potential energy of the fluid, +Return the problem's (`prob`) domain-averaged potential energy of the fluid, + ```math \\int \\frac1{2} \\frac{ψ²}{ℓ²} \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} \\frac{|ψ̂|²}{ℓ²} . ``` """ +@inline potential_energy(prob) = potential_energy(prob.sol, prob.vars, prob.params, prob.grid) + +@inline potential_energy(sol, vars, params::BarotropicQGParams, grid) = 0 + function potential_energy(sol, vars, params::EquivalentBarotropicQGParams, grid) streamfunctionfrompv!(vars.ψh, sol, params, grid) + return 1 / params.deformation_radius^2 * parsevalsum2(vars.ψh, grid) / (2 * grid.Lx * grid.Ly) end -@inline potential_energy(sol, vars, params::BarotropicQGParams, grid) = 0 - -@inline potential_energy(prob) = potential_energy(prob.sol, prob.vars, prob.params, prob.grid) - """ energy(prob) - energy(sol, grid, vars, params) -Return the domain-averaged total energy of the fluid, that is, the kinetic energy for a -pure barotropic flow or the sum of kinetic and potential energies for an equivalent barotropic -flow. +Return the problem's (`prob`) domain-averaged total energy of the fluid, that is, the kinetic +energy for a pure barotropic flow *or* the sum of kinetic and potential energies for an equivalent +barotropic flow. """ @inline energy(prob) = energy(prob.sol, prob.vars, prob.params, prob.grid) @@ -481,9 +484,9 @@ flow. """ enstrophy(prob) - enstrophy(sol, vars, params, grid) -Return the domain-averaged enstrophy +Return the problem's (`prob`) domain-averaged enstrophy + ```math \\int \\frac1{2} (q + η)² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} |q̂ + η̂|² . ``` @@ -497,10 +500,11 @@ end """ energy_dissipation(prob) - energy_dissipation(sol, vars, params, grid) -Return the domain-averaged energy dissipation rate. nν must be >= 1. +Return the problem's (`prob`) domain-averaged energy dissipation rate. """ +@inline energy_dissipation(prob) = energy_dissipation(prob.sol, prob.vars, prob.params, prob.grid) + @inline function energy_dissipation(sol, vars, params::BarotropicQGParams, grid) energy_dissipationh = vars.uh # use vars.uh as scratch variable @@ -511,14 +515,13 @@ end energy_dissipation(sol, vars, params::EquivalentBarotropicQGParams, grid) = error("not implemented for finite deformation radius") -@inline energy_dissipation(prob) = energy_dissipation(prob.sol, prob.vars, prob.params, prob.grid) - """ enstrophy_dissipation(prob) - enstrophy_dissipation(sol, vars, params, grid) -Return the domain-averaged enstrophy dissipation rate. nν must be >= 1. +Return the problem's (`prob`) domain-averaged enstrophy dissipation rate. """ +@inline enstrophy_dissipation(prob) = enstrophy_dissipation(prob.sol, prob.vars, prob.params, prob.grid) + @inline function enstrophy_dissipation(sol, vars, params::BarotropicQGParams, grid) enstrophy_dissipationh = vars.uh # use vars.uh as scratch variable @@ -529,14 +532,13 @@ end @inline enstrophy_dissipation(sol, vars, params::EquivalentBarotropicQGParams, grid) = error("not implemented for finite deformation radius") -@inline enstrophy_dissipation(prob) = enstrophy_dissipation(prob.sol, prob.vars, prob.params, prob.grid) - """ energy_work(prob) - energy_work(sol, vars, params, grid) -Return the domain-averaged rate of work of energy by the forcing `Fh`. +Return the problem's (`prob`) domain-averaged rate of work of energy by the forcing `F`. """ +@inline energy_work(prob) = energy_work(prob.sol, prob.vars, prob.params, prob.grid) + @inline function energy_work(sol, vars::ForcedVars, params::BarotropicQGParams, grid) energy_workh = vars.uh # use vars.uh as scratch variable @@ -553,14 +555,13 @@ end @inline energy_work(sol, vars, params::EquivalentBarotropicQGParams, grid) = error("not implemented for finite deformation radius") -@inline energy_work(prob) = energy_work(prob.sol, prob.vars, prob.params, prob.grid) - """ enstrophy_work(prob) - enstrophy_work(sol, vars, params, grid) -Return the domain-averaged rate of work of enstrophy by the forcing `Fh`. +Return the problem's (`prob`) domain-averaged rate of work of enstrophy by the forcing ``F``. """ +@inline enstrophy_work(prob) = enstrophy_work(prob.sol, prob.vars, prob.params, prob.grid) + @inline function enstrophy_work(sol, vars::ForcedVars, params::BarotropicQGParams, grid) enstrophy_workh = vars.uh # use vars.uh as scratch variable @@ -577,14 +578,13 @@ end @inline enstrophy_work(sol, vars, params::EquivalentBarotropicQGParams, grid) = error("not implemented for finite deformation radius") -@inline enstrophy_work(prob) = enstrophy_work(prob.sol, prob.vars, prob.params, prob.grid) - """ energy_drag(prob) - energy_drag(sol, vars, params, grid) -Return the extraction of domain-averaged energy by drag μ. +Return the problem's (`prob`) extraction of domain-averaged energy by drag ``μ``. """ +@inline energy_drag(prob) = energy_drag(prob.sol, prob.vars, prob.params, prob.grid) + @inline function energy_drag(sol, vars, params::BarotropicQGParams, grid) energy_dragh = vars.uh # use vars.uh as scratch variable @@ -595,14 +595,13 @@ end @inline energy_drag(sol, vars, params::EquivalentBarotropicQGParams, grid) = error("not implemented for finite deformation radius") -@inline energy_drag(prob) = energy_drag(prob.sol, prob.vars, prob.params, prob.grid) - """ enstrophy_drag(prob) - enstrophy_drag(sol, vars, params, grid) -Return the extraction of domain-averaged enstrophy by drag/hypodrag μ. +Return the problem's (`prob`) extraction of domain-averaged enstrophy by drag/hypodrag ``μ``. """ +@inline enstrophy_drag(prob) = enstrophy_drag(prob.sol, prob.vars, prob.params, prob.grid) + @inline function enstrophy_drag(sol, vars, params::BarotropicQGParams, grid) enstrophy_dragh = vars.uh # use vars.uh as scratch variable @@ -613,6 +612,4 @@ end @inline enstrophy_drag(sol, vars, params::EquivalentBarotropicQGParams, grid) = error("not implemented for finite deformation radius") -@inline enstrophy_drag(prob) = enstrophy_drag(prob.sol, prob.vars, prob.params, prob.grid) - end # module diff --git a/src/twodnavierstokes.jl b/src/twodnavierstokes.jl index be691a7a..ff80ed92 100644 --- a/src/twodnavierstokes.jl +++ b/src/twodnavierstokes.jl @@ -83,9 +83,9 @@ function Problem(dev::Device=CPU(); aliased_fraction = 1/3, T = Float64) - grid = TwoDGrid(dev, nx, Lx, ny, Ly; aliased_fraction=aliased_fraction, T=T) + grid = TwoDGrid(dev, nx, Lx, ny, Ly; aliased_fraction=aliased_fraction, T) - params = Params{T}(ν, nν, μ, nμ, calcF) + params = Params(T(ν), nν, T(μ), nμ, calcF) vars = calcF == nothingfunction ? DecayingVars(dev, grid) : (stochastic ? StochasticForcedVars(dev, grid) : ForcedVars(dev, grid)) @@ -100,9 +100,9 @@ end # ---------- """ - Params{T}(ν, nν, μ, nμ, calcF!) + struct Params{T} <: AbstractParams -A struct containing the parameters for the two-dimensional Navier-Stokes. Included are: +The parameters for a two-dimensional Navier-Stokes problem: $(TYPEDFIELDS) """ @@ -137,7 +137,7 @@ hypo-viscocity of order ``n_μ`` with coefficient ``μ``, L = - ν |𝐤|^{2 n_ν} - μ |𝐤|^{2 n_μ} . ``` -Plain old viscocity corresponds to ``n_ν=1`` while ``n_μ=0`` corresponds to linear drag. +Plain-old viscocity corresponds to ``n_ν = 1`` while ``n_μ = 0`` corresponds to linear drag. The nonlinear term is computed via the function `calcN!`. """ @@ -156,9 +156,9 @@ end abstract type TwoDNavierStokesVars <: AbstractVars end """ - Vars{Aphys, Atrans, F, P}(ζ, u, v, ζh, uh, vh, Fh, prevsol) + struct Vars{Aphys, Atrans, F, P} <: TwoDNavierStokesVars -The variables for two-dimensional Navier-Stokes: +The variables for two-dimensional Navier-Stokes problem: $(FIELDS) """ @@ -188,7 +188,7 @@ const StochasticForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractAr """ DecayingVars(dev, grid) -Return the `vars` for unforced two-dimensional Navier-Stokes problem on device `dev` and +Return the variables `vars` for unforced two-dimensional Navier-Stokes problem on device `dev` and with `grid`. """ function DecayingVars(::Dev, grid::AbstractGrid) where Dev @@ -203,7 +203,7 @@ end """ ForcedVars(dev, grid) -Return the `vars` for forced two-dimensional Navier-Stokes on device `dev` and with `grid`. +Return the variables `vars` for forced two-dimensional Navier-Stokes on device `dev` and with `grid`. """ function ForcedVars(dev::Dev, grid::AbstractGrid) where Dev T = eltype(grid) @@ -217,7 +217,7 @@ end """ StochasticForcedVars(dev, grid) -Return the `vars` for stochastically forced two-dimensional Navier-Stokes on device `dev` and +Return the variables `vars` for stochastically forced two-dimensional Navier-Stokes on device `dev` and with `grid`. """ function StochasticForcedVars(dev::Dev, grid::AbstractGrid) where Dev @@ -290,7 +290,7 @@ end """ addforcing!(N, sol, t, clock, vars, params, grid) -When the problem includes forcing, calculate the forcing term ``F̂`` and add it to the +When the problem includes forcing, calculate the forcing term ``F̂`` and add it to the nonlinear term ``N``. """ addforcing!(N, sol, t, clock, vars::DecayingVars, params, grid) = nothing @@ -320,7 +320,7 @@ end """ updatevars!(prob) -Update variables in `vars` with solution in `sol`. +Update problem's variables in `prob.vars` using the state in `prob.sol`. """ function updatevars!(prob) vars, grid, sol = prob.vars, prob.grid, prob.sol @@ -341,7 +341,7 @@ end """ set_ζ!(prob, ζ) -Set the solution `sol` as the transform of `ζ` and then update variables in `vars`. +Set the solution `sol` as the transform of `ζ` and then update variables in `prob.vars`. """ function set_ζ!(prob, ζ) mul!(prob.sol, prob.grid.rfftplan, ζ) @@ -360,8 +360,10 @@ Return the domain-averaged kinetic energy. Since ``u² + v² = |{\\bf ∇} ψ|² kinetic energy is ```math -\\int \\frac1{2} |{\\bf ∇} ψ|² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} |𝐤|² |ψ̂|² . +\\int \\frac1{2} |{\\bf ∇} ψ|² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} |𝐤|² |ψ̂|² , ``` + +where ``ψ`` is the streamfunction. """ @inline function energy(prob) sol, vars, grid = prob.sol, prob.vars, prob.grid @@ -374,16 +376,15 @@ end """ enstrophy(prob) -Returns the domain-averaged enstrophy, +Returns the problem's (`prob`) domain-averaged enstrophy, ```math -\\int \\frac1{2} ζ² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} |ζ̂|² . +\\int \\frac1{2} ζ² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} |ζ̂|² , ``` + +where ``ζ`` is the relative vorticity. """ -@inline function enstrophy(prob) - sol, grid = prob.sol, prob.grid - return 1 / (2 * grid.Lx * grid.Ly) * parsevalsum(abs2.(sol), grid) -end +@inline enstrophy(prob) = 1 / (2 * prob.grid.Lx * prob.grid.Ly) * parsevalsum(abs2.(prob.sol), prob.grid) """ energy_dissipation(prob, ξ, νξ) @@ -406,23 +407,23 @@ where ``ξ`` and ``nξ`` could be either the (hyper)-viscosity coefficient ``ν` end """ - energy_dissipation_hyperviscosity(prob, ξ, νξ) + energy_dissipation_hyperviscosity(prob) -Return the domain-averaged energy dissipation rate done by the ``ν`` (hyper)-viscosity. +Return the problem's (`prob`) domain-averaged energy dissipation rate done by the ``ν`` (hyper)-viscosity. """ energy_dissipation_hyperviscosity(prob) = energy_dissipation(prob, prob.params.ν, prob.params.nν) """ - energy_dissipation_hypoviscosity(prob, ξ, νξ) + energy_dissipation_hypoviscosity(prob) -Return the domain-averaged energy dissipation rate done by the ``μ`` (hypo)-viscosity. +Return the problem's (`prob`) domain-averaged energy dissipation rate done by the ``μ`` (hypo)-viscosity. """ energy_dissipation_hypoviscosity(prob) = energy_dissipation(prob, prob.params.μ, prob.params.nμ) """ enstrophy_dissipation(prob, ξ, νξ) -Return the domain-averaged enstrophy dissipation rate done by the viscous term, +Return the problem's (`prob`) domain-averaged enstrophy dissipation rate done by the viscous term, ```math ξ (-1)^{n_ξ+1} \\int ζ ∇^{2n_ξ} ζ \\frac{𝖽x 𝖽y}{L_x L_y} = - ξ \\sum_{𝐤} |𝐤|^{2n_ξ} |ζ̂|² , @@ -441,29 +442,32 @@ where ``ξ`` and ``nξ`` could be either the (hyper)-viscosity coefficient ``ν` end """ - enstrophy_dissipation_hyperviscosity(prob, ξ, νξ) + enstrophy_dissipation_hyperviscosity(prob) -Return the domain-averaged enstrophy dissipation rate done by the ``ν`` (hyper)-viscosity. +Return the problem's (`prob`) domain-averaged enstrophy dissipation rate done by the ``ν`` (hyper)-viscosity. """ enstrophy_dissipation_hyperviscosity(prob) = enstrophy_dissipation(prob, prob.params.ν, prob.params.nν) """ - enstrophy_dissipation_hypoviscosity(prob, ξ, νξ) + enstrophy_dissipation_hypoviscosity(prob) -Return the domain-averaged enstrophy dissipation rate done by the ``μ`` (hypo)-viscosity. +Return the problem's (`prob`) domain-averaged enstrophy dissipation rate done by the ``μ`` (hypo)-viscosity. """ enstrophy_dissipation_hypoviscosity(prob) = enstrophy_dissipation(prob, prob.params.μ, prob.params.nμ) """ energy_work(prob) - energy_work(sol, vars, grid) -Return the domain-averaged rate of work of energy by the forcing ``F``, +Return the problem's (`prob`) domain-averaged rate of work of energy by the forcing ``F``, ```math - \\int ψ F \\frac{𝖽x 𝖽y}{L_x L_y} = - \\sum_{𝐤} ψ̂ F̂^* . ``` + +where ``ψ`` is the stream flow. """ +@inline energy_work(prob) = energy_work(prob.sol, prob.vars, prob.grid) + @inline function energy_work(sol, vars::ForcedVars, grid) energy_workh = vars.uh # use vars.uh as scratch variable @@ -478,18 +482,19 @@ end return 1 / (grid.Lx * grid.Ly) * parsevalsum(energy_workh, grid) end -@inline energy_work(prob) = energy_work(prob.sol, prob.vars, prob.grid) - """ enstrophy_work(prob) - enstrophy_work(sol, vars, grid) -Return the domain-averaged rate of work of enstrophy by the forcing ``F``, +Return the problem's (`prob`) domain-averaged rate of work of enstrophy by the forcing ``F``, ```math -\\int ζ F \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} ζ̂ F̂^* . +\\int ζ F \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} ζ̂ F̂^* , ``` + +where ``ζ`` is the relative vorticity. """ +@inline enstrophy_work(prob) = enstrophy_work(prob.sol, prob.vars, prob.grid) + @inline function enstrophy_work(sol, vars::ForcedVars, grid) enstrophy_workh = vars.uh # use vars.uh as scratch variable @@ -499,11 +504,9 @@ end @inline function enstrophy_work(sol, vars::StochasticForcedVars, grid) enstrophy_workh = vars.uh # use vars.uh as scratch variable - + @. enstrophy_workh = (vars.prevsol + sol) / 2 * conj(vars.Fh) return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_workh, grid) end -@inline enstrophy_work(prob) = enstrophy_work(prob.sol, prob.vars, prob.grid) - end # module diff --git a/test/test_twodnavierstokes.jl b/test/test_twodnavierstokes.jl index 09a287ab..12a2b07e 100644 --- a/test/test_twodnavierstokes.jl +++ b/test/test_twodnavierstokes.jl @@ -1,6 +1,6 @@ function test_twodnavierstokes_lambdipole(n, dt, dev::Device=CPU(); L=2π, Ue=1, Re=L/20, ν=0.0, nν=1, ti=L/Ue*0.01, nm=3) nt = round(Int, ti/dt) - prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, dt=dt, stepper="FilteredRK4") + prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν, nν, dt, stepper="FilteredRK4") x, y = gridpoints(prob.grid) ζ = prob.vars.ζ