From 39144486527486a0537e70cbe3da38030bc02a69 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Thu, 20 Jul 2023 19:30:19 -0400 Subject: [PATCH 01/18] introducing stable upwinding? --- src/SpeedyWeather.jl | 1 + src/abstract_types.jl | 3 + src/dynamics/diagnostic_variables.jl | 10 +++- src/dynamics/models.jl | 3 + src/dynamics/tendencies_dynamics.jl | 82 ++----------------------- src/dynamics/vertical_advection.jl | 90 ++++++++++++++++++++++++++++ 6 files changed, 110 insertions(+), 79 deletions(-) create mode 100644 src/dynamics/vertical_advection.jl diff --git a/src/SpeedyWeather.jl b/src/SpeedyWeather.jl index 74777bc97..14dd3dff8 100644 --- a/src/SpeedyWeather.jl +++ b/src/SpeedyWeather.jl @@ -155,6 +155,7 @@ include("dynamics/forcing.jl") include("dynamics/geopotential.jl") include("dynamics/initial_conditions.jl") include("dynamics/horizontal_diffusion.jl") +include("dynamics/vertical_advection.jl") include("dynamics/implicit.jl") include("dynamics/scaling.jl") include("dynamics/tendencies.jl") diff --git a/src/abstract_types.jl b/src/abstract_types.jl index 5f7dea0c6..fb730a348 100644 --- a/src/abstract_types.jl +++ b/src/abstract_types.jl @@ -26,6 +26,9 @@ abstract type AbstractColumnVariables{NF} end # FORCING (Barotropic and ShallowWaterModel) abstract type AbstractForcing{NF} end +# VERTICAL ADVECTION (PrimitiveEquation) +abstract type VerticalAdvection{NF} end + # PARAMETERIZATIONS abstract type AbstractParameterization{NF} end abstract type BoundaryLayerDrag{NF} <: AbstractParameterization{NF} end diff --git a/src/dynamics/diagnostic_variables.jl b/src/dynamics/diagnostic_variables.jl index 141ba6484..37b1d89d6 100644 --- a/src/dynamics/diagnostic_variables.jl +++ b/src/dynamics/diagnostic_variables.jl @@ -46,11 +46,14 @@ struct GridVariables{NF<:AbstractFloat,Grid<:AbstractGrid{NF}} vor_grid ::Grid # vorticity div_grid ::Grid # divergence temp_grid ::Grid # absolute temperature [K] + temp_grid_prev ::Grid temp_virt_grid ::Grid # virtual tempereature [K] humid_grid ::Grid # specific_humidity geopot_grid ::Grid # geopotential (is that needed?) u_grid ::Grid # zonal velocity *coslat [m/s] v_grid ::Grid # meridional velocity *coslat [m/s] + u_grid_prev ::Grid # zonal velocity *coslat [m/s] + v_grid_prev ::Grid # meridional velocity *coslat [m/s] end function Base.zeros(::Type{GridVariables},SG::SpectralGrid) @@ -59,14 +62,17 @@ function Base.zeros(::Type{GridVariables},SG::SpectralGrid) vor_grid = zeros(Grid{NF},nlat_half) # vorticity div_grid = zeros(Grid{NF},nlat_half) # divergence temp_grid = zeros(Grid{NF},nlat_half) # absolute temperature + temp_grid_prev = zeros(Grid{NF},nlat_half) # absolute temperature temp_virt_grid = zeros(Grid{NF},nlat_half) # virtual temperature humid_grid = zeros(Grid{NF},nlat_half) # specific humidity geopot_grid = zeros(Grid{NF},nlat_half) # geopotential u_grid = zeros(Grid{NF},nlat_half) # zonal velocity *coslat v_grid = zeros(Grid{NF},nlat_half) # meridonal velocity *coslat + u_grid_prev = zeros(Grid{NF},nlat_half) # zonal velocity *coslat + v_grid_prev = zeros(Grid{NF},nlat_half) # meridonal velocity *coslat - return GridVariables( vor_grid,div_grid,temp_grid,temp_virt_grid,humid_grid,geopot_grid, - u_grid,v_grid) + return GridVariables( vor_grid,div_grid,temp_grid,temp_grid_prev,temp_virt_grid,humid_grid,geopot_grid, + u_grid,v_grid,u_grid_prev,v_grid_prev) end """ diff --git a/src/dynamics/models.jl b/src/dynamics/models.jl index 63b143998..0bd1ec87e 100644 --- a/src/dynamics/models.jl +++ b/src/dynamics/models.jl @@ -139,6 +139,7 @@ Base.@kwdef struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr boundary_layer_drag::BoundaryLayerDrag{NF} = LinearDrag(spectral_grid) temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid) static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid) + vertical_advection::VerticalAdvection{NF} = FirstOrderUpwind(spectral_grid) # vertical_diffusion::VerticalDiffusion{NF} = VerticalLaplacian(spectral_grid) # NUMERICS @@ -206,6 +207,8 @@ Base.@kwdef struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid) static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid) large_scale_condensation::AbstractCondensation{NF} = SpeedyCondensation(spectral_grid) + vertical_advection::VerticalAdvection{NF} = FirstOrderUpwind(spectral_grid) + # vertical_diffusion::VerticalDiffusion{NF} = VerticalLaplacian(spectral_grid) # NUMERICS diff --git a/src/dynamics/tendencies_dynamics.jl b/src/dynamics/tendencies_dynamics.jl index 9a38c0d72..3c4e815da 100644 --- a/src/dynamics/tendencies_dynamics.jl +++ b/src/dynamics/tendencies_dynamics.jl @@ -162,83 +162,6 @@ function vertical_velocity!( (uv∇lnp_sum_above .+ Δσₖ*uv∇lnp) # and level k σₖ-weighted uv∇lnp here end -function vertical_advection!( layer::DiagnosticVariablesLayer, - diagn::DiagnosticVariables, - model::PrimitiveEquation) - - (; k ) = layer # which layer are we on? - - wet_core = model isa PrimitiveWet - (; σ_levels_thick, nlev ) = model.geometry - - # for k==1 "above" term is 0, for k==nlev "below" term is zero - # avoid out-of-bounds indexing with k_above, k_below as follows - k_above = max(1,k-1) # just saturate, because M_1/2 = 0 (which zeros that term) - k_below = min(k+1,nlev) # just saturate, because M_nlev+1/2 = 0 (which zeros that term) - - # mass fluxes, M_1/2 = M_nlev+1/2 = 0, but k=1/2 isn't explicitly stored - σ_tend_above = diagn.layers[k_above].dynamics_variables.σ_tend - σ_tend_below = layer.dynamics_variables.σ_tend - - # layer thickness Δσ on level k - Δσₖ = σ_levels_thick[k] - - u_tend = layer.tendencies.u_tend_grid # zonal velocity - u_above = diagn.layers[k_above].grid_variables.u_grid - u = layer.grid_variables.u_grid - u_below = diagn.layers[k_below].grid_variables.u_grid - - _vertical_advection!(u_tend,σ_tend_above,σ_tend_below,u_above,u,u_below,Δσₖ) - - v_tend = layer.tendencies.v_tend_grid # meridional velocity - v_above = diagn.layers[k_above].grid_variables.v_grid - v = layer.grid_variables.v_grid - v_below = diagn.layers[k_below].grid_variables.v_grid - - _vertical_advection!(v_tend,σ_tend_above,σ_tend_below,v_above,v,v_below,Δσₖ) - - T_tend = layer.tendencies.temp_tend_grid # temperature - T_above = diagn.layers[k_above].grid_variables.temp_grid - T = layer.grid_variables.temp_grid - T_below = diagn.layers[k_below].grid_variables.temp_grid - - _vertical_advection!(T_tend,σ_tend_above,σ_tend_below,T_above,T,T_below,Δσₖ) - - if wet_core - q_tend = layer.tendencies.humid_tend_grid # humidity - q_above = diagn.layers[k_above].grid_variables.humid_grid - q = layer.grid_variables.humid_grid - q_below = diagn.layers[k_below].grid_variables.humid_grid - - _vertical_advection!(q_tend,σ_tend_above,σ_tend_below,q_above,q,q_below,Δσₖ) - end -end - -# MULTI THREADED VERSION only writes into layer k -function _vertical_advection!( ξ_tend::Grid, # tendency of quantity ξ at k - σ_tend_above::Grid, # vertical velocity at k-1/2 - σ_tend_below::Grid, # vertical velocity at k+1/2 - ξ_above::Grid, # quantity ξ at k-1 - ξ::Grid, # quantity ξ at k - ξ_below::Grid, # quantity ξ at k+1 - Δσₖ::NF # layer thickness on σ levels - ) where {NF<:AbstractFloat,Grid<:AbstractGrid{NF}} - # Δσₖ⁻¹ = 1/Δσₖ # precompute - - # # += as the tendencies already contain the parameterizations - # for ij in eachgridpoint(ξ_tend,σ_tend_above,σ_tend_below,ξ_above,ξ,ξ_below) - # # 1st order upwind scheme - # σ̇⁺ = max(σ_tend_above[ij],0) # velocity into layer k from above - # σ̇⁻ = min(σ_tend_below[ij],0) # velocity into layer k from below - # ξ_tend[ij] -= Δσₖ⁻¹ * (σ̇⁺*(ξ[ij] - ξ_above[ij]) + σ̇⁻*(ξ_below[ij] - ξ[ij])) - # end - - # centred scheme - Δσₖ2⁻¹ = -1/2Δσₖ # precompute - # += as the tendencies already contain the parameterizations - @. ξ_tend += Δσₖ2⁻¹ * (σ_tend_below*(ξ_below - ξ) + σ_tend_above*(ξ - ξ_above)) -end - """ $(TYPEDSIGNATURES) Function barrier to unpack `model`.""" @@ -750,9 +673,14 @@ function SpeedyTransforms.gridded!( diagn::DiagnosticVariablesLayer, (; vor_grid, div_grid, u_grid, v_grid ) = diagn.grid_variables (; temp_grid, humid_grid ) = diagn.grid_variables + (; temp_grid_prev, u_grid_prev, v_grid_prev) = diagn.grid_variables U = diagn.dynamics_variables.a # reuse work arrays for velocities spectral V = diagn.dynamics_variables.b # U = u*coslat, V=v*coslat + @. temp_grid_prev = temp_grid + @. u_grid_prev = u_grid + @. v_grid_prev = v_grid + S = model.spectral_transform wet_core = model isa PrimitiveWet diff --git a/src/dynamics/vertical_advection.jl b/src/dynamics/vertical_advection.jl new file mode 100644 index 000000000..bb6bd8f67 --- /dev/null +++ b/src/dynamics/vertical_advection.jl @@ -0,0 +1,90 @@ +# Dispersive and diffusive advection schemes `NF` is the type, `B` the half-stencil size +abstract type DiffusiveVerticalAdvection{NF, B} <: VerticalAdvection{NF} end +abstract type DispersiveVerticalAdvection{NF, B} <: VerticalAdvection{NF} end + +struct FirstOrderUpwind{NF} <: DiffusiveVerticalAdvection{NF, 1} end +struct ThirdOrderUpwind{NF} <: DiffusiveVerticalAdvection{NF, 2} end +struct WENO3{NF} <: DiffusiveVerticalAdvection{NF, 2} end + +struct SecondOrderCentered{NF} <: DispersiveVerticalAdvection{NF, 1} end +struct FourthOrderCentered{NF} <: DispersiveVerticalAdvection{NF, 2} end + +for T in (:FirstOrderUpwind, :ThirdOrderUpwind, :WENO3, :SecondOrderCentered, :FourthOrderCentered) + @eval $T(spectral_grid::SpectralGrid) = $T{spectral_grid.NF}() +end + +@inline retrieve_time_step(::DiffusiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid_prev)) +@inline retrieve_time_step(::DispersiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid_prev)) + +function vertical_advection!( layer::DiagnosticVariablesLayer, + diagn::DiagnosticVariables, + model::PrimitiveEquation) + + (; k ) = layer # which layer are we on? + + wet_core = model isa PrimitiveWet + (; σ_levels_thick, nlev ) = model.geometry + + vertical_advection = model.vertical_advection + + # for k==1 "above" term is 0, for k==nlev "below" term is zero + # avoid out-of-bounds indexing with k_above, k_below as follows + k_above = max(1,k-1) # just saturate, because M_1/2 = 0 (which zeros that term) + k_below = min(k+1,nlev) # just saturate, because M_nlev+1/2 = 0 (which zeros that term) + + # mass fluxes, M_1/2 = M_nlev+1/2 = 0, but k=1/2 isn't explicitly stored + σ_tend_above = diagn.layers[k_above].dynamics_variables.σ_tend + σ_tend_below = layer.dynamics_variables.σ_tend + + # layer thickness Δσ on level k + Δσₖ = σ_levels_thick[k] + + for var in (:u, :v, :temp) + ξ_tend = getproperty(layer.tendencies, Symbol(var, :_tend_grid)) + ξ_above = retrieve_time_step(vertical_advection, diagn.layers[k_above].grid_variables, var) + ξ = retrieve_time_step(vertical_advection, layer.grid_variables, var) + ξ_below = retrieve_time_step(vertical_advection, diagn.layers[k_below].grid_variables, var) + + _vertical_advection!(ξ_tend,σ_tend_above,σ_tend_below,ξ_above,ξ,ξ_below,Δσₖ,vertical_advection) + end + + if wet_core + ξ_tend = getproperty(layer.tendencies, :humid_tend_grid) + ξ_above = retrieve_time_step(vertical_advection, diagn.layers[k_above].grid_variables, :humid) + ξ = retrieve_time_step(vertical_advection, layer.grid_variables, :humid) + ξ_below = retrieve_time_step(vertical_advection, diagn.layers[k_below].grid_variables, :humid) + + _vertical_advection!(q_tend,σ_tend_above,σ_tend_below,q_above,q,q_below,Δσₖ,vertical_advection) + end +end + +@inline reconstructed_at_face(::FirstOrderUpwind, u, c⁻, c⁺) = ifelse(u > 0, c⁻, c⁺) +@inline reconstructed_at_face(::SecondOrderCentered, u, c⁻, c⁺) = (c⁻ + c⁺) / 2 + +# MULTI THREADED VERSION only writes into layer k +function _vertical_advection!( ξ_tend::Grid, # tendency of quantity ξ at k + σ_tend_above::Grid, # vertical velocity at k-1/2 + σ_tend_below::Grid, # vertical velocity at k+1/2 + ξ_above::Grid, # quantity ξ at k-1 + ξ::Grid, # quantity ξ at k + ξ_below::Grid, # quantity ξ at k+1 + Δσₖ::NF, # layer thickness on σ levels + vertical_advection # vertical advection scheme + ) where {NF<:AbstractFloat,Grid<:AbstractGrid{NF}} + Δσₖ⁻¹ = 1/Δσₖ # precompute + + # += as the tendencies already contain the parameterizations + for ij in eachgridpoint(ξ_tend,σ_tend_above,σ_tend_below,ξ_above,ξ,ξ_below) + σ̇⁻ = σ_tend_above[ij] # velocity into layer k from above + σ̇⁺ = σ_tend_below[ij] # velocity out of layer k to below + + ξₖ₋₁ = ξ_above[ij] + ξₖ = ξ[ij] + ξₖ₊₁ = ξ_below[ij] + + ξ⁺ = reconstructed_at_face(vertical_advection, σ̇⁺, ξₖ, ξₖ₊₁) + ξ⁻ = reconstructed_at_face(vertical_advection, σ̇⁻, ξₖ₋₁, ξₖ) + + ξ_tend[ij] -= Δσₖ⁻¹ * (σ̇⁺ * ξ⁺ - σ̇⁻ * ξ⁻ - ξₖ * (σ̇⁺ - σ̇⁻)) + end +end From 65c0d2e59b5bb381fbef74c709f805e7561fd080 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Thu, 20 Jul 2023 20:03:11 -0400 Subject: [PATCH 02/18] added vertical advection --- src/abstract_types.jl | 2 +- src/dynamics/vertical_advection.jl | 72 ++++++++++++++++++++---------- 2 files changed, 49 insertions(+), 25 deletions(-) diff --git a/src/abstract_types.jl b/src/abstract_types.jl index fb730a348..432ee309d 100644 --- a/src/abstract_types.jl +++ b/src/abstract_types.jl @@ -27,7 +27,7 @@ abstract type AbstractColumnVariables{NF} end abstract type AbstractForcing{NF} end # VERTICAL ADVECTION (PrimitiveEquation) -abstract type VerticalAdvection{NF} end +abstract type VerticalAdvection{NF,B} end # PARAMETERIZATIONS abstract type AbstractParameterization{NF} end diff --git a/src/dynamics/vertical_advection.jl b/src/dynamics/vertical_advection.jl index bb6bd8f67..886a80fe1 100644 --- a/src/dynamics/vertical_advection.jl +++ b/src/dynamics/vertical_advection.jl @@ -1,6 +1,6 @@ # Dispersive and diffusive advection schemes `NF` is the type, `B` the half-stencil size -abstract type DiffusiveVerticalAdvection{NF, B} <: VerticalAdvection{NF} end -abstract type DispersiveVerticalAdvection{NF, B} <: VerticalAdvection{NF} end +abstract type DiffusiveVerticalAdvection{NF, B} <: VerticalAdvection{NF, B} end +abstract type DispersiveVerticalAdvection{NF, B} <: VerticalAdvection{NF, B} end struct FirstOrderUpwind{NF} <: DiffusiveVerticalAdvection{NF, 1} end struct ThirdOrderUpwind{NF} <: DiffusiveVerticalAdvection{NF, 2} end @@ -14,7 +14,9 @@ for T in (:FirstOrderUpwind, :ThirdOrderUpwind, :WENO3, :SecondOrderCentered, :F end @inline retrieve_time_step(::DiffusiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid_prev)) -@inline retrieve_time_step(::DispersiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid_prev)) +@inline retrieve_time_step(::DispersiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid)) + +@inline boundary_buffer(::VerticalAdvection{NF, B}) where {NF, B} = B function vertical_advection!( layer::DiagnosticVariablesLayer, diagn::DiagnosticVariables, @@ -29,11 +31,20 @@ function vertical_advection!( layer::DiagnosticVariablesLayer, # for k==1 "above" term is 0, for k==nlev "below" term is zero # avoid out-of-bounds indexing with k_above, k_below as follows - k_above = max(1,k-1) # just saturate, because M_1/2 = 0 (which zeros that term) - k_below = min(k+1,nlev) # just saturate, because M_nlev+1/2 = 0 (which zeros that term) + # b = boundary_buffer(vertical_advection) + # @. k_stencil = max(min(k-b:k+b, nlev), 1) + + # k_above = max(1,k-1) # just saturate, because M_1/2 = 0 (which zeros that term) + # k_below = min(k+1,nlev) # just saturate, because M_nlev+1/2 = 0 (which zeros that term) + + k⁻ = max(1,k-1) # just saturate, because M_1/2 = 0 (which zeros that term) + k⁺ = min(k+1,nlev) # just saturate, because M_nlev+1/2 = 0 (which zeros that term) + + k⁻⁻ = max(1,k-2) # just saturate, because M_1/2 = 0 (which zeros that term) + k⁺⁺ = min(k+2,nlev) # just saturate, because M_nlev+1/2 = 0 (which zeros that term) # mass fluxes, M_1/2 = M_nlev+1/2 = 0, but k=1/2 isn't explicitly stored - σ_tend_above = diagn.layers[k_above].dynamics_variables.σ_tend + σ_tend_above = diagn.layers[k⁻].dynamics_variables.σ_tend σ_tend_below = layer.dynamics_variables.σ_tend # layer thickness Δσ on level k @@ -41,50 +52,63 @@ function vertical_advection!( layer::DiagnosticVariablesLayer, for var in (:u, :v, :temp) ξ_tend = getproperty(layer.tendencies, Symbol(var, :_tend_grid)) - ξ_above = retrieve_time_step(vertical_advection, diagn.layers[k_above].grid_variables, var) - ξ = retrieve_time_step(vertical_advection, layer.grid_variables, var) - ξ_below = retrieve_time_step(vertical_advection, diagn.layers[k_below].grid_variables, var) + ξ₋₂ = retrieve_time_step(vertical_advection, diagn.layers[k⁻⁻].grid_variables, var) + ξ₋₁ = retrieve_time_step(vertical_advection, diagn.layers[k⁻].grid_variables, var) + ξ = retrieve_time_step(vertical_advection, layer.grid_variables, var) + ξ₊₁ = retrieve_time_step(vertical_advection, diagn.layers[k⁺].grid_variables, var) + ξ₊₂ = retrieve_time_step(vertical_advection, diagn.layers[k⁺⁺].grid_variables, var) - _vertical_advection!(ξ_tend,σ_tend_above,σ_tend_below,ξ_above,ξ,ξ_below,Δσₖ,vertical_advection) + _vertical_advection!(ξ_tend,σ_tend_above,σ_tend_below,ξ₋₂,ξ₋₁,ξ,ξ₊₁,ξ₊₂,Δσₖ,vertical_advection) end if wet_core ξ_tend = getproperty(layer.tendencies, :humid_tend_grid) - ξ_above = retrieve_time_step(vertical_advection, diagn.layers[k_above].grid_variables, :humid) - ξ = retrieve_time_step(vertical_advection, layer.grid_variables, :humid) - ξ_below = retrieve_time_step(vertical_advection, diagn.layers[k_below].grid_variables, :humid) - - _vertical_advection!(q_tend,σ_tend_above,σ_tend_below,q_above,q,q_below,Δσₖ,vertical_advection) + ξ₋₂ = retrieve_time_step(vertical_advection, diagn.layers[k⁻⁻].grid_variables, :humid) + ξ₋₁ = retrieve_time_step(vertical_advection, diagn.layers[k⁻].grid_variables, :humid) + ξ = retrieve_time_step(vertical_advection, layer.grid_variables, :humid) + ξ₊₁ = retrieve_time_step(vertical_advection, diagn.layers[k⁺].grid_variables, :humid) + ξ₊₂ = retrieve_time_step(vertical_advection, diagn.layers[k⁺⁺].grid_variables, :humid) + + _vertical_advection!(ξ_tend,σ_tend_above,σ_tend_below,ξ₋₂,ξ₋₁,ξ,ξ₊₁,ξ₊₂,Δσₖ,vertical_advection) end end -@inline reconstructed_at_face(::FirstOrderUpwind, u, c⁻, c⁺) = ifelse(u > 0, c⁻, c⁺) -@inline reconstructed_at_face(::SecondOrderCentered, u, c⁻, c⁺) = (c⁻ + c⁺) / 2 +@inline reconstructed_at_face(::FirstOrderUpwind, u, c⁻⁻, c⁻, c⁺, c⁺⁺) = ifelse(u > 0, c⁻, c⁺) +@inline reconstructed_at_face(::ThirdOrderUpwind, u, c⁻⁻, c⁻, c⁺, c⁺⁺) = ifelse(u > 0, (2c⁻⁻ + 5c⁻ - c⁺ ) / 6, + (-c⁻ + 5c⁺ + 2c⁺⁺) / 6) + +@inline reconstructed_at_face(::SecondOrderCentered, u, c⁻⁻, c⁻, c⁺, c⁺⁺) = (c⁻ + c⁺) / 2 +@inline reconstructed_at_face(::FourthOrderCentered, u, c⁻⁻, c⁻, c⁺, c⁺⁺) = (- c⁻⁻ + 7c⁻ + 7c⁺ - c⁺⁺) / 12 # MULTI THREADED VERSION only writes into layer k function _vertical_advection!( ξ_tend::Grid, # tendency of quantity ξ at k σ_tend_above::Grid, # vertical velocity at k-1/2 σ_tend_below::Grid, # vertical velocity at k+1/2 - ξ_above::Grid, # quantity ξ at k-1 + ξ₋₂::Grid, # quantity ξ at k-2 + ξ₋₁::Grid, # quantity ξ at k-1 ξ::Grid, # quantity ξ at k - ξ_below::Grid, # quantity ξ at k+1 + ξ₊₁::Grid, # quantity ξ at k+1 + ξ₊₂::Grid, # quantity ξ at k+2 Δσₖ::NF, # layer thickness on σ levels vertical_advection # vertical advection scheme ) where {NF<:AbstractFloat,Grid<:AbstractGrid{NF}} Δσₖ⁻¹ = 1/Δσₖ # precompute # += as the tendencies already contain the parameterizations - for ij in eachgridpoint(ξ_tend,σ_tend_above,σ_tend_below,ξ_above,ξ,ξ_below) + for ij in eachgridpoint(ξ_tend) σ̇⁻ = σ_tend_above[ij] # velocity into layer k from above σ̇⁺ = σ_tend_below[ij] # velocity out of layer k to below - ξₖ₋₁ = ξ_above[ij] + ξₖ₋₂ = ξ₋₂[ij] + ξₖ₋₁ = ξ₋₁[ij] ξₖ = ξ[ij] - ξₖ₊₁ = ξ_below[ij] + ξₖ₊₁ = ξ₊₁[ij] + ξₖ₊₂ = ξ₊₂[ij] - ξ⁺ = reconstructed_at_face(vertical_advection, σ̇⁺, ξₖ, ξₖ₊₁) - ξ⁻ = reconstructed_at_face(vertical_advection, σ̇⁻, ξₖ₋₁, ξₖ) + ξ⁺ = reconstructed_at_face(vertical_advection, σ̇⁺, ξₖ₋₁, ξₖ, ξₖ₊₁, ξₖ₊₂) + ξ⁻ = reconstructed_at_face(vertical_advection, σ̇⁻, ξₖ₋₂, ξₖ₋₁, ξₖ, ξₖ₊₁) ξ_tend[ij] -= Δσₖ⁻¹ * (σ̇⁺ * ξ⁺ - σ̇⁻ * ξ⁻ - ξₖ * (σ̇⁺ - σ̇⁻)) end end + From efde3ec000606e539c86fe3098780cb9e621c21e Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Fri, 21 Jul 2023 09:52:35 -0400 Subject: [PATCH 03/18] cleaner --- src/dynamics/vertical_advection.jl | 27 ++++++++++++--------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/src/dynamics/vertical_advection.jl b/src/dynamics/vertical_advection.jl index 886a80fe1..c74eb8dd2 100644 --- a/src/dynamics/vertical_advection.jl +++ b/src/dynamics/vertical_advection.jl @@ -2,21 +2,18 @@ abstract type DiffusiveVerticalAdvection{NF, B} <: VerticalAdvection{NF, B} end abstract type DispersiveVerticalAdvection{NF, B} <: VerticalAdvection{NF, B} end -struct FirstOrderUpwind{NF} <: DiffusiveVerticalAdvection{NF, 1} end -struct ThirdOrderUpwind{NF} <: DiffusiveVerticalAdvection{NF, 2} end -struct WENO3{NF} <: DiffusiveVerticalAdvection{NF, 2} end +struct UpwindVerticalAdvection{NF, B} <: DiffusiveVerticalAdvection{NF, B} end +struct WENOVerticalAdvection{NF, B} <: DiffusiveVerticalAdvection{NF, B} end +struct CenteredVerticalAdvection{NF, B} <: DispersiveVerticalAdvection{NF, B} end -struct SecondOrderCentered{NF} <: DispersiveVerticalAdvection{NF, 1} end -struct FourthOrderCentered{NF} <: DispersiveVerticalAdvection{NF, 2} end - -for T in (:FirstOrderUpwind, :ThirdOrderUpwind, :WENO3, :SecondOrderCentered, :FourthOrderCentered) - @eval $T(spectral_grid::SpectralGrid) = $T{spectral_grid.NF}() -end +CenteredVerticalAdvection(spectral_grid; order = 2) = CenteredVerticalAdvection{spectral_grid.NF, order}() +UpwindVerticalAdvection(spectral_grid; order = 3) = UpwindVerticalAdvection{spectral_grid.NF, order}() +WENOVerticalAdvection(spectral_grid; order = 3) = WENOVerticalAdvection{spectral_grid.NF, order}() @inline retrieve_time_step(::DiffusiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid_prev)) @inline retrieve_time_step(::DispersiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid)) -@inline boundary_buffer(::VerticalAdvection{NF, B}) where {NF, B} = B +@inline order(::VerticalAdvection{NF, B}) where {NF, B} = B function vertical_advection!( layer::DiagnosticVariablesLayer, diagn::DiagnosticVariables, @@ -73,12 +70,12 @@ function vertical_advection!( layer::DiagnosticVariablesLayer, end end -@inline reconstructed_at_face(::FirstOrderUpwind, u, c⁻⁻, c⁻, c⁺, c⁺⁺) = ifelse(u > 0, c⁻, c⁺) -@inline reconstructed_at_face(::ThirdOrderUpwind, u, c⁻⁻, c⁻, c⁺, c⁺⁺) = ifelse(u > 0, (2c⁻⁻ + 5c⁻ - c⁺ ) / 6, - (-c⁻ + 5c⁺ + 2c⁺⁺) / 6) +@inline reconstructed_at_face(::UpwindVerticalAdvection{NF, 1}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = ifelse(u > 0, c⁻, c⁺) +@inline reconstructed_at_face(::UpwindVerticalAdvection{NF, 3}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = ifelse(u > 0, (2c⁻⁻ + 5c⁻ - c⁺ ) / 6, + (-c⁻ + 5c⁺ + 2c⁺⁺) / 6) -@inline reconstructed_at_face(::SecondOrderCentered, u, c⁻⁻, c⁻, c⁺, c⁺⁺) = (c⁻ + c⁺) / 2 -@inline reconstructed_at_face(::FourthOrderCentered, u, c⁻⁻, c⁻, c⁺, c⁺⁺) = (- c⁻⁻ + 7c⁻ + 7c⁺ - c⁺⁺) / 12 +@inline reconstructed_at_face(::CenteredVerticalAdvection{NF, 2}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) = (c⁻ + c⁺) / 2 +@inline reconstructed_at_face(::CenteredVerticalAdvection{NF, 3}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) = (- c⁻⁻ + 7c⁻ + 7c⁺ - c⁺⁺) / 12 # MULTI THREADED VERSION only writes into layer k function _vertical_advection!( ξ_tend::Grid, # tendency of quantity ξ at k From df7eac27c1e547218bd0e768fc6a51c1fba7e082 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Fri, 21 Jul 2023 11:12:34 -0400 Subject: [PATCH 04/18] bugfix --- src/dynamics/vertical_advection.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dynamics/vertical_advection.jl b/src/dynamics/vertical_advection.jl index c74eb8dd2..8b9e7ec38 100644 --- a/src/dynamics/vertical_advection.jl +++ b/src/dynamics/vertical_advection.jl @@ -74,8 +74,8 @@ end @inline reconstructed_at_face(::UpwindVerticalAdvection{NF, 3}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = ifelse(u > 0, (2c⁻⁻ + 5c⁻ - c⁺ ) / 6, (-c⁻ + 5c⁺ + 2c⁺⁺) / 6) -@inline reconstructed_at_face(::CenteredVerticalAdvection{NF, 2}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) = (c⁻ + c⁺) / 2 -@inline reconstructed_at_face(::CenteredVerticalAdvection{NF, 3}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) = (- c⁻⁻ + 7c⁻ + 7c⁺ - c⁺⁺) / 12 +@inline reconstructed_at_face(::CenteredVerticalAdvection{NF, 2}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = (c⁻ + c⁺) / 2 +@inline reconstructed_at_face(::CenteredVerticalAdvection{NF, 3}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = (- c⁻⁻ + 7c⁻ + 7c⁺ - c⁺⁺) / 12 # MULTI THREADED VERSION only writes into layer k function _vertical_advection!( ξ_tend::Grid, # tendency of quantity ξ at k From 4e8ccd8dae83a8b393a4e776057c975d0d5505a7 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Fri, 21 Jul 2023 18:16:24 -0400 Subject: [PATCH 05/18] remove WENO3 for the moment --- src/dynamics/vertical_advection.jl | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/src/dynamics/vertical_advection.jl b/src/dynamics/vertical_advection.jl index 8b9e7ec38..5e1fa9ff9 100644 --- a/src/dynamics/vertical_advection.jl +++ b/src/dynamics/vertical_advection.jl @@ -3,12 +3,10 @@ abstract type DiffusiveVerticalAdvection{NF, B} <: VerticalAdvection{NF, B} end abstract type DispersiveVerticalAdvection{NF, B} <: VerticalAdvection{NF, B} end struct UpwindVerticalAdvection{NF, B} <: DiffusiveVerticalAdvection{NF, B} end -struct WENOVerticalAdvection{NF, B} <: DiffusiveVerticalAdvection{NF, B} end struct CenteredVerticalAdvection{NF, B} <: DispersiveVerticalAdvection{NF, B} end CenteredVerticalAdvection(spectral_grid; order = 2) = CenteredVerticalAdvection{spectral_grid.NF, order}() UpwindVerticalAdvection(spectral_grid; order = 3) = UpwindVerticalAdvection{spectral_grid.NF, order}() -WENOVerticalAdvection(spectral_grid; order = 3) = WENOVerticalAdvection{spectral_grid.NF, order}() @inline retrieve_time_step(::DiffusiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid_prev)) @inline retrieve_time_step(::DispersiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid)) @@ -70,13 +68,6 @@ function vertical_advection!( layer::DiagnosticVariablesLayer, end end -@inline reconstructed_at_face(::UpwindVerticalAdvection{NF, 1}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = ifelse(u > 0, c⁻, c⁺) -@inline reconstructed_at_face(::UpwindVerticalAdvection{NF, 3}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = ifelse(u > 0, (2c⁻⁻ + 5c⁻ - c⁺ ) / 6, - (-c⁻ + 5c⁺ + 2c⁺⁺) / 6) - -@inline reconstructed_at_face(::CenteredVerticalAdvection{NF, 2}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = (c⁻ + c⁺) / 2 -@inline reconstructed_at_face(::CenteredVerticalAdvection{NF, 3}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = (- c⁻⁻ + 7c⁻ + 7c⁺ - c⁺⁺) / 12 - # MULTI THREADED VERSION only writes into layer k function _vertical_advection!( ξ_tend::Grid, # tendency of quantity ξ at k σ_tend_above::Grid, # vertical velocity at k-1/2 @@ -109,3 +100,9 @@ function _vertical_advection!( ξ_tend::Grid, # tendency of quantity end end +@inline reconstructed_at_face(::UpwindVerticalAdvection{NF, 1}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = ifelse(u > 0, c⁻, c⁺) +@inline reconstructed_at_face(::UpwindVerticalAdvection{NF, 3}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = ifelse(u > 0, (2c⁻⁻ + 5c⁻ - c⁺ ) / 6, + (-c⁻ + 5c⁺ + 2c⁺⁺) / 6) + +@inline reconstructed_at_face(::CenteredVerticalAdvection{NF, 2}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = (c⁻ + c⁺) / 2 +@inline reconstructed_at_face(::CenteredVerticalAdvection{NF, 3}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = (- c⁻⁻ + 7c⁻ + 7c⁺ - c⁺⁺) / 12 From 4cede889435bf053d38615093252912a01a339e5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Mon, 24 Jul 2023 19:07:50 -0400 Subject: [PATCH 06/18] generalized schemes --- src/dynamics/diagnostic_variables.jl | 6 +- src/dynamics/models.jl | 2 +- src/dynamics/vertical_advection.jl | 87 +++++++++++----------------- 3 files changed, 39 insertions(+), 56 deletions(-) diff --git a/src/dynamics/diagnostic_variables.jl b/src/dynamics/diagnostic_variables.jl index 37b1d89d6..a256c4d64 100644 --- a/src/dynamics/diagnostic_variables.jl +++ b/src/dynamics/diagnostic_variables.jl @@ -46,14 +46,14 @@ struct GridVariables{NF<:AbstractFloat,Grid<:AbstractGrid{NF}} vor_grid ::Grid # vorticity div_grid ::Grid # divergence temp_grid ::Grid # absolute temperature [K] - temp_grid_prev ::Grid + temp_grid_prev ::Grid # absolute temperature of previous time step [K] temp_virt_grid ::Grid # virtual tempereature [K] humid_grid ::Grid # specific_humidity geopot_grid ::Grid # geopotential (is that needed?) u_grid ::Grid # zonal velocity *coslat [m/s] v_grid ::Grid # meridional velocity *coslat [m/s] - u_grid_prev ::Grid # zonal velocity *coslat [m/s] - v_grid_prev ::Grid # meridional velocity *coslat [m/s] + u_grid_prev ::Grid # zonal velocity *coslat of previous time step [m/s] + v_grid_prev ::Grid # meridional velocity *coslat of previous time step [m/s] end function Base.zeros(::Type{GridVariables},SG::SpectralGrid) diff --git a/src/dynamics/models.jl b/src/dynamics/models.jl index 0bd1ec87e..0badb8216 100644 --- a/src/dynamics/models.jl +++ b/src/dynamics/models.jl @@ -207,7 +207,7 @@ Base.@kwdef struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid) static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid) large_scale_condensation::AbstractCondensation{NF} = SpeedyCondensation(spectral_grid) - vertical_advection::VerticalAdvection{NF} = FirstOrderUpwind(spectral_grid) + vertical_advection::VerticalAdvection{NF} = UpwindVerticalAdvection(spectral_grid) # vertical_diffusion::VerticalDiffusion{NF} = VerticalLaplacian(spectral_grid) diff --git a/src/dynamics/vertical_advection.jl b/src/dynamics/vertical_advection.jl index 5e1fa9ff9..1f34f8ee2 100644 --- a/src/dynamics/vertical_advection.jl +++ b/src/dynamics/vertical_advection.jl @@ -5,14 +5,20 @@ abstract type DispersiveVerticalAdvection{NF, B} <: VerticalAdvection{NF, B} end struct UpwindVerticalAdvection{NF, B} <: DiffusiveVerticalAdvection{NF, B} end struct CenteredVerticalAdvection{NF, B} <: DispersiveVerticalAdvection{NF, B} end -CenteredVerticalAdvection(spectral_grid; order = 2) = CenteredVerticalAdvection{spectral_grid.NF, order}() -UpwindVerticalAdvection(spectral_grid; order = 3) = UpwindVerticalAdvection{spectral_grid.NF, order}() +CenteredVerticalAdvection(spectral_grid; order = 2) = CenteredVerticalAdvection{spectral_grid.NF, order ÷ 2}() +UpwindVerticalAdvection(spectral_grid; order = 5) = UpwindVerticalAdvection{spectral_grid.NF, (order + 1) ÷ 2}() @inline retrieve_time_step(::DiffusiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid_prev)) @inline retrieve_time_step(::DispersiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid)) @inline order(::VerticalAdvection{NF, B}) where {NF, B} = B +@inline function retrieve_stencil(k, layers, var, nlev, vertical_advection::VerticalAdvection{NF, B}) where {NF, B} + k_stencil = max.(min.(nlev, k-B:k+B), 1) + ξ_stencil = Tuple(retrieve_time_step(vertical_advection, layers[k].grid_variables, var) for k in k_stencil) + return ξ_stencil +end + function vertical_advection!( layer::DiagnosticVariablesLayer, diagn::DiagnosticVariables, model::PrimitiveEquation) @@ -26,17 +32,7 @@ function vertical_advection!( layer::DiagnosticVariablesLayer, # for k==1 "above" term is 0, for k==nlev "below" term is zero # avoid out-of-bounds indexing with k_above, k_below as follows - # b = boundary_buffer(vertical_advection) - # @. k_stencil = max(min(k-b:k+b, nlev), 1) - - # k_above = max(1,k-1) # just saturate, because M_1/2 = 0 (which zeros that term) - # k_below = min(k+1,nlev) # just saturate, because M_nlev+1/2 = 0 (which zeros that term) - - k⁻ = max(1,k-1) # just saturate, because M_1/2 = 0 (which zeros that term) - k⁺ = min(k+1,nlev) # just saturate, because M_nlev+1/2 = 0 (which zeros that term) - - k⁻⁻ = max(1,k-2) # just saturate, because M_1/2 = 0 (which zeros that term) - k⁺⁺ = min(k+2,nlev) # just saturate, because M_nlev+1/2 = 0 (which zeros that term) + k⁻ = max(1,k-1) # just saturate, because M_1/2 = 0 (which zeros that term) # mass fluxes, M_1/2 = M_nlev+1/2 = 0, but k=1/2 isn't explicitly stored σ_tend_above = diagn.layers[k⁻].dynamics_variables.σ_tend @@ -46,40 +42,31 @@ function vertical_advection!( layer::DiagnosticVariablesLayer, Δσₖ = σ_levels_thick[k] for var in (:u, :v, :temp) - ξ_tend = getproperty(layer.tendencies, Symbol(var, :_tend_grid)) - ξ₋₂ = retrieve_time_step(vertical_advection, diagn.layers[k⁻⁻].grid_variables, var) - ξ₋₁ = retrieve_time_step(vertical_advection, diagn.layers[k⁻].grid_variables, var) - ξ = retrieve_time_step(vertical_advection, layer.grid_variables, var) - ξ₊₁ = retrieve_time_step(vertical_advection, diagn.layers[k⁺].grid_variables, var) - ξ₊₂ = retrieve_time_step(vertical_advection, diagn.layers[k⁺⁺].grid_variables, var) - - _vertical_advection!(ξ_tend,σ_tend_above,σ_tend_below,ξ₋₂,ξ₋₁,ξ,ξ₊₁,ξ₊₂,Δσₖ,vertical_advection) + ξ_tend = getproperty(layer.tendencies, Symbol(var, :_tend_grid)) + ξ_sten = retrieve_stencil(k, diagn.layers, var, nlev, vertical_advection) + ξ = retrieve_time_step(vertical_advection, layer.grid_variables, var) + + _vertical_advection!(ξ_tend,σ_tend_above,σ_tend_below,ξ_sten,ξ,Δσₖ,vertical_advection) end if wet_core ξ_tend = getproperty(layer.tendencies, :humid_tend_grid) - ξ₋₂ = retrieve_time_step(vertical_advection, diagn.layers[k⁻⁻].grid_variables, :humid) - ξ₋₁ = retrieve_time_step(vertical_advection, diagn.layers[k⁻].grid_variables, :humid) - ξ = retrieve_time_step(vertical_advection, layer.grid_variables, :humid) - ξ₊₁ = retrieve_time_step(vertical_advection, diagn.layers[k⁺].grid_variables, :humid) - ξ₊₂ = retrieve_time_step(vertical_advection, diagn.layers[k⁺⁺].grid_variables, :humid) + ξ_sten = retrieve_stencil(k, diagn.layers, :humid, nlev, vertical_advection) + ξ = retrieve_time_step(vertical_advection, layer.grid_variables, :humid) - _vertical_advection!(ξ_tend,σ_tend_above,σ_tend_below,ξ₋₂,ξ₋₁,ξ,ξ₊₁,ξ₊₂,Δσₖ,vertical_advection) + _vertical_advection!(ξ_tend,σ_tend_above,σ_tend_below,ξ_sten,ξ,Δσₖ,vertical_advection) end end # MULTI THREADED VERSION only writes into layer k -function _vertical_advection!( ξ_tend::Grid, # tendency of quantity ξ at k - σ_tend_above::Grid, # vertical velocity at k-1/2 - σ_tend_below::Grid, # vertical velocity at k+1/2 - ξ₋₂::Grid, # quantity ξ at k-2 - ξ₋₁::Grid, # quantity ξ at k-1 - ξ::Grid, # quantity ξ at k - ξ₊₁::Grid, # quantity ξ at k+1 - ξ₊₂::Grid, # quantity ξ at k+2 - Δσₖ::NF, # layer thickness on σ levels - vertical_advection # vertical advection scheme - ) where {NF<:AbstractFloat,Grid<:AbstractGrid{NF}} +function _vertical_advection!( ξ_tend::Grid, # tendency of quantity ξ at k + σ_tend_above::Grid, # vertical velocity at k-1/2 + σ_tend_below::Grid, # vertical velocity at k+1/2 + ξ_sten, # ξ stencil for vertical advection (from k-B to k+B) + ξ::Grid, # ξ at level k + Δσₖ::NF, # layer thickness on σ levels + adv::VerticalAdvection{NF, B} # vertical advection scheme + ) where {NF<:AbstractFloat,Grid<:AbstractGrid{NF}, B} Δσₖ⁻¹ = 1/Δσₖ # precompute # += as the tendencies already contain the parameterizations @@ -87,22 +74,18 @@ function _vertical_advection!( ξ_tend::Grid, # tendency of quantity σ̇⁻ = σ_tend_above[ij] # velocity into layer k from above σ̇⁺ = σ_tend_below[ij] # velocity out of layer k to below - ξₖ₋₂ = ξ₋₂[ij] - ξₖ₋₁ = ξ₋₁[ij] - ξₖ = ξ[ij] - ξₖ₊₁ = ξ₊₁[ij] - ξₖ₊₂ = ξ₊₂[ij] - - ξ⁺ = reconstructed_at_face(vertical_advection, σ̇⁺, ξₖ₋₁, ξₖ, ξₖ₊₁, ξₖ₊₂) - ξ⁻ = reconstructed_at_face(vertical_advection, σ̇⁻, ξₖ₋₂, ξₖ₋₁, ξₖ, ξₖ₊₁) + ξ⁺ = reconstructed_at_face(ij, adv, σ̇⁺, ξ_sten[2:end]) + ξ⁻ = reconstructed_at_face(ij, adv, σ̇⁻, ξ_sten[1:end-1]) - ξ_tend[ij] -= Δσₖ⁻¹ * (σ̇⁺ * ξ⁺ - σ̇⁻ * ξ⁻ - ξₖ * (σ̇⁺ - σ̇⁻)) + ξ_tend[ij] -= Δσₖ⁻¹ * (σ̇⁺ * ξ⁺ - σ̇⁻ * ξ⁻ - ξ[ij] * (σ̇⁺ - σ̇⁻)) end end -@inline reconstructed_at_face(::UpwindVerticalAdvection{NF, 1}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = ifelse(u > 0, c⁻, c⁺) -@inline reconstructed_at_face(::UpwindVerticalAdvection{NF, 3}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = ifelse(u > 0, (2c⁻⁻ + 5c⁻ - c⁺ ) / 6, - (-c⁻ + 5c⁺ + 2c⁺⁺) / 6) +@inline reconstructed_at_face(ij, ::UpwindVerticalAdvection{NF, 1}, u, ξ) where NF = ifelse(u > 0, ξ[1][ij], ξ[2][ij]) +@inline reconstructed_at_face(ij, ::UpwindVerticalAdvection{NF, 2}, u, ξ) where NF = ifelse(u > 0, (2ξ[1][ij] + 5ξ[2][ij] - ξ[3][ij]) / 6, + (-ξ[2][ij] + 5ξ[3][ij] + 2ξ[4][ij]) / 6) +@inline reconstructed_at_face(ij, ::UpwindVerticalAdvection{NF, 3}, u, ξ) where NF = ifelse(u > 0, ( 2ξ[1][ij] - 13ξ[2][ij] + 47ξ[3][ij] + 27ξ[4][ij] - 3ξ[5][ij]) / 60, + (-3ξ[2][ij] + 27ξ[3][ij] + 47ξ[4][ij] - 13ξ[5][ij] + 2ξ[6][ij]) / 60) -@inline reconstructed_at_face(::CenteredVerticalAdvection{NF, 2}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = (c⁻ + c⁺) / 2 -@inline reconstructed_at_face(::CenteredVerticalAdvection{NF, 3}, u, c⁻⁻, c⁻, c⁺, c⁺⁺) where NF = (- c⁻⁻ + 7c⁻ + 7c⁺ - c⁺⁺) / 12 +@inline reconstructed_at_face(ij, ::CenteredVerticalAdvection{NF, 1}, u, ξ) where NF = ( ξ[1][ij] + ξ[2][ij]) / 2 +@inline reconstructed_at_face(ij, ::CenteredVerticalAdvection{NF, 2}, u, ξ) where NF = (-ξ[1][ij] + 7ξ[2][ij] + 7ξ[3][ij] - ξ[4][ij]) / 12 From 346c51bd76ffa89aca9bef93a8e90c70b466a626 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Mon, 24 Jul 2023 19:25:01 -0400 Subject: [PATCH 07/18] fixed test --- .gitignore | 3 +++ src/dynamics/models.jl | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index e84627101..3cc979bbb 100644 --- a/.gitignore +++ b/.gitignore @@ -30,3 +30,6 @@ Manifest.toml # no model output folders run_*/ + +# no video outputs +*.mp4 diff --git a/src/dynamics/models.jl b/src/dynamics/models.jl index 0badb8216..41b04807e 100644 --- a/src/dynamics/models.jl +++ b/src/dynamics/models.jl @@ -139,7 +139,7 @@ Base.@kwdef struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr boundary_layer_drag::BoundaryLayerDrag{NF} = LinearDrag(spectral_grid) temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid) static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid) - vertical_advection::VerticalAdvection{NF} = FirstOrderUpwind(spectral_grid) + vertical_advection::VerticalAdvection{NF} = UpwindVerticalAdvection(spectral_grid) # vertical_diffusion::VerticalDiffusion{NF} = VerticalLaplacian(spectral_grid) # NUMERICS From 17481c01103227f2bd74b8ad76b4540f98370da9 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Mon, 24 Jul 2023 20:21:33 -0400 Subject: [PATCH 08/18] Weno fifth order --- src/dynamics/vertical_advection.jl | 42 ++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/src/dynamics/vertical_advection.jl b/src/dynamics/vertical_advection.jl index 1f34f8ee2..bc64f1b0d 100644 --- a/src/dynamics/vertical_advection.jl +++ b/src/dynamics/vertical_advection.jl @@ -3,10 +3,12 @@ abstract type DiffusiveVerticalAdvection{NF, B} <: VerticalAdvection{NF, B} end abstract type DispersiveVerticalAdvection{NF, B} <: VerticalAdvection{NF, B} end struct UpwindVerticalAdvection{NF, B} <: DiffusiveVerticalAdvection{NF, B} end +struct WENOVerticalAdvection{NF} <: DiffusiveVerticalAdvection{NF, 3} end struct CenteredVerticalAdvection{NF, B} <: DispersiveVerticalAdvection{NF, B} end CenteredVerticalAdvection(spectral_grid; order = 2) = CenteredVerticalAdvection{spectral_grid.NF, order ÷ 2}() UpwindVerticalAdvection(spectral_grid; order = 5) = UpwindVerticalAdvection{spectral_grid.NF, (order + 1) ÷ 2}() +WENOVerticalAdvection(spectral_grid) = WENOVerticalAdvection{spectral_grid.NF}() @inline retrieve_time_step(::DiffusiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid_prev)) @inline retrieve_time_step(::DispersiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid)) @@ -89,3 +91,43 @@ end @inline reconstructed_at_face(ij, ::CenteredVerticalAdvection{NF, 1}, u, ξ) where NF = ( ξ[1][ij] + ξ[2][ij]) / 2 @inline reconstructed_at_face(ij, ::CenteredVerticalAdvection{NF, 2}, u, ξ) where NF = (-ξ[1][ij] + 7ξ[2][ij] + 7ξ[3][ij] - ξ[4][ij]) / 12 + +const ε = 1e-6 +const d₀ = 3/10 +const d₁ = 3/5 +const d₂ = 1/10 + +@inline weight_β₀(S, NF) = NF(13/12) * (S[1] - 2S[2] + S[3])^2 + NF(1/4) * (3S[1] - 4S[2] + S[3])^2 +@inline weight_β₁(S, NF) = NF(13/12) * (S[1] - 2S[2] + S[3])^2 + NF(1/4) * ( S[1] - S[3])^2 +@inline weight_β₂(S, NF) = NF(13/12) * (S[1] - 2S[2] + S[3])^2 + NF(1/4) * ( S[1] - 4S[2] + 3S[3])^2 + +@inline p₀(S) = (2S[1] + 5S[2] - S[3]) / 6 +@inline p₁(S) = (-S[1] + 5S[2] + 2S[3]) / 6 +@inline p₂(S) = (2S[1] - 7S[2] + 11S[3]) / 6 + +@inline function weno_reconstruction(S₀, S₁, S₂) + β₀ = weight_β₀(S₀, NF) + β₁ = weight_β₁(S₁, NF) + β₂ = weight_β₂(S₂, NF) + + w₀ = NF(d₀) / (β₀ + NF(ε))^2 + w₁ = NF(d₁) / (β₁ + NF(ε))^2 + w₂ = NF(d₂) / (β₂ + NF(ε))^2 + + w₀, w₁, w₂ = (w₀, w₁, w₂) ./ (w₀ + w₁ + w₂) + + return p₀(S₀) * w₀ + p₁(S₁) * w₁ + p₂(S₂) * w₂ +end + +@inline function reconstructed_at_face(ij, ::WENOVerticalAdvection{NF}, u, ξ) where NF + if u > 0 + S₀ = (ξ[1][ij], ξ[2][ij], ξ[3][ij]) + S₁ = (ξ[2][ij], ξ[3][ij], ξ[4][ij]) + S₂ = (ξ[3][ij], ξ[4][ij], ξ[5][ij]) + else + S₀ = (ξ[6][ij], ξ[5][ij], ξ[4][ij]) + S₁ = (ξ[5][ij], ξ[4][ij], ξ[3][ij]) + S₂ = (ξ[4][ij], ξ[3][ij], ξ[2][ij]) + end + return weno_reconstruction(S₀, S₁, S₂) +end From 680ad53b134745aa79245886cf26b7d37ce7a011 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Mon, 24 Jul 2023 20:21:41 -0400 Subject: [PATCH 09/18] weno 5th order --- src/dynamics/vertical_advection.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dynamics/vertical_advection.jl b/src/dynamics/vertical_advection.jl index bc64f1b0d..4cc4df6af 100644 --- a/src/dynamics/vertical_advection.jl +++ b/src/dynamics/vertical_advection.jl @@ -105,7 +105,7 @@ const d₂ = 1/10 @inline p₁(S) = (-S[1] + 5S[2] + 2S[3]) / 6 @inline p₂(S) = (2S[1] - 7S[2] + 11S[3]) / 6 -@inline function weno_reconstruction(S₀, S₁, S₂) +@inline function weno_reconstruction(S₀, S₁, S₂, NF) β₀ = weight_β₀(S₀, NF) β₁ = weight_β₁(S₁, NF) β₂ = weight_β₂(S₂, NF) @@ -129,5 +129,5 @@ end S₁ = (ξ[5][ij], ξ[4][ij], ξ[3][ij]) S₂ = (ξ[4][ij], ξ[3][ij], ξ[2][ij]) end - return weno_reconstruction(S₀, S₁, S₂) + return weno_reconstruction(S₀, S₁, S₂, NF) end From c69d139640d65db9e5203c6b8bf467010e19a371 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Mon, 24 Jul 2023 20:22:40 -0400 Subject: [PATCH 10/18] test fixed --- src/dynamics/diagnostic_variables.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/dynamics/diagnostic_variables.jl b/src/dynamics/diagnostic_variables.jl index a256c4d64..b18323739 100644 --- a/src/dynamics/diagnostic_variables.jl +++ b/src/dynamics/diagnostic_variables.jl @@ -49,6 +49,7 @@ struct GridVariables{NF<:AbstractFloat,Grid<:AbstractGrid{NF}} temp_grid_prev ::Grid # absolute temperature of previous time step [K] temp_virt_grid ::Grid # virtual tempereature [K] humid_grid ::Grid # specific_humidity + humid_grid_prev ::Grid # specific_humidity of previous time step geopot_grid ::Grid # geopotential (is that needed?) u_grid ::Grid # zonal velocity *coslat [m/s] v_grid ::Grid # meridional velocity *coslat [m/s] @@ -65,13 +66,14 @@ function Base.zeros(::Type{GridVariables},SG::SpectralGrid) temp_grid_prev = zeros(Grid{NF},nlat_half) # absolute temperature temp_virt_grid = zeros(Grid{NF},nlat_half) # virtual temperature humid_grid = zeros(Grid{NF},nlat_half) # specific humidity + humid_grid_prev = zeros(Grid{NF},nlat_half) # specific humidity geopot_grid = zeros(Grid{NF},nlat_half) # geopotential u_grid = zeros(Grid{NF},nlat_half) # zonal velocity *coslat v_grid = zeros(Grid{NF},nlat_half) # meridonal velocity *coslat u_grid_prev = zeros(Grid{NF},nlat_half) # zonal velocity *coslat v_grid_prev = zeros(Grid{NF},nlat_half) # meridonal velocity *coslat - return GridVariables( vor_grid,div_grid,temp_grid,temp_grid_prev,temp_virt_grid,humid_grid,geopot_grid, + return GridVariables( vor_grid,div_grid,temp_grid,temp_grid_prev,temp_virt_grid,humid_grid,humid_grid_prev,geopot_grid, u_grid,v_grid,u_grid_prev,v_grid_prev) end From 7d3bc7817be8bd24b2f180436874b1fe8c56d4bb Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Mon, 24 Jul 2023 20:38:48 -0400 Subject: [PATCH 11/18] corrected WENO --- src/dynamics/vertical_advection.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/dynamics/vertical_advection.jl b/src/dynamics/vertical_advection.jl index 4cc4df6af..7fb86aad8 100644 --- a/src/dynamics/vertical_advection.jl +++ b/src/dynamics/vertical_advection.jl @@ -121,13 +121,13 @@ end @inline function reconstructed_at_face(ij, ::WENOVerticalAdvection{NF}, u, ξ) where NF if u > 0 - S₀ = (ξ[1][ij], ξ[2][ij], ξ[3][ij]) + S₀ = (ξ[3][ij], ξ[4][ij], ξ[5][ij]) S₁ = (ξ[2][ij], ξ[3][ij], ξ[4][ij]) - S₂ = (ξ[3][ij], ξ[4][ij], ξ[5][ij]) + S₂ = (ξ[1][ij], ξ[2][ij], ξ[3][ij]) else - S₀ = (ξ[6][ij], ξ[5][ij], ξ[4][ij]) + S₀ = (ξ[4][ij], ξ[3][ij], ξ[2][ij]) S₁ = (ξ[5][ij], ξ[4][ij], ξ[3][ij]) - S₂ = (ξ[4][ij], ξ[3][ij], ξ[2][ij]) + S₂ = (ξ[6][ij], ξ[5][ij], ξ[4][ij]) end return weno_reconstruction(S₀, S₁, S₂, NF) end From 6fcde0648ee2bf3114b5fc3a72a447ec24bcf8e8 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Tue, 25 Jul 2023 11:42:57 -0400 Subject: [PATCH 12/18] cleaner --- src/dynamics/vertical_advection.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/dynamics/vertical_advection.jl b/src/dynamics/vertical_advection.jl index 7fb86aad8..f6fd851b5 100644 --- a/src/dynamics/vertical_advection.jl +++ b/src/dynamics/vertical_advection.jl @@ -8,7 +8,7 @@ struct CenteredVerticalAdvection{NF, B} <: DispersiveVerticalAdvection{NF, B} en CenteredVerticalAdvection(spectral_grid; order = 2) = CenteredVerticalAdvection{spectral_grid.NF, order ÷ 2}() UpwindVerticalAdvection(spectral_grid; order = 5) = UpwindVerticalAdvection{spectral_grid.NF, (order + 1) ÷ 2}() -WENOVerticalAdvection(spectral_grid) = WENOVerticalAdvection{spectral_grid.NF}() +WENOVerticalAdvection(spectral_grid) = WENOVerticalAdvection{spectral_grid.NF}() @inline retrieve_time_step(::DiffusiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid_prev)) @inline retrieve_time_step(::DispersiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid)) @@ -76,18 +76,18 @@ function _vertical_advection!( ξ_tend::Grid, # tendency of qu σ̇⁻ = σ_tend_above[ij] # velocity into layer k from above σ̇⁺ = σ_tend_below[ij] # velocity out of layer k to below - ξ⁺ = reconstructed_at_face(ij, adv, σ̇⁺, ξ_sten[2:end]) - ξ⁻ = reconstructed_at_face(ij, adv, σ̇⁻, ξ_sten[1:end-1]) + ξᶠ⁺ = reconstructed_at_face(ij, adv, σ̇⁺, ξ_sten[2:end]) + ξᶠ⁻ = reconstructed_at_face(ij, adv, σ̇⁻, ξ_sten[1:end-1]) - ξ_tend[ij] -= Δσₖ⁻¹ * (σ̇⁺ * ξ⁺ - σ̇⁻ * ξ⁻ - ξ[ij] * (σ̇⁺ - σ̇⁻)) + ξ_tend[ij] -= Δσₖ⁻¹ * (σ̇⁺ * ξᶠ⁺ - σ̇⁻ * ξᶠ⁻ - ξ[ij] * (σ̇⁺ - σ̇⁻)) end end @inline reconstructed_at_face(ij, ::UpwindVerticalAdvection{NF, 1}, u, ξ) where NF = ifelse(u > 0, ξ[1][ij], ξ[2][ij]) -@inline reconstructed_at_face(ij, ::UpwindVerticalAdvection{NF, 2}, u, ξ) where NF = ifelse(u > 0, (2ξ[1][ij] + 5ξ[2][ij] - ξ[3][ij]) / 6, - (-ξ[2][ij] + 5ξ[3][ij] + 2ξ[4][ij]) / 6) -@inline reconstructed_at_face(ij, ::UpwindVerticalAdvection{NF, 3}, u, ξ) where NF = ifelse(u > 0, ( 2ξ[1][ij] - 13ξ[2][ij] + 47ξ[3][ij] + 27ξ[4][ij] - 3ξ[5][ij]) / 60, - (-3ξ[2][ij] + 27ξ[3][ij] + 47ξ[4][ij] - 13ξ[5][ij] + 2ξ[6][ij]) / 60) +@inline reconstructed_at_face(ij, ::UpwindVerticalAdvection{NF, 2}, u, ξ) where NF = ifelse(u > 0, (2ξ[1][ij] + 5ξ[2][ij] - ξ[3][ij]) / 6, + (2ξ[4][ij] + 5ξ[3][ij] - ξ[2][ij]) / 6) +@inline reconstructed_at_face(ij, ::UpwindVerticalAdvection{NF, 3}, u, ξ) where NF = ifelse(u > 0, (2ξ[1][ij] - 13ξ[2][ij] + 47ξ[3][ij] + 27ξ[4][ij] - 3ξ[5][ij]) / 60, + (2ξ[6][ij] - 13ξ[5][ij] + 47ξ[4][ij] + 27ξ[3][ij] - 3ξ[2][ij]) / 60) @inline reconstructed_at_face(ij, ::CenteredVerticalAdvection{NF, 1}, u, ξ) where NF = ( ξ[1][ij] + ξ[2][ij]) / 2 @inline reconstructed_at_face(ij, ::CenteredVerticalAdvection{NF, 2}, u, ξ) where NF = (-ξ[1][ij] + 7ξ[2][ij] + 7ξ[3][ij] - ξ[4][ij]) / 12 From 7945d835b62494c6fb4532fa12056ad36cb55b2a Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Thu, 3 Aug 2023 20:30:28 -0400 Subject: [PATCH 13/18] cleared humid_grid_prev --- src/dynamics/diagnostic_variables.jl | 4 +--- src/dynamics/vertical_advection.jl | 36 ++++++++++++++++++---------- 2 files changed, 24 insertions(+), 16 deletions(-) diff --git a/src/dynamics/diagnostic_variables.jl b/src/dynamics/diagnostic_variables.jl index b18323739..a256c4d64 100644 --- a/src/dynamics/diagnostic_variables.jl +++ b/src/dynamics/diagnostic_variables.jl @@ -49,7 +49,6 @@ struct GridVariables{NF<:AbstractFloat,Grid<:AbstractGrid{NF}} temp_grid_prev ::Grid # absolute temperature of previous time step [K] temp_virt_grid ::Grid # virtual tempereature [K] humid_grid ::Grid # specific_humidity - humid_grid_prev ::Grid # specific_humidity of previous time step geopot_grid ::Grid # geopotential (is that needed?) u_grid ::Grid # zonal velocity *coslat [m/s] v_grid ::Grid # meridional velocity *coslat [m/s] @@ -66,14 +65,13 @@ function Base.zeros(::Type{GridVariables},SG::SpectralGrid) temp_grid_prev = zeros(Grid{NF},nlat_half) # absolute temperature temp_virt_grid = zeros(Grid{NF},nlat_half) # virtual temperature humid_grid = zeros(Grid{NF},nlat_half) # specific humidity - humid_grid_prev = zeros(Grid{NF},nlat_half) # specific humidity geopot_grid = zeros(Grid{NF},nlat_half) # geopotential u_grid = zeros(Grid{NF},nlat_half) # zonal velocity *coslat v_grid = zeros(Grid{NF},nlat_half) # meridonal velocity *coslat u_grid_prev = zeros(Grid{NF},nlat_half) # zonal velocity *coslat v_grid_prev = zeros(Grid{NF},nlat_half) # meridonal velocity *coslat - return GridVariables( vor_grid,div_grid,temp_grid,temp_grid_prev,temp_virt_grid,humid_grid,humid_grid_prev,geopot_grid, + return GridVariables( vor_grid,div_grid,temp_grid,temp_grid_prev,temp_virt_grid,humid_grid,geopot_grid, u_grid,v_grid,u_grid_prev,v_grid_prev) end diff --git a/src/dynamics/vertical_advection.jl b/src/dynamics/vertical_advection.jl index f6fd851b5..8ec81d022 100644 --- a/src/dynamics/vertical_advection.jl +++ b/src/dynamics/vertical_advection.jl @@ -10,17 +10,27 @@ CenteredVerticalAdvection(spectral_grid; order = 2) = CenteredVerticalAdvection{ UpwindVerticalAdvection(spectral_grid; order = 5) = UpwindVerticalAdvection{spectral_grid.NF, (order + 1) ÷ 2}() WENOVerticalAdvection(spectral_grid) = WENOVerticalAdvection{spectral_grid.NF}() -@inline retrieve_time_step(::DiffusiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid_prev)) -@inline retrieve_time_step(::DispersiveVerticalAdvection, variables, var) = getproperty(variables, Symbol(var, :_grid)) +@inline retrieve_previous_time_step(variables, var) = getproperty(variables, Symbol(var, :_grid_prev)) +@inline retrieve_current_time_step(variables, var) = getproperty(variables, Symbol(var, :_grid)) -@inline order(::VerticalAdvection{NF, B}) where {NF, B} = B +@inline retrieve_time_step(::DiffusiveVerticalAdvection, variables, var) = retrieve_previous_time_step(variables, var) +@inline retrieve_time_step(::DispersiveVerticalAdvection, variables, var) = retrieve_current_time_step(variables, var) -@inline function retrieve_stencil(k, layers, var, nlev, vertical_advection::VerticalAdvection{NF, B}) where {NF, B} +@inline function retrieve_current_stencil(k, layers, var, nlev, ::VerticalAdvection{NF, B}) where {NF, B} k_stencil = max.(min.(nlev, k-B:k+B), 1) - ξ_stencil = Tuple(retrieve_time_step(vertical_advection, layers[k].grid_variables, var) for k in k_stencil) + ξ_stencil = Tuple(retrieve_current_time_step(layers[k].grid_variables, var) for k in k_stencil) return ξ_stencil end +@inline function retrieve_previous_stencil(k, layers, var, nlev, ::VerticalAdvection{NF, B}) where {NF, B} + k_stencil = max.(min.(nlev, k-B:k+B), 1) + ξ_stencil = Tuple(retrieve_current_time_step(layers[k].grid_variables, var) for k in k_stencil) + return ξ_stencil +end + +@inline retrieve_stencil(k, layers, var, nlev, scheme::DiffusiveVerticalAdvection) = retrieve_previous_stencil(k, layers, var, nlev, scheme) +@inline retrieve_stencil(k, layers, var, nlev, scheme::DispersiveVerticalAdvection) = retrieve_current_stencil(k, layers, var, nlev, scheme) + function vertical_advection!( layer::DiagnosticVariablesLayer, diagn::DiagnosticVariables, model::PrimitiveEquation) @@ -30,7 +40,7 @@ function vertical_advection!( layer::DiagnosticVariablesLayer, wet_core = model isa PrimitiveWet (; σ_levels_thick, nlev ) = model.geometry - vertical_advection = model.vertical_advection + scheme = model.vertical_advection # for k==1 "above" term is 0, for k==nlev "below" term is zero # avoid out-of-bounds indexing with k_above, k_below as follows @@ -45,18 +55,18 @@ function vertical_advection!( layer::DiagnosticVariablesLayer, for var in (:u, :v, :temp) ξ_tend = getproperty(layer.tendencies, Symbol(var, :_tend_grid)) - ξ_sten = retrieve_stencil(k, diagn.layers, var, nlev, vertical_advection) - ξ = retrieve_time_step(vertical_advection, layer.grid_variables, var) + ξ_sten = retrieve_stencil(k, diagn.layers, var, nlev, scheme) + ξ = retrieve_time_step(scheme, layer.grid_variables, var) - _vertical_advection!(ξ_tend,σ_tend_above,σ_tend_below,ξ_sten,ξ,Δσₖ,vertical_advection) + _vertical_advection!(ξ_tend,σ_tend_above,σ_tend_below,ξ_sten,ξ,Δσₖ,scheme) end if wet_core - ξ_tend = getproperty(layer.tendencies, :humid_tend_grid) - ξ_sten = retrieve_stencil(k, diagn.layers, :humid, nlev, vertical_advection) - ξ = retrieve_time_step(vertical_advection, layer.grid_variables, :humid) + ξ_tend = getproperty(layer.tendencies, :humid_tend_grid) + ξ_sten = retrieve_current_stencil(k, diagn.layers, :humid, nlev, scheme) + ξ = retrieve_current_time_step(layer.grid_variables, :humid) - _vertical_advection!(ξ_tend,σ_tend_above,σ_tend_below,ξ_sten,ξ,Δσₖ,vertical_advection) + _vertical_advection!(ξ_tend,σ_tend_above,σ_tend_below,ξ_sten,ξ,Δσₖ,scheme) end end From 0a796f4b3a3aa34806e55182b6c07027d63c7226 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Mon, 7 Aug 2023 16:31:25 -0400 Subject: [PATCH 14/18] add ZWENO formulation --- src/dynamics/vertical_advection.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/dynamics/vertical_advection.jl b/src/dynamics/vertical_advection.jl index 8ec81d022..a8ee7c767 100644 --- a/src/dynamics/vertical_advection.jl +++ b/src/dynamics/vertical_advection.jl @@ -111,18 +111,20 @@ const d₂ = 1/10 @inline weight_β₁(S, NF) = NF(13/12) * (S[1] - 2S[2] + S[3])^2 + NF(1/4) * ( S[1] - S[3])^2 @inline weight_β₂(S, NF) = NF(13/12) * (S[1] - 2S[2] + S[3])^2 + NF(1/4) * ( S[1] - 4S[2] + 3S[3])^2 -@inline p₀(S) = (2S[1] + 5S[2] - S[3]) / 6 -@inline p₁(S) = (-S[1] + 5S[2] + 2S[3]) / 6 -@inline p₂(S) = (2S[1] - 7S[2] + 11S[3]) / 6 +@inline p₀(S) = (2S[1] + 5S[2] - S[3]) / 6 # downind stencil +@inline p₁(S) = (-S[1] + 5S[2] + 2S[3]) / 6 # upwind stencil +@inline p₂(S) = (2S[1] - 7S[2] + 11S[3]) / 6 # extrapolating stencil + +@inline τ₅(β₀, β₁, β₂) = abs(β₂ - β₀) @inline function weno_reconstruction(S₀, S₁, S₂, NF) β₀ = weight_β₀(S₀, NF) β₁ = weight_β₁(S₁, NF) β₂ = weight_β₂(S₂, NF) - w₀ = NF(d₀) / (β₀ + NF(ε))^2 - w₁ = NF(d₁) / (β₁ + NF(ε))^2 - w₂ = NF(d₂) / (β₂ + NF(ε))^2 + w₀ = NF(d₀) * (1 + (τ₅(β₀, β₁, β₂) / (β₀ + NF(ε)))^2) + w₁ = NF(d₁) * (1 + (τ₅(β₀, β₁, β₂) / (β₁ + NF(ε)))^2) + w₂ = NF(d₂) * (1 + (τ₅(β₀, β₁, β₂) / (β₂ + NF(ε)))^2) w₀, w₁, w₂ = (w₀, w₁, w₂) ./ (w₀ + w₁ + w₂) From 90d72339ceac82d51df6d88fd7282237943417e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Milan=20Kl=C3=B6wer?= Date: Tue, 15 Aug 2023 11:27:28 -0400 Subject: [PATCH 15/18] CenteredVerticalAdvection default --- src/dynamics/models.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dynamics/models.jl b/src/dynamics/models.jl index 41b04807e..d0d86c998 100644 --- a/src/dynamics/models.jl +++ b/src/dynamics/models.jl @@ -139,7 +139,7 @@ Base.@kwdef struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr boundary_layer_drag::BoundaryLayerDrag{NF} = LinearDrag(spectral_grid) temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid) static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid) - vertical_advection::VerticalAdvection{NF} = UpwindVerticalAdvection(spectral_grid) + vertical_advection::VerticalAdvection{NF} = CenteredVerticalAdvection(spectral_grid) # vertical_diffusion::VerticalDiffusion{NF} = VerticalLaplacian(spectral_grid) # NUMERICS @@ -207,7 +207,7 @@ Base.@kwdef struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid) static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid) large_scale_condensation::AbstractCondensation{NF} = SpeedyCondensation(spectral_grid) - vertical_advection::VerticalAdvection{NF} = UpwindVerticalAdvection(spectral_grid) + vertical_advection::VerticalAdvection{NF} = CenteredVerticalAdvection(spectral_grid) # vertical_diffusion::VerticalDiffusion{NF} = VerticalLaplacian(spectral_grid) From 6db585d0161b5b6f632e14e1c19d39f1b4979cba Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Mon, 28 Aug 2023 10:48:17 -0400 Subject: [PATCH 16/18] add some tests? --- test/vertical_advection.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 test/vertical_advection.jl diff --git a/test/vertical_advection.jl b/test/vertical_advection.jl new file mode 100644 index 000000000..acd333bd9 --- /dev/null +++ b/test/vertical_advection.jl @@ -0,0 +1,15 @@ +using SpeedyWeather: WENOVerticalAdvection, CenteredVerticalAdvection, UpwindVerticalAdvection +using SpeedyWeather: vertical_advection! + +@testset "Vertical advection runs" begin + spectral_grid = SpectralGrid() + model_types = (PrimitiveDry, PrimitiveWet) + advection_schems = (WENOVerticalAdvection, CenteredVerticalAdvection, UpwindVerticalAdvection) + + for Model in model_types, VerticalAdvection in advection_schems + model = Model(; spectral_grid, vertical_advection = VerticalAdvection(spectral_grid)) + simulation = initialize!(model) + run!(simulation,n_days=10) + @test simulation.model.feedback.nars_detected == false + end +end \ No newline at end of file From b510f8a8ebc5800a728be0fd83574e7313652310 Mon Sep 17 00:00:00 2001 From: Simone Silvestri <33547697+simone-silvestri@users.noreply.github.com> Date: Mon, 28 Aug 2023 11:02:23 -0400 Subject: [PATCH 17/18] bugfix --- test/vertical_advection.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/vertical_advection.jl b/test/vertical_advection.jl index acd333bd9..b00f699f2 100644 --- a/test/vertical_advection.jl +++ b/test/vertical_advection.jl @@ -3,7 +3,7 @@ using SpeedyWeather: vertical_advection! @testset "Vertical advection runs" begin spectral_grid = SpectralGrid() - model_types = (PrimitiveDry, PrimitiveWet) + model_types = (PrimitiveDryModel, PrimitiveWetModel) advection_schems = (WENOVerticalAdvection, CenteredVerticalAdvection, UpwindVerticalAdvection) for Model in model_types, VerticalAdvection in advection_schems From cc3dfceb733aa6ae180bcdb47b066bf9a064a351 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Milan=20Kl=C3=B6wer?= Date: Mon, 28 Aug 2023 16:02:27 -0400 Subject: [PATCH 18/18] add vertical advection tests to runtests.jl --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 4ddbc9acf..36fcbb080 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,6 +18,7 @@ include("spectral_gradients.jl") # DYNAMICS include("diffusion.jl") include("time_stepping.jl") +include("vertical_advection.jl") # VERTICAL LEVELS include("vertical_levels.jl")