diff --git a/docs/src/API.md b/docs/src/API.md index 0a6914445..eb305e226 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -67,7 +67,7 @@ Microphysics2M ```@autodocs Modules = [Parameters] -Pages = ["Microphysics2M.jl", "Microphysics2MParams.jl"] +Pages = ["Microphysics2M.jl", "Microphysics2MParams.jl", "ActivationSchemes.jl"] ``` ## Size distributions @@ -136,6 +136,7 @@ BulkMicrophysicsTendencies.Instantaneous BulkMicrophysicsTendencies.InstantaneousVerbose BulkMicrophysicsTendencies.LinearizedAverage BulkMicrophysicsTendencies.bulk_microphysics_tendencies +BulkMicrophysicsTendencies.activation_source ``` # P3 scheme diff --git a/src/BulkMicrophysicsTendencies.jl b/src/BulkMicrophysicsTendencies.jl index 86265d455..f35a39a2c 100644 --- a/src/BulkMicrophysicsTendencies.jl +++ b/src/BulkMicrophysicsTendencies.jl @@ -663,6 +663,38 @@ end # --- 2-Moment Microphysics Helper Functions --- +""" + activation_source(scheme, tps, ρ, T, q_tot, q_lcl, q_ice, n_lcl, w, p) + +Return the cloud-droplet number source `dn_lcl_activation_dt` [kg⁻¹ air / s] +implied by the activation `scheme` at the given local state. + +All arguments are passed by value so the function is pure and suitable for +both broadcast and GPU kernels; inputs not needed by a particular scheme +(e.g. `w`, `p` for [`CMP.DiagnosticNc`](@ref)) are simply ignored. + +# Arguments +- `scheme`: activation scheme; dispatch target. +- `tps`: thermodynamics parameters. +- `ρ`, `T`, `q_tot`, `q_lcl`, `q_ice`: local thermodynamic state. +- `n_lcl`: cloud-droplet number per mass of air [kg⁻¹]. +- `w`: vertical velocity [m/s] (for supersaturation-driven schemes). +- `p`: pressure [Pa] (for supersaturation-driven schemes). +""" +@inline function activation_source( + ::CMP.NoActivation, tps, ρ, T, q_tot, q_lcl, q_ice, n_lcl, w, p, +) + return zero(n_lcl) +end + +@inline function activation_source( + scheme::CMP.DiagnosticNc, tps, ρ, T, q_tot, q_lcl, q_ice, n_lcl, w, p, +) + FT = UT.promote_typeof(n_lcl, q_lcl) + target = ifelse(q_lcl > FT(scheme.q_thresh), FT(scheme.N_c), zero(FT)) + return (target - n_lcl) / FT(scheme.τ_relax) +end + """ warm_rain_tendencies_2m(sb, q_lcl, q_rai, ρ, n_lcl, n_rai) @@ -708,7 +740,13 @@ Used by both warm-only and warm+ice dispatch methods to reduce code duplication. dn_rai_dt = zero(FT) # --- Aerosol activation --- - dn_lcl_activation_dt = zero(FT) + # `w`, `p` are the per-cell ambient inputs that supersaturation-driven + # schemes will need; positional so the call site can `@.`-broadcast. + # `NoActivation`/`DiagnosticNc` ignore them. + dn_lcl_activation_dt = activation_source( + warm_rain.activation_scheme, tps, ρ, T, q_tot, q_lcl, q_ice, n_lcl, w, p, + ) + dn_lcl_dt += dn_lcl_activation_dt # --- Condensation of vapor / evaporation of cloud liquid water --- micro_mock = (; q_tot, q_lcl, q_icl = q_ice, q_rai, q_sno = zero(q_ice)) diff --git a/src/parameters/ActivationSchemes.jl b/src/parameters/ActivationSchemes.jl new file mode 100644 index 000000000..98797e1fc --- /dev/null +++ b/src/parameters/ActivationSchemes.jl @@ -0,0 +1,63 @@ +export AbstractActivationScheme, NoActivation, DiagnosticNc + +""" + AbstractActivationScheme + +Abstract super-type for aerosol-activation closures used inside +`BulkMicrophysicsTendencies`. + +Every concrete subtype is a pure, immutable parameter struct that fully +specifies one activation closure. BMT dispatches on the scheme type to +compute the cloud-droplet number source `dn_lcl_activation_dt` [kg⁻¹ s⁻¹] +given the local thermodynamic state — the same options-pattern the 1-moment +scheme uses for its prescribed-droplet-number autoconversion +([`PrescribedNd`](@ref)). + +Currently provided: [`NoActivation`](@ref) (the null source, default) and +[`DiagnosticNc`](@ref) (relaxation toward a prescribed droplet number). +Supersaturation-driven closures (Twomey, Abdul-Razzak–Ghan) are deferred to +a broader activation API. +""" +abstract type AbstractActivationScheme <: ParametersType end + +""" + NoActivation() + +Singleton scheme that suppresses aerosol activation: the activation source +is zero regardless of state. The default — used when the host model supplies +its own activation, and for tests that isolate other processes. +""" +struct NoActivation <: AbstractActivationScheme end + +""" + DiagnosticNc(; N_c, q_thresh = 1e-7, τ_relax = 60) + +Diagnostic-`N_c` activation: relax the droplet number toward a prescribed +target. The RCEMIP-I default and the pre-ARG GCM standard. + +- While `q_lcl > q_thresh`, the droplet target is `N_c` [kg⁻¹ air] and the + tendency is `(N_c − n_lcl) / τ_relax`. +- Otherwise the target is zero (cloud has evaporated) and the tendency is + `-n_lcl / τ_relax`. + +# Design notes +- `τ_relax` is a FIXED physical timescale, NOT the integrator timestep: + relaxation by `1/dt` produces an infinitely-stiff source as adaptive + solvers refine the step. The 60 s default keeps `n_lcl` tracking cloud + mass within one macro-step while remaining integrable. +- `q_thresh` avoids activating into essentially vapor-only cells, where + numerical noise in `q_lcl` would otherwise create phantom droplets. + +# Fields +$(DocStringExtensions.FIELDS) +""" +@kwdef struct DiagnosticNc{FT} <: AbstractActivationScheme + "Droplet number target [kg⁻¹ air] while cloud mass is present" + N_c::FT + "Cloud-existence threshold on `q_lcl` [kg/kg]. Default: 1e-7." + q_thresh::FT = oftype(N_c, 1e-7) + "Fixed relaxation timescale [s]. Default: 60." + τ_relax::FT = oftype(N_c, 60) +end + +Base.broadcastable(x::AbstractActivationScheme) = tuple(x) diff --git a/src/parameters/Microphysics2MParams.jl b/src/parameters/Microphysics2MParams.jl index a8e059270..1fcc298d3 100644 --- a/src/parameters/Microphysics2MParams.jl +++ b/src/parameters/Microphysics2MParams.jl @@ -10,20 +10,27 @@ Parameters for 2-moment warm rain processes (Seifert-Beheng 2006). - `air_properties::AP`: AirProperties — air properties for evaporation - `condevap::CE`: MM2015 cond-evap relaxation timescale - `subdep::SD`: MM2015 sub-dep relaxation timescale +- `activation_scheme::AS`: aerosol → cloud-droplet activation closure + (defaults to [`NoActivation`](@ref); set to e.g. [`DiagnosticNc`](@ref) + to enable activation in the fused 2M tendency). """ -@kwdef struct WarmRainParams2M{SB, AP, CE, SD} <: ParametersType +@kwdef struct WarmRainParams2M{SB, AP, CE, SD, AS <: AbstractActivationScheme} <: ParametersType seifert_beheng::SB air_properties::AP condevap::CE subdep::SD + activation_scheme::AS = NoActivation() end # Construct WarmRainParams2M from a ClimaParams TOML dictionary -WarmRainParams2M(toml_dict::CP.ParamDict; is_limited = true) = +WarmRainParams2M(toml_dict::CP.ParamDict; is_limited = true, + activation_scheme::AbstractActivationScheme = NoActivation(), +) = WarmRainParams2M(; seifert_beheng = SB2006(toml_dict; is_limited), air_properties = AirProperties(toml_dict), condevap = CondEvap2M(toml_dict), subdep = SubDep2M(toml_dict), + activation_scheme, ) Base.show(io::IO, mime::MIME"text/plain", x::WarmRainParams2M) = @@ -146,10 +153,11 @@ Create a `Microphysics2MParams` object from a ClimaParams TOML dictionary. Microphysics2MParams(toml_dict::CP.ParamDict; with_ice = false, is_limited = true, quadrature = QUAD.ChebyshevGauss(100), + activation_scheme = NoActivation(), inp_depletion_model = NIceProxyDepletion(τ_act = 300), ) = Microphysics2MParams(; # Warm rain parameters (always present) - warm_rain = WarmRainParams2M(toml_dict; is_limited), + warm_rain = WarmRainParams2M(toml_dict; is_limited, activation_scheme), # Optional ice phase parameters ice = with_ice ? P3IceParams(toml_dict; is_limited, quadrature, inp_depletion_model) : diff --git a/src/parameters/Parameters.jl b/src/parameters/Parameters.jl index 3dd738a25..e6bc98c37 100644 --- a/src/parameters/Parameters.jl +++ b/src/parameters/Parameters.jl @@ -49,6 +49,7 @@ include("MicrophysicsP3.jl") include("TerminalVelocity.jl") # Unified parameter containers +include("ActivationSchemes.jl") include("Microphysics0MParams.jl") include("Microphysics1MOptions.jl") include("Microphysics1MParams.jl") diff --git a/test/activation_schemes_tests.jl b/test/activation_schemes_tests.jl new file mode 100644 index 000000000..dde6451d8 --- /dev/null +++ b/test/activation_schemes_tests.jl @@ -0,0 +1,74 @@ +using Test + +import ClimaParams as CP +import CloudMicrophysics as CM +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.BulkMicrophysicsTendencies as BMT +import CloudMicrophysics.ThermodynamicsInterface as TDI +import ForwardDiff as FD + +function test_activation_schemes(FT) + tps = TDI.TD.Parameters.ThermodynamicsParameters(FT) + args = (tps, FT(1.05), FT(288), FT(0.012), FT(2e-3), FT(0), FT(5e7), FT(0), FT(0)) + splice(a; q_lcl = a[5], n_lcl = a[7]) = (a[1:4]..., q_lcl, a[6], n_lcl, a[8:9]...) + + @testset "NoActivation is the null source ($FT)" begin + @test BMT.activation_source(CMP.NoActivation(), args...) === zero(FT) + end + + @testset "DiagnosticNc relaxation ($FT)" begin + scheme = CMP.DiagnosticNc(; N_c = FT(1e8)) + # cloudy: relax up toward the target + s = BMT.activation_source(scheme, args...) + @test s ≈ (FT(1e8) - FT(5e7)) / FT(60) + # above target: relax down + s_dn = BMT.activation_source(scheme, splice(args; n_lcl = FT(2e8))...) + @test s_dn ≈ (FT(1e8) - FT(2e8)) / FT(60) + # cloud-free: target zero (decay of leftover droplets) + s0 = BMT.activation_source(scheme, splice(args; q_lcl = FT(0))...) + @test s0 ≈ -FT(5e7) / FT(60) + # custom timescale and threshold + s2 = BMT.activation_source( + CMP.DiagnosticNc(; N_c = FT(1e8), q_thresh = FT(1e-2), τ_relax = FT(30)), args..., + ) + @test s2 ≈ -FT(5e7) / FT(30) # q_lcl below the raised threshold + end + + @testset "activation in the fused 2M tendency ($FT)" begin + mp_on = CMP.Microphysics2MParams( + FT; activation_scheme = CMP.DiagnosticNc(; N_c = FT(1e8)), + ) + mp_off = CMP.Microphysics2MParams(FT) + @test mp_off.warm_rain.activation_scheme isa CMP.NoActivation + t_on = BMT.bulk_microphysics_tendencies( + BMT.Microphysics2Moment(), mp_on, tps, + FT(1.05), FT(288), FT(0.012), FT(2e-3), FT(5e7), FT(1e-4), FT(4e4), + ) + t_off = BMT.bulk_microphysics_tendencies( + BMT.Microphysics2Moment(), mp_off, tps, + FT(1.05), FT(288), FT(0.012), FT(2e-3), FT(5e7), FT(1e-4), FT(4e4), + ) + @test t_off.dn_lcl_activation_dt === zero(FT) + @test t_on.dn_lcl_activation_dt ≈ (FT(1e8) - FT(5e7)) / FT(60) + # the activation source is folded into the total droplet tendency + @test t_on.dn_lcl_dt - t_off.dn_lcl_dt ≈ t_on.dn_lcl_activation_dt + end + + @testset "activation source is concretely typed under mixed args ($FT)" begin + D = FD.Dual{Nothing, FT, 8} + scheme = CMP.DiagnosticNc(; N_c = FT(1e8)) + for sig in ( + (FT, FT, FT, D, FT, D, FT, FT), # Dual q_lcl, n_lcl + (FT, FT, FT, FT, FT, D, FT, FT), # Dual n_lcl only + (FT, FT, FT, D, FT, FT, FT, FT), # Dual q_lcl only + ) + rts = Base.return_types( + BMT.activation_source, Tuple{typeof(scheme), typeof(tps), FT, sig...}, + ) + @test all(isconcretetype, rts) + end + end +end + +test_activation_schemes(Float64) +test_activation_schemes(Float32) diff --git a/test/runtests.jl b/test/runtests.jl index 20c6c78de..7839e4599 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -27,6 +27,7 @@ TT.@testset "All tests" begin include("ventilation_tests.jl") include("ad_compat_tests.jl") include("return_type_tests.jl") + include("activation_schemes_tests.jl") include("DistributionTools_tests.jl") include("unrolled_logsumexp.jl") end