From aa279a02c78f7ef935c3b9685723802c5f4b600d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Milan=20Kl=C3=B6wer?= Date: Fri, 6 Dec 2024 13:59:53 +0000 Subject: [PATCH] AbstractSurfacePerturbation (#631) --- CHANGELOG.md | 1 + docs/src/stochastic_physics.md | 6 +- src/physics/convection.jl | 92 ++++++++++++++++++++++++++----- src/physics/longwave_radiation.jl | 7 ++- test/runtests.jl | 1 + test/stochastic_physics.jl | 76 +++++++++++++++++++++++++ 6 files changed, 164 insertions(+), 19 deletions(-) create mode 100644 test/stochastic_physics.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index c1d83c9e6..3fe7a32c1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Unreleased +- AbstractSurfacePerturbation introduced [#631](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/631) - AbstractStochasticPhysics and SPPT implemented [#629](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/629) - Introduced seperate extended tests that are not run on every commit [#628](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/628) - One-band longwave radiation [#624](https://github.com/SpeedyWeather/SpeedyWeather.jl/pull/624) diff --git a/docs/src/stochastic_physics.md b/docs/src/stochastic_physics.md index 468b58205..5c128128c 100644 --- a/docs/src/stochastic_physics.md +++ b/docs/src/stochastic_physics.md @@ -3,9 +3,9 @@ Stochastic physics introduces stochasticity into the parameterizations of physical process. There is conceptually several classes of how this can be done - - Stochastic perturbations of the tendencies - - Stochastic parameter perturbations - - Stochastic perturbations to the inputs of parameterizations +- Stochastic perturbations of the tendencies +- Stochastic parameter perturbations +- Stochastic perturbations to the inputs of parameterizations all of these use random numbers created from some random processes (white noise or with some autocorrelation in space or time or sampled from some other distribution) diff --git a/src/physics/convection.jl b/src/physics/convection.jl index 758aacf95..44c57421b 100644 --- a/src/physics/convection.jl +++ b/src/physics/convection.jl @@ -1,6 +1,39 @@ +abstract type AbstractSurfacePerturbation <: AbstractParameterization end + +# subtypes don't require an initialize! method defined by default +initialize!(::AbstractSurfacePerturbation, ::PrimitiveEquation) = nothing + +export NoSurfacePerturbation + +"""Returns the surface temperature and humidity without +perturbation from the lowermost layer.""" +struct NoSurfacePerturbation <: AbstractSurfacePerturbation end + +"""$(TYPEDSIGNATURES) +Returns the surface temperature and humidity without +perturbation from the lowermost layer. Used as a functor.""" +function (SP::NoSurfacePerturbation)( + column::ColumnVariables, + model::PrimitiveEquation, +) + (; temp, humid, nlayers) = column + return temp[nlayers], humid[nlayers] +end + +# only return temperature for primitive dry model +function (SP::NoSurfacePerturbation)( + column::ColumnVariables, + model::PrimitiveDry, +) + (; temp, nlayers) = column + return temp[nlayers] +end + abstract type AbstractConvection <: AbstractParameterization end export NoConvection + +"""Dummy type to disable convection.""" struct NoConvection <: AbstractConvection end NoConvection(::SpectralGrid) = NoConvection() initialize!(::NoConvection, ::PrimitiveEquation) = nothing @@ -8,12 +41,11 @@ convection!(::ColumnVariables, ::NoConvection, ::PrimitiveEquation) = nothing export SimplifiedBettsMiller -""" -The simplified Betts-Miller convection scheme from Frierson, 2007, +"""The simplified Betts-Miller convection scheme from Frierson, 2007, https://doi.org/10.1175/JAS3935.1. This implements the qref-formulation in their paper. Fields and options are $(TYPEDFIELDS)""" -@kwdef struct SimplifiedBettsMiller{NF} <: AbstractConvection +@kwdef struct SimplifiedBettsMiller{NF, SP} <: AbstractConvection "number of vertical layers" nlayers::Int @@ -22,10 +54,23 @@ $(TYPEDFIELDS)""" "[OPTION] Relative humidity for reference profile" relative_humidity::NF = 0.7 + + "[OPTION] Surface perturbation of temp, humid to calculate the moist pseudo adiabat" + surface_temp_humid::SP # don't specify default here but in generator below +end + +# generator function +function SimplifiedBettsMiller( + SG::SpectralGrid; + surface_temp_humid = NoSurfacePerturbation(), + kwargs...) + # type inference happens here as @kwdef only defines no/all types specified + return SimplifiedBettsMiller{SG.NF, typeof(surface_temp_humid)}( + nlayers=SG.nlayers; surface_temp_humid, kwargs...) end -SimplifiedBettsMiller(SG::SpectralGrid; kwargs...) = SimplifiedBettsMiller{SG.NF}(nlayers=SG.nlayers; kwargs...) -initialize!(::SimplifiedBettsMiller, ::PrimitiveEquation) = nothing +# nothing to initialize but maybe the surface temp/humid functor? +initialize!(SBM::SimplifiedBettsMiller, model::PrimitiveEquation) = initialize!(SBM.surface_temp_humid, model) # function barrier for all AbstractConvection function convection!( @@ -42,7 +87,7 @@ function convection!( model::PrimitiveWet, ) convection!(column, scheme, model.clausius_clapeyron, - model.geometry, model.planet, model.atmosphere, model.time_stepping) + model.geometry, model.planet, model.atmosphere, model.time_stepping, model) end """ @@ -60,6 +105,7 @@ function convection!( planet::AbstractPlanet, atmosphere::AbstractAtmosphere, time_stepping::AbstractTimeStepper, + model::PrimitiveWet, ) where NF σ = geometry.σ_levels_full @@ -74,8 +120,8 @@ function convection!( humid_ref_profile = column.b # specific humidity [kg/kg] profile to adjust to # CONVECTIVE CRITERIA AND FIRST GUESS RELAXATION - temp_parcel = temp[nlayers] - humid_parcel = humid[nlayers] + # force conversion to NF here + temp_parcel, humid_parcel = convert.(NF, SBM.surface_temp_humid(column, model)) level_zero_buoyancy = pseudo_adiabat!(temp_ref_profile, temp_parcel, humid_parcel, temp_virt, geopot, pₛ, σ, @@ -248,16 +294,28 @@ The simplified Betts-Miller convection scheme from Frierson, 2007, https://doi.org/10.1175/JAS3935.1 but with humidity set to zero. Fields and options are $(TYPEDFIELDS)""" -@kwdef struct DryBettsMiller{NF} <: AbstractConvection +@kwdef struct DryBettsMiller{NF, SP} <: AbstractConvection "number of vertical layers/levels" nlayers::Int "[OPTION] Relaxation time for profile adjustment" time_scale::Second = Hour(4) + + "[OPTION] Surface perturbation of temp to calculate the dry adiabat" + surface_temp::SP # don't set default here but in generator below end -DryBettsMiller(SG::SpectralGrid; kwargs...) = DryBettsMiller{SG.NF}(nlayers=SG.nlayers; kwargs...) -initialize!(::DryBettsMiller, ::PrimitiveEquation) = nothing +# generator function +function DryBettsMiller( + SG::SpectralGrid; + surface_temp = NoSurfacePerturbation(), + kwargs...) + # infer type here as @kwdef only defines constructors for no/all types specified + return DryBettsMiller{SG.NF, typeof(surface_temp)}(nlayers=SG.nlayers; surface_temp, kwargs...) +end + +# nothing to initialize but maybe the surface temp functor? +initialize!(DBM::DryBettsMiller, model::PrimitiveEquation) = initialize!(DBM.surface_temp, model) # function barrier to unpack model function convection!( @@ -265,7 +323,8 @@ function convection!( scheme::DryBettsMiller, model::PrimitiveEquation, ) - convection!(column, scheme, model.geometry, model.atmosphere) + # TODO check whether unpacking makes a performance difference? (it shouldn't?) + convection!(column, scheme, model.geometry, model.atmosphere, model) end """ @@ -281,6 +340,7 @@ function convection!( DBM::DryBettsMiller, geometry::Geometry, atmosphere::AbstractAtmosphere, + model::PrimitiveEquation, ) where NF σ = geometry.σ_levels_full @@ -292,8 +352,12 @@ function convection!( temp_ref_profile = column.a # temperature [K] reference profile to adjust to # CONVECTIVE CRITERIA AND FIRST GUESS RELAXATION + # force conversion to NF here, add "," to ignore additional returns like humidity + temp_parcel, = convert.(NF, DBM.surface_temp(column, model)) level_zero_buoyancy = dry_adiabat!(temp_ref_profile, - temp, σ, + temp, + temp_parcel, + σ, atmosphere) local PT::NF = 0 # precipitation due to coolinga @@ -331,6 +395,7 @@ set to NaN instead and should be skipped in the relaxation.""" function dry_adiabat!( temp_ref_profile::AbstractVector, temp_environment::AbstractVector, + temp_parcel::Real, σ::AbstractVector, atmosphere::AbstractAtmosphere, ) @@ -342,7 +407,6 @@ function dry_adiabat!( length(σ) == length(temp_environment) || throw(BoundsError) nlayers = length(temp_ref_profile) # number of vertical levels - temp_parcel = temp_environment[nlayers] # parcel is at lowermost layer temperature temp_ref_profile .= NaN # reset profile from any previous calculation temp_ref_profile[nlayers] = temp_parcel # start profile with parcel temperature diff --git a/src/physics/longwave_radiation.jl b/src/physics/longwave_radiation.jl index a5079d84e..e888b4e5a 100644 --- a/src/physics/longwave_radiation.jl +++ b/src/physics/longwave_radiation.jl @@ -16,7 +16,7 @@ end export UniformCooling """ -Uniform cooling following Pauluis and Garner, 2006. JAS. https://doi.org/10.1175/JAS3705.1 +Uniform cooling following Paulius and Garner, 2006. JAS. https://doi.org/10.1175/JAS3705.1 imposing a default temperature tendency of -1.5K/day (=1K/16hours for a `time_scale` of 16 hours) on every level except for the stratosphere (diagnosed as `temp < temp_min`) where a relaxation term with `time_scale_stratosphere` towards `temp_stratosphere` is applied. @@ -199,9 +199,12 @@ function longwave_radiation!( column.flux_temp_upward[k] += U # accumulate that flux end + # store outgoing longwave radiation (OLR) for diagnostics + column.outgoing_longwave_radiation = U + # DOWNWARD flux D local D::NF = 0 # top boundary condition of longwave flux - # non need to accumulate 0 at top downward flux + # no need to accumulate 0 at top downward flux @inbounds for k in 1:nlayers # Frierson et al. 2006, equation 7 D += dτ[k]*(B[k] - D) diff --git a/test/runtests.jl b/test/runtests.jl index 6b2bec67d..b590bc4e0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -48,6 +48,7 @@ include("optical_depth.jl") include("longwave_radiation.jl") # include("shortwave_radiation.jl") include("random_process.jl") +include("stochastic_physics.jl") # INITIALIZATION AND INTEGRATION include("run_speedy.jl") diff --git a/test/stochastic_physics.jl b/test/stochastic_physics.jl new file mode 100644 index 000000000..899a4abb1 --- /dev/null +++ b/test/stochastic_physics.jl @@ -0,0 +1,76 @@ +@testset "Surface perturbation wet/dry model" begin + @kwdef struct MySurfacePerturbation <: SpeedyWeather.AbstractSurfacePerturbation + noise_scale_temp::Float32 = 1 # standard deviation of perturbation, in Kelvin + end + + # implement the actual perturbation as a functor + function (SP::MySurfacePerturbation)( + column::ColumnVariables, + model::PrimitiveEquation, + ) + (; temp, humid, nlayers) = column + T, q = temp[nlayers], humid[nlayers] # unperturbed lowermost layer variables + + # perturb temperature additive only + r = column.random_value + T += SP.noise_scale_temp*r + + return T, q + end + + spectral_grid = SpectralGrid() + my_surface_perturbation = MySurfacePerturbation(noise_scale_temp=1) + + # construct convection with that perturbation + convection = SimplifiedBettsMiller(spectral_grid, surface_temp_humid=my_surface_perturbation) + + # don't forget the random process otherwise your `column.random_value` is always 0! + random_process = SpectralAR1Process(spectral_grid, seed=1) + + # construct model + model = PrimitiveWetModel(spectral_grid; random_process, convection) + simulation = initialize!(model) + + run!(simulation, period=Day(5)) + + k = spectral_grid.nlayers + surface_humid = simulation.diagnostic_variables.grid.humid_grid[:, k] + + # reconstruct model with different seed + random_process = SpectralAR1Process(spectral_grid, seed=2) + model = PrimitiveWetModel(spectral_grid; random_process, convection) + # reinitialize, rerun check whether it's different! + simulation = initialize!(model) # choose a different random seed for random_process + run!(simulation, period=Day(5)) + surface_humid2 = simulation.diagnostic_variables.grid.humid_grid[:, k] + + for ij in eachindex(surface_humid, surface_humid2) + if surface_humid[ij] != 0 + @test surface_humid[ij] != surface_humid2[ij] + end + end + + # for the dry model too + convection = DryBettsMiller(spectral_grid, surface_temp=my_surface_perturbation) + model = PrimitiveDryModel(spectral_grid; random_process, convection) + simulation = initialize!(model) + run!(simulation, period=Day(5)) + + k = spectral_grid.nlayers + surface_temp = simulation.diagnostic_variables.grid.temp_grid[:, k] + + # reinitialize, rerun check whether it's different! + simulation = initialize!(model) # choose a different random seed for random_process + run!(simulation, period=Day(5)) + surface_temp2 = simulation.diagnostic_variables.grid.temp_grid[:, k] + + # count the number of grid points that are different in temperature + N = length(surface_temp) + local n::Int64 = 0 + for ij in eachindex(surface_temp, surface_temp2) + n += surface_temp[ij] != surface_temp2[ij] + end + + # 75% of grid points must be different + @test n/N > 3/4 +end \ No newline at end of file