From 09d92b0d9a877774aba16b82246f10b9a782f25a Mon Sep 17 00:00:00 2001 From: Haakon Ludvig Langeland Ervik <45243236+haakon-e@users.noreply.github.com> Date: Mon, 15 Jun 2026 17:04:15 -0400 Subject: [PATCH] rft: cleanup from #681 - rename `f23_deposition_rate` -> `deposition_rate` - rename `f23_immersion_limit_rate` -> `immersion_limit_rate` - delete `AbstractINPDepletion` abstract type - removed nonsensical comments - trimmed verbose "docs-like" code comments. Follow-up: proper user/developer docs - removed references to internal notes - added TODOs for follow-up work - clean up some docstrings - clean up collision closed form code --- docs/src/API.md | 4 +- src/BulkMicrophysicsTendencies.jl | 150 +++-------- src/IceNucleation.jl | 120 ++++----- src/Microphysics2M.jl | 12 - src/P3_particle_properties.jl | 82 ++---- src/P3_processes.jl | 277 ++++++--------------- src/P3_size_distribution.jl | 111 +-------- src/P3_terminal_velocity.jl | 19 +- src/parameters/IceNucleation.jl | 32 +-- test/gpu_tests.jl | 18 +- test/heterogeneous_ice_nucleation_tests.jl | 26 +- test/p3_tests.jl | 235 +++++------------ test/performance_tests.jl | 6 +- 13 files changed, 275 insertions(+), 817 deletions(-) diff --git a/docs/src/API.md b/docs/src/API.md index c8ca29251..349326449 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -333,8 +333,8 @@ HetIceNucleation.P3_deposition_N_i HetIceNucleation.P3_het_N_i HetIceNucleation.INP_concentration_frequency HetIceNucleation.INP_concentration_mean -HetIceNucleation.f23_deposition_rate -HetIceNucleation.f23_immersion_limit_rate +HetIceNucleation.deposition_rate +HetIceNucleation.immersion_limit_rate HetIceNucleation.liquid_freezing_rate HetIceNucleation.n_active ``` diff --git a/src/BulkMicrophysicsTendencies.jl b/src/BulkMicrophysicsTendencies.jl index 8a8abebac..958887851 100644 --- a/src/BulkMicrophysicsTendencies.jl +++ b/src/BulkMicrophysicsTendencies.jl @@ -666,7 +666,7 @@ end """ warm_rain_tendencies_2m(sb, q_lcl, q_rai, ρ, n_lcl, n_rai) -Internal helper function that computes core 2M warm rain processes: +Internal helper function that computes 2M warm rain processes: autoconversion, self-collection, accretion, and rain breakup. Used by both warm-only and warm+ice dispatch methods to reduce code duplication. @@ -680,7 +680,7 @@ Used by both warm-only and warm+ice dispatch methods to reduce code duplication. - `n_rai`: Rain number per kg air (1/kg) # Returns -`NamedTuple` with core warm rain tendencies: +`NamedTuple` with warm rain tendencies: - `dq_lcl_dt`: Cloud liquid mass tendency (kg/kg/s) - `dq_rai_dt`: Rain mass tendency (kg/kg/s) - `dn_lcl_dt`: Cloud number tendency (1/kg/s) @@ -708,10 +708,6 @@ Used by both warm-only and warm+ice dispatch methods to reduce code duplication. dn_rai_dt = zero(FT) # --- Aerosol activation --- - # No activation source: the active schemes (e.g. diagnostic-Nc) - # will appear in a follow-up PR. - # The `w`, `p` ambient inputs (which the active schemes needed) are kept as - # trailing defaults for that re-add. dn_lcl_activation_dt = zero(FT) # --- Condensation of vapor / evaporation of cloud liquid water --- @@ -747,12 +743,12 @@ Used by both warm-only and warm+ice dispatch methods to reduce code duplication. dn_lcl_dt += accr.dN_lcl_dt / ρ # --- Rain self-collection --- - ∂ₜN_rai_sc = CM2.rain_self_collection(sb.pdf_r, sb.self, q_rai, ρ, N_rai) - dn_rai_dt += ∂ₜN_rai_sc / ρ + ∂ρn_rai_sc_∂t = CM2.rain_self_collection(sb.pdf_r, sb.self, q_rai, ρ, N_rai) + dn_rai_dt += ∂ρn_rai_sc_∂t / ρ # --- Rain breakup --- - ∂ₜN_rai_br = CM2.rain_breakup(sb.pdf_r, sb.brek, q_rai, ρ, N_rai, ∂ₜN_rai_sc) - dn_rai_dt += ∂ₜN_rai_br / ρ + ∂ρn_rai_br_∂t = CM2.rain_breakup(sb.pdf_r, sb.brek, q_rai, ρ, N_rai, ∂ρn_rai_sc_∂t) + dn_rai_dt += ∂ρn_rai_br_∂t / ρ # --- Number adjustment for mass limits --- # Cloud liquid @@ -803,7 +799,7 @@ For warm rain + P3 ice, see the method that accepts `Microphysics2MParams{FT, WR - `dq_rim_dt`: Rime mass tendency (always zero for warm-only) - `db_rim_dt`: Rime volume tendency (always zero for warm-only) """ -@inline function bulk_microphysics_tendencies( +@inline function bulk_microphysics_tendencies( # TODO: Delete this function ::Microphysics2Moment, mp::CMP.Microphysics2MParams{WR, Nothing}, tps, ρ, T, q_tot, q_lcl, n_lcl, q_rai, n_rai, q_ice = zero(ρ), n_ice = zero(ρ), q_rim = zero(ρ), b_rim = zero(ρ), logλ = zero(ρ), @@ -827,7 +823,7 @@ For warm rain + P3 ice, see the method that accepts `Microphysics2MParams{FT, WR dq_rim_dt = zero(ρ) db_rim_dt = zero(ρ) - # --- Core Warm Rain Processes (shared helper, includes activation) --- + # --- Warm Rain Processes warm = warm_rain_tendencies_2m(mp.warm_rain, tps, T, q_tot, q_lcl, q_rai, q_ice, ρ, n_lcl, n_rai, w, p) dq_lcl_dt = warm.dq_lcl_dt dn_lcl_dt = warm.dn_lcl_dt @@ -914,12 +910,6 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch. N_ice = n_ice * ρ # [1 / m³ air] L_rim = q_rim * ρ # [kg rim / m³ air] B_rim = b_rim * ρ # [m³ rim / m³ air] - - # `state_from_prognostic` handles the F_rim/ρ_rim regularisation - # and domain clamps internally. The regularised ratios - # (`UT.rime_mass_fraction`, `UT.rime_density`) smoothly blend to - # zero as their denominators shrink, avoiding the hard - # discontinuity at `q_ice = ϵ` / `b_rim = ϵ`. state = CMP3.state_from_prognostic(mp.ice.scheme, L_ice, N_ice, L_rim, B_rim) # Unpack warm rain parameters @@ -932,45 +922,27 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch. dq_rim_dt = zero(ρ) db_rim_dt = zero(ρ) - # --- Core Warm Rain Processes (includes activation, see warm_rain_tendencies_2m) --- + # --- Warm Rain Processes warm = warm_rain_tendencies_2m(mp.warm_rain, tps, T, q_tot, q_lcl, q_rai, q_ice, ρ, n_lcl, n_rai, w, p) dq_lcl_dt = warm.dq_lcl_dt dn_lcl_dt = warm.dn_lcl_dt dq_rai_dt = warm.dq_rai_dt dn_rai_dt = warm.dn_rai_dt dn_lcl_activation_dt = warm.dn_lcl_activation_dt - # NOTE on latent-heat coupling: per-process phase-change rates - # (`S_cond`, `S_dep`, `S_frz_net`) are no longer returned. Hosts that - # need a combined LH rate can compute it from the prognostic mass - # tendencies via the identity - # ḣ_lh = Lv(T) · (dq_lcl/dt + dq_rai/dt) + Ls(T) · dq_ice/dt - # which holds because Ls = Lv + Lf (sublimation = vaporisation + - # fusion) makes the L_f · S_frz term cancel when expressed in - # species coordinates. This keeps BMT dt-agnostic and the host's LH - # rate self-consistent with whatever clipping the host applies to - # the species tendencies. - - # --- P3 Ice Processes --- + + # --- P3 Ice Processes p3 = mp.ice.scheme vel = mp.ice.terminal_velocity pdf_c = mp.ice.cloud_pdf pdf_r = mp.ice.rain_pdf ice_nucleation = mp.ice.ice_nucleation inp_depletion_model = mp.ice.inp_depletion_model - - # Quadrature for the P3 size-distribution integrals. The rule is built - # once (host-side, from `quadrature_order`) and stored on `P3IceParams`, - # so here it is just read — no per-cell construction inside this (GPU) - # kernel. GaussLegendre is selected for the orders where it is - # meaningfully more accurate than ChebyshevGauss on the smooth integrands - # (≈20× lower error on the dominant ice-rain collision integral; see - # `Quadrature.build_quadrature` / `GaussLegendre`), otherwise ChebyshevGauss. quad = mp.ice.quad # Only compute ice processes if there is ice mass/number present if q_ice > ϵₘ && n_ice > ϵₙ - # --- Liquid-ice collisions --- + # --- Liquid-ice collisions coll = CMP3.bulk_liquid_ice_collision_sources( state, logλ, pdf_c, pdf_r, L_lcl, N_lcl, L_rai, N_rai, aps, tps, vel, ρ, T; quad, @@ -983,7 +955,7 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch. dq_rim_dt += coll.∂ₜL_rim / ρ db_rim_dt += coll.∂ₜB_rim / ρ - # --- Ice self-collection (aggregation) --- + # --- Ice self-collection (aggregation) S_ice_agg = CMP3.ice_self_collection(state, logλ, vel, ρ; quad) dn_ice_dt -= S_ice_agg.dNdt / ρ @@ -1001,87 +973,38 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch. dn_rai_dt += ∂ₜn_ice_melt # Melted ice becomes rain drops dq_ice_dt -= ∂ₜq_ice_melt dn_ice_dt -= ∂ₜn_ice_melt # Ice particles consumed by melting - # Rim mass and rim volume drain proportionally to ice mass during - # melting (uniform-melt assumption — F_rim and ρ_rim are - # preserved through the melt) + # Rim mass and rim volume drain proportionally to ice mass during melting dq_rim_dt -= ∂ₜq_ice_melt * state.F_rim db_rim_dt -= ifelse(state.ρ_rim > 0, ∂ₜq_ice_melt * state.F_rim / state.ρ_rim, zero(FT)) end - # --- ----------------------------- --- - # --- Ice nucleation (F23 + Bigg) --- - # --- ----------------------------- --- - # - # Two parallel heterogeneous-freezing channels, MM15-aligned: - # - # Pathway 1 — F23 deposition / condensation-freezing nucleation - # (analog of Fortran `qinuc`): vapor → pristine ice on - # an INP. F_rim = 0 at genesis; subsequent growth via - # Pathway 5 (MM2015 vapor deposition). - # - # Pathway 2 — F23-bounded Bigg immersion freezing of cloud drops - # (analog of Fortran `qcheti`): drop + INP → fully-rimed - # ice (embryo graupel). F_rim = 1 at genesis; mass and - # number drained from q_lcl, n_lcl. - # + # --- Ice nucleation (F23 + Bigg) τ_act = inp_depletion_model.τ_act - # Vapor deposition nucleation size + # Vapor deposition nucleation size. TODO: put into ClimaParams. D_nuc = FT(10e-6) # 10 μm nascent crystal - small-D tail of the P3 m_nuc = p3.ρ_i * CO.volume_sphere_D(D_nuc) - # F23 INP-activation depletion proxy. With `NIceProxyDepletion` - # (the only model currently provided) this returns `n_ice` — the - # legacy always-on proxy. (A prognostic activation-memory tracer is - # deferred to a follow-up PR.) + # F23 INP-activation depletion proxy. n_active = CM_HetIce.n_active(inp_depletion_model, n_ice) - # ---- Pathway 1: F23 deposition nucleation (vapor → pristine ice) ---- - # - # Target-relaxation form, structurally identical to Cooper-style qinuc - # with F23's log-normal target substituted for Cooper's exponential. - # We route through `CM_HetIce.f23_deposition_rate`, which uses the - # strict-MM15 doc conditions `T < T_freeze − 15 K ∧ S_i > 0.05` by default. - f23_dep = CM_HetIce.f23_deposition_rate( + # --- deposition nucleation (vapor → pristine ice) + dep = CM_HetIce.deposition_rate( ice_nucleation, tps, T, ρ, q_tot, q_lcl + q_rai, q_ice, n_active; m_nuc, τ_act, inpc_log_shift, ) - dn_ice_dt += f23_dep.∂ₜn_frz - dq_ice_dt += f23_dep.∂ₜq_frz - # NO contribution to q_rim, b_rim — pristine deposition crystals have F_rim = 0. - - # ---- Pathway 2: F23-bounded Bigg immersion freezing of cloud drops ---- - # - # Bigg (1953) volume-weighted kinetics over the SB2006 cloud-drop PSD - # (vanilla MM15 / qcheti, both rates ≥ 0 by construction): - # - # ∂ₜn_imm^Bigg = J_bigg(T) · (π/6) · M_D³(N_lcl, λ_c, ν_c, μ_c) - # ∂ₜq_imm^Bigg = J_bigg(T) · ρ_w · (π/6)² · M_D⁶(N_lcl, λ_c, ν_c, μ_c) - # - # F23 INPC imposes an upper bound on the activation rate (the new - # piece for clean-air / OU-SIF regimes; ≥ 0): - # - # ∂ₜn_imm^INPC = INPC / τ_act - # - # Realised rate is the smaller of the two; mass tracks the limiting - # branch via the size-weighted Bigg mass: - # - # ∂ₜn_imm = min(∂ₜn_imm^Bigg, ∂ₜn_imm^INPC) - # ∂ₜq_imm = ∂ₜq_imm^Bigg · ∂ₜn_imm / ∂ₜn_imm^Bigg - # - # Per MM15 (Fortran `qcheti`) frozen cloud drops are fully rimed - # (embryo graupel) → mass adds to BOTH qitot and qirim/birim. + dn_ice_dt += dep.∂ₜn_frz + dq_ice_dt += dep.∂ₜq_frz + # No contribution to q_rim, b_rim — pristine deposition crystals have F_rim = 0. + + # --- F23-bounded Bigg immersion freezing of cloud drops cld_bigg = CM_HetIce.liquid_freezing_rate( mp.ice.rain_freezing, pdf_c, tps, q_lcl, ρ, N_lcl, T, ) - cld_cap = CM_HetIce.f23_immersion_limit_rate( + cld_cap = CM_HetIce.immersion_limit_rate( ice_nucleation, T, ρ; τ = τ_act, inpc_log_shift, n_active, ) ∂ₜn_imm = min(cld_bigg.∂ₜn_frz, cld_cap.∂ₜn_frz) - # Mass = (size-weighted Bigg mass) × (rate ratio). When Bigg is - # silent (off-gate or zero N_lcl/q_lcl), `cld_bigg.∂ₜn_frz` is exactly - # zero — and so is `∂ₜn_imm` — so any safe ratio gives ∂ₜq_imm = 0. - # Use 0 as the safe value to avoid 0/0 NaN. ∂ₜq_imm = ifelse(cld_bigg.∂ₜn_frz > 0, cld_bigg.∂ₜq_frz * ∂ₜn_imm / cld_bigg.∂ₜn_frz, zero(FT)) # Drain liquid: @@ -1093,7 +1016,7 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch. dq_rim_dt += ∂ₜq_imm # F_rim = 1 (frozen drop) db_rim_dt += ∂ₜq_imm / p3.ρ_i # solid-ice rime volume - # --- Ice Sublimation / Deposition --- + # --- Ice Sublimation / Deposition n_per_q_ice = ifelse(q_ice > ϵₘ, n_ice / q_ice, zero(n_ice)) # Deposition/sublimation of cloud ice micro_mock = (; q_tot, q_lcl, q_icl = q_ice, q_rai, q_sno = zero(q_ice)) @@ -1108,27 +1031,13 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch. ∂ₜn_ice_dep = ifelse(∂ₜq_ice_dep < 0, n_per_q_ice * ∂ₜq_ice_dep, zero(∂ₜq_ice_dep)) dq_ice_dt += ∂ₜq_ice_dep dn_ice_dt += ∂ₜn_ice_dep - # Rim mass and rim volume drain proportionally during *sublimation* - # (∂ₜq_ice_dep < 0). Same uniform-melt-style assumption as the melt - # branch: F_rim and ρ_rim are preserved through the phase change, so - # q_rim drains at F_rim · ∂ₜq_ice_sub and b_rim drains at - # (F_rim/ρ_rim) · ∂ₜq_ice_sub. Per MM15 Eqs. for S_qrim and S_Br, - # `−(qrim/qi)·QISUB` and `−(qrim/(ρ_r·qi))·QISUB` (matching the - # `−QIMLT` rim drains we already have). Without this branch, the - # cirrus-level F_rim ≈ 1 band (where ice nucleates, partially - # sublimates, and leaves rim mass behind) would persist - # indefinitely — same root cause as the melt-side bug, just on the - # other side of the freezing line. Deposition (∂ₜq_ice_dep > 0) - # adds pristine ice mass with F_rim = 0, so q_rim/b_rim are - # unaffected; we only drain on the sublimation branch. ∂ₜq_ice_sub = min(∂ₜq_ice_dep, 0) # ≤ 0; zero on the deposition branch dq_rim_dt += ∂ₜq_ice_sub * state.F_rim db_rim_dt += ifelse(state.ρ_rim > 0, ∂ₜq_ice_sub * state.F_rim / state.ρ_rim, zero(FT)) - # --- Ice number adjustment for mass limits --- - # Number adjustment for ice mass limits (Horn 2012, DOI: 10.5194/gmd-5-345-2012). + # --- Ice number adjustment for mass limits # Nudges n_ice toward [q_ice / x_max, q_ice / x_min] over timescale τ. - numadj = (; + numadj = (; # TODO: put into ClimaParams τ = FT(100), x_min = FT(1e-12), # min mean ice particle mass [kg] (~10 μm crystal) x_max = FT(1e-5), # max mean ice particle mass [kg] (~5 mm aggregate) @@ -1136,8 +1045,7 @@ to be non-Nothing, eliminating runtime type checks and dynamic dispatch. ∂ₜn_ice_numadj = CM2.number_tendency_from_mass_limits(numadj, q_ice, n_ice) dn_ice_dt += ∂ₜn_ice_numadj - # --- Rain Heterogeneous Freezing (Bigg 1953) --- - # Immersion freezing of rain drops, as in Morrison & Milbrandt (2015). + # --- Rain Heterogeneous Freezing (Bigg 1953) rain_frz = CM_HetIce.liquid_freezing_rate(mp.ice.rain_freezing, pdf_r, tps, q_rai, ρ, N_rai, T) # Rain → ice (frozen rain is fully rimed, per MM15) diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl index cad2a20ec..bda69cda8 100644 --- a/src/IceNucleation.jl +++ b/src/IceNucleation.jl @@ -20,8 +20,8 @@ export P3_het_N_i export INP_concentration_frequency export INP_concentration_mean export liquid_freezing_rate -export f23_immersion_limit_rate -export f23_deposition_rate +export immersion_limit_rate +export deposition_rate export n_active """ @@ -248,12 +248,12 @@ function INP_concentration_mean((; T_freeze)::CMP.Frostenberg2023, T) end """ - liquid_freezing_rate(parameterization, pdf, tps, q_rai, ρ, N_rai, T) + liquid_freezing_rate(opt, pdf, tps, q_rai, ρ, N_rai, T) Compute the rate of liquid water freezing into ice. # Arguments - - `parameterization`: The [`CMP.RainFreezing`](@ref) parameterization. + - `opt`: The [`CMP.RainFreezing`](@ref) parameterization. - `pdf`: The liquid water particle size distribution (PSD) PDF. - `tps`: Thermodynamics parameters. - `q`: Liquid water specific content [kg(water) kg⁻¹(air)]. @@ -266,7 +266,7 @@ Compute the rate of liquid water freezing into ice. + `∂ₜn_frz`: Specific number freezing rate [kg⁻¹(air) s⁻¹]. + `∂ₜq_frz`: Specific mass freezing rate [kg(water) kg⁻¹(air) s⁻¹]. """ -function liquid_freezing_rate(parameterization::CMP.RainFreezing, pdf, tps, q, ρ, N, T) +function liquid_freezing_rate(opt::CMP.RainFreezing, pdf, tps, q, ρ, N, T) FT = eltype(q) T_freeze = TDI.TD.Parameters.T_freeze(tps) (; ρw) = pdf # [kg(water) m⁻³(water)] @@ -276,7 +276,7 @@ function liquid_freezing_rate(parameterization::CMP.RainFreezing, pdf, tps, q, (; Dr_mean) = CM2.pdf_rain_parameters(pdf, q, ρ, N) # Bigg (1953) volumetric freezing rate: - J_bigg = parameterization(T, T_freeze) # [m⁻³(water) s⁻¹] + J_bigg = opt(T, T_freeze) # [m⁻³(water) s⁻¹] # Diameter PSD moments via exponential_Mⁿ: # M_D^k = ∫ D^k n(D) dD @@ -306,7 +306,10 @@ function liquid_freezing_rate(parameterization::CMP.RainFreezing, pdf, tps, q, end """ - liquid_freezing_rate(parameterization, pdf::CMP.CloudParticlePDF_SB2006, tps, q, ρ, N, T) + liquid_freezing_rate( + opt::CMP.RainFreezing, pdf::CMP.CloudParticlePDF_SB2006, + tps, q, ρ, N, T, + ) Compute the rate of cloud-droplet immersion freezing into ice using the same Bigg (1953) kinetics as the rain version, but integrated over the @@ -322,16 +325,15 @@ Bigg's per-drop freezing probability is `J_bigg(T) · (π/6) · D³`. Integratin against the PSD gives closed-form number- and mass-freezing rates: ``` -∂ₜn_frz = J_bigg · (π/6) · M_D³(N₀c, λc, νcD, μcD) +∂ₜn_frz = J_bigg · (π/6) · M_D³(N₀c, λc, νcD, μcD) ∂ₜq_frz = J_bigg · ρw · (π/6)² · M_D⁶(N₀c, λc, νcD, μcD) ``` -with `M_Dᵏ` the kth diameter moment computed by -`generalized_gamma_Mⁿ`. Volume-weighting → bigger drops freeze -first, captured analytically; no numerical PSD integration required. +with `M_Dᵏ` the kth diameter moment computed by `generalized_gamma_Mⁿ`. +Volume-weighting → bigger drops freeze first # Arguments - - `parameterization`: The [`CMP.RainFreezing`](@ref) parameterization + - `opt`: The [`CMP.RainFreezing`](@ref) parameterization (Bigg / Barklie-Gokhale parameters; the `Rain` in the name is historical — the kinetics apply to any liquid-drop PSD). - `pdf`: The [`CMP.CloudParticlePDF_SB2006`](@ref) cloud-droplet PSD. @@ -347,7 +349,7 @@ first, captured analytically; no numerical PSD integration required. + `∂ₜq_frz`: Specific mass freezing rate [kg(water) kg⁻¹(air) s⁻¹]. """ function liquid_freezing_rate( - parameterization::CMP.RainFreezing, pdf::CMP.CloudParticlePDF_SB2006, + opt::CMP.RainFreezing, pdf::CMP.CloudParticlePDF_SB2006, tps, q, ρ, N, T, ) FT = eltype(q) @@ -359,12 +361,11 @@ function liquid_freezing_rate( (; λc, νcD, μcD) = CM2.pdf_cloud_parameters(pdf, q, ρ, N) # Bigg (1953) volumetric freezing rate per unit drop water volume. - J_bigg = parameterization(T, T_freeze) # [m⁻³(water) s⁻¹] + J_bigg = opt(T, T_freeze) # [m⁻³(water) s⁻¹] # Diameter moments via the closed-form generalized-gamma formula. # When N or L is essentially zero, `pdf_cloud_parameters` returns - # `λc = Inf`, which makes `B^(-k/μ) → 0`, so the moments vanish — no - # extra guard needed. + # `λc = Inf`, which makes `B^(-k/μ) → 0`, so the moments vanish M_D³ = DT.generalized_gamma_Mⁿ(νcD, μcD, λc, n, 3) # [m³ · kg⁻¹(air)] M_D⁶ = DT.generalized_gamma_Mⁿ(νcD, μcD, λc, n, 6) # [m⁶ · kg⁻¹(air)] @@ -372,9 +373,9 @@ function liquid_freezing_rate( ∂ₜn_frz = J_bigg * V_1 * M_D³ # [kg⁻¹(air) s⁻¹] ∂ₜq_frz = J_bigg * ρw * V_1^2 * M_D⁶ # [kg(water) kg⁻¹(air) s⁻¹] - # Same gates as the rain branch: non-trivial number/mass and T < -4 °C. + # Check non-trivial number/mass and that T < -4 °C. ϵₘ, ϵₙ = UT.ϵ_numerics_2M_M(FT), UT.ϵ_numerics_2M_N(FT) - cond = (n > ϵₙ) & (q > ϵₘ) & (T < T_freeze - 4) + cond = (n > ϵₙ) & (q > ϵₘ) & (T < T_freeze - 4) # TODO: Make a parameter. ∂ₜn_frz = ifelse(cond, ∂ₜn_frz, zero(n)) ∂ₜq_frz = ifelse(cond, ∂ₜq_frz, zero(q)) @@ -382,10 +383,9 @@ function liquid_freezing_rate( end """ - f23_immersion_limit_rate(f23_params, T, ρ; τ, inpc_log_shift, n_active) + immersion_limit_rate(opt::CMP.Frostenberg2023, T, ρ; τ, inpc_log_shift, n_active) -Compute the **F23-INPC-imposed upper limit** on the cloud-droplet immersion -freezing number rate. +Compute the F23-INPC-imposed upper limit on the cloud-droplet immersion freezing number rate. The Frostenberg 2023 climatology specifies a target INP concentration `INPC(T)` in air [m⁻³(air)]. Treating that target as a budget that should be activated @@ -396,49 +396,41 @@ kg of air per second is ∂ₜn_lim = max(0, INPC(T)/ρ - n_active) / τ. ``` -`n_active` is the depletion proxy supplied by the host (see -`AbstractINPDepletion` and [`n_active`](@ref)): with the -`NIceProxyDepletion` model the host passes `n_ice` (or pass zero to -recover the no-depletion form). - -This is the "INPC-only" leg of Pathway 2 — the cap that forces the realized -immersion-freezing rate to fall below pure Bigg kinetics in clean-air or -low-q_lcl regimes. Combined with [`liquid_freezing_rate`](@ref) on the cloud -PSD, the realized Pathway-2 rate is `min(bigg, ∂ₜn_lim)`. +`n_active` is the depletion proxy; see [`n_active`](@ref). # Arguments - - `f23_params`: The [`CMP.Frostenberg2023`](@ref) parameters. + - `opt`: The [`CMP.Frostenberg2023`](@ref) parameters. - `T`: Air temperature [K]. - `ρ`: Air density [kg(air) m⁻³(air)]. # Keyword arguments - `τ`: Relaxation timescale [s] (default `300`). - - `inpc_log_shift`: Additive shift to `log(INPC)` (e.g. an OU-SIF stochastic - excursion). Default `0`. - - `n_active`: Depletion proxy [kg⁻¹(air)] (default `0`, i.e. no depletion; - recovers the no-memory cap). + - `inpc_log_shift`: Additive shift to `log(INPC)` Default `0`. + - `n_active`: Depletion proxy [kg⁻¹(air)] # Returns - A `NamedTuple` `(; ∂ₜn_frz)` — the specific number freezing-rate cap [kg⁻¹(air) s⁻¹]. Zero when `T ≥ T_freeze`. """ -function f23_immersion_limit_rate( - f23_params::CMP.Frostenberg2023, T, ρ; +function immersion_limit_rate( + opt::CMP.Frostenberg2023, T, ρ; τ = oftype(T, 300), inpc_log_shift = zero(T), n_active = zero(T), ) - T ≥ f23_params.T_freeze && return (; ∂ₜn_frz = zero(T)) - log_inpc = INP_concentration_mean(f23_params, T) + inpc_log_shift + T ≥ opt.T_freeze && return (; ∂ₜn_frz = zero(T)) + log_inpc = INP_concentration_mean(opt, T) + inpc_log_shift INPC_per_kg = exp(log_inpc) / ρ # [kg⁻¹(air)] ∂ₜn_frz = max(zero(T), INPC_per_kg - n_active) / τ # [kg⁻¹(air) s⁻¹] return (; ∂ₜn_frz) end """ - f23_deposition_rate(f23_params, tps, T, ρ, q_vap, n_ice; - τ, inpc_log_shift, T_thresh, S_i_thresh, ρ_i, D_nuc) + deposition_rate( + opt::CMP.Frostenberg2023, tps, T, ρ, q_vap, n_ice; + τ, inpc_log_shift, T_thresh, S_i_thresh, ρ_i, D_nuc, + ) -Compute the **F23 deposition nucleation rate** (Pathway 1). +Compute the Frostenberg 2023 deposition nucleation rate The rate is the Frostenberg 2023 INP concentration (treated as a budget) relaxed toward depletion at `n_ice` over timescale `τ`: @@ -455,31 +447,22 @@ m_starter = (π/6) · ρ_i · D_nuc³ (default ≈ 4.8e-13 kg for a 10 μm solid-ice crystal). The mass tendency is the implied mass injection, capped by half the local vapor excess over -ice saturation per relaxation window — this prevents the channel from -sourcing more vapor than physically available when OU-SIF excursions push -INPC above the local supply: +ice saturation per relaxation window: ``` q_excess = max(0, q_vap - q_sat_ice) ∂ₜq_frz = min(m_starter · ∂ₜn_frz, ½ q_excess / τ) ``` -`q_sat_ice` is computed internally from `tps` via -`TDI.saturation_vapor_specific_content_over_ice(tps, T, ρ)`. - -Two physical gates close both rates outside the activation window: +Conditions for activation: - `T < T_thresh` (default `T_freeze - 15 K`, i.e. -15 °C): below this temperature deposition nucleation is active. - `S_i ≡ q_vap/q_sat_ice - 1 > S_i_thresh` (default 5 %): ice supersaturation is required for vapor to nucleate onto INPs. -The `n_ice` argument is a proxy for already-activated INPs; in a Phillips- -style framework it would be replaced by an explicit prognostic -`N_INP_active`. - # Arguments - - `f23_params`: The [`CMP.Frostenberg2023`](@ref) parameters. + - `opt`: The [`CMP.Frostenberg2023`](@ref) parameters. - `tps`: Thermodynamics parameters (used for the ice-saturation curve). - `T`: Air temperature [K]. - `ρ`: Air density [kg(air) m⁻³(air)]. @@ -490,34 +473,31 @@ style framework it would be replaced by an explicit prognostic # Keyword arguments - `τ`: Relaxation timescale [s] (default `300`). - `inpc_log_shift`: Additive shift to `log(INPC)` (default `0`). - - `T_thresh`: Activation temperature threshold [K] (default - `f23_params.T_freeze - 15`). + - `T_thresh`: Activation temperature threshold [K] (default `opt.T_freeze - 15`). - `S_i_thresh`: Activation ice-supersaturation threshold (default `0.05`). - - `ρ_i`: Solid-ice density used for the starter mass [kg/m³] (default - `916.7`). - - `D_nuc`: Nascent crystal diameter [m] (default `10e-6`, the small-D - tail of the P3 distribution). + - `ρ_i`: Solid-ice density used for the starter mass [kg/m³] (default `916.7`). + - `D_nuc`: Nascent crystal diameter [m] (default `10e-6`, the small-D tail of the P3 distribution). # Returns - A `NamedTuple` `(; ∂ₜn_frz, ∂ₜq_frz)` with the specific number rate [kg⁻¹(air) s⁻¹] and specific mass rate [kg(ice) kg⁻¹(air) s⁻¹]. Zero outside the activation window. """ -function f23_deposition_rate( - f23_params::CMP.Frostenberg2023, tps, T, ρ, q_tot, q_liq, q_ice, n_ice; m_nuc, - T_thresh = f23_params.T_freeze - 15, S_i_thresh = oftype(f23_params.T_freeze, 0.05), - τ_act = 300, inpc_log_shift = 0, +function deposition_rate( + opt::CMP.Frostenberg2023, tps, T, ρ, q_tot, q_liq, q_ice, n_ice; m_nuc, + T_thresh = opt.T_freeze - 15, S_i_thresh = oftype(opt.T_freeze, 0.05), + τ_act = 300, inpc_log_shift = 0, # TODO: Put these in ClimaParams ) q_sat_ice = TDI.saturation_vapor_specific_content_over_ice(tps, T, ρ) q_vap = TDI.q_vap(q_tot, q_liq, q_ice) S_i = q_vap / q_sat_ice - 1 - # Activation gates. - gated = (T < T_thresh) & (S_i > S_i_thresh) + # Activation conditions + cond = (T < T_thresh) & (S_i > S_i_thresh) # Nucleation rate, limited by ambient INP availability and relaxation time τ_act. - log_inpc = INP_concentration_mean(f23_params, T) + inpc_log_shift + log_inpc = INP_concentration_mean(opt, T) + inpc_log_shift INPC_per_kg = exp(log_inpc) / ρ ∂ₜn_frz = max(0, INPC_per_kg - n_ice) / τ_act - ∂ₜn_frz = ifelse(gated, ∂ₜn_frz, zero(∂ₜn_frz)) + ∂ₜn_frz = ifelse(cond, ∂ₜn_frz, zero(∂ₜn_frz)) # Vapor-excess cap on the implied mass injection. q_excess = max(0, q_vap - q_sat_ice) # Implied mass injection, capped by half the local vapor excess per relaxation window. @@ -530,10 +510,10 @@ end # --------------------------------------------------------------------------- """ - n_active(model::CMP.AbstractINPDepletion, n_ice) + n_active(model::CMP.NIceProxyDepletion, n_ice) Return the depletion proxy `n_active` to subtract from the F23 INPC -target in [`f23_deposition_rate`](@ref) and any analogous INPC-budgeted +target in [`deposition_rate`](@ref) and any analogous INPC-budgeted rate. For `NIceProxyDepletion` (the only model currently provided) this is the in-cell ice number `n_ice`. (A prognostic activation-memory model returning a host-supplied tracer is deferred to a follow-up PR.) diff --git a/src/Microphysics2M.jl b/src/Microphysics2M.jl index 4c83c6937..986326ab0 100644 --- a/src/Microphysics2M.jl +++ b/src/Microphysics2M.jl @@ -67,18 +67,6 @@ Return the parameters of the rain drop diameter distribution """ function pdf_rain_parameters(pdf_r::CMP.RainParticlePDF_SB2006_notlimited, qᵣ, ρₐ, Nᵣ) FT = eltype(qᵣ) - # Degenerate-input guard. The unlimited PDF inverts `xr_mean = Lᵣ/Nᵣ` - # and `Dr_mean = 1/cbrt(π·ρw/xr_mean)`; if EITHER `Nᵣ` or `qᵣ` is - # non-positive (which explicit RK sub-stages transiently produce when - # number is depleted faster than mass, or vice-versa) the result is a - # *negative* `Dr_mean` — finite and non-zero, so it slips past the - # downstream `iszero(Dr_mean)` checks and makes `exponential_quantile` - # throw `DomainError("Mean parameter must be positive")`. Return zero - # params instead (consumers already special-case `iszero(N₀r)`). The old - # guard only fired when BOTH were zero (`&&`); broaden to `||` with the - # same `ϵ` thresholds as the cloud PDF's `log_pdf_cloud_parameters_mass` - # guard (these `ϵ = eps(FT) > 0`, so all non-positive inputs are caught). - # Strict no-op for physically valid (above-ϵ) inputs. (Nᵣ < UT.ϵ_numerics_2M_N(FT) || qᵣ < UT.ϵ_numerics_2M_M(FT)) && return (; N₀r = zero(Nᵣ), Dr_mean = zero(qᵣ), xr_mean = zero(qᵣ)) (; ρw) = pdf_r diff --git a/src/P3_particle_properties.jl b/src/P3_particle_properties.jl index 6bb271938..3ae12f8c8 100644 --- a/src/P3_particle_properties.jl +++ b/src/P3_particle_properties.jl @@ -10,20 +10,9 @@ This struct bundles the P3 parameterizations `params`, the provided rime state # Construction -Two positional, broadcast-friendly, GPU-clean entry points; both wrap the -6-field inner constructor and auto-compute `thresholds` so the cached -derivatives stay consistent with the rime state. - - - `P3State(params, L_ice, N_ice, F_rim, ρ_rim)` — direct construction - from explicit `(F_rim, ρ_rim)`. Used by tests, audit prototypes, and - documentation plots that want exact rime-state inputs. **Performs no - validation.** Out-of-range inputs (`F_rim ∉ [0, 1)`, `ρ_rim ∉ [0, - ρ_l]`) produce a state whose downstream evaluation typically yields - NaN or non-physical values. - - [`state_from_prognostic`](@ref) — production / upstream-model entry - point. Accepts the volumetric prognostic variables `(ρq_ice, - ρn_ice, ρq_rim, ρb_rim)`, regularises them into `(F_rim, ρ_rim)`, - and returns the constructed state. + - [`state_from_prognostic`](@ref): Main entry point. + Accepts the volumetric prognostic variables `(ρq_ice, ρn_ice, ρq_rim, ρb_rim)`, + regularises them into `(F_rim, ρ_rim)`, and returns the constructed state. # Fields $(FIELDS) @@ -51,14 +40,6 @@ struct P3State{FT, PARAMS <: CMP.ParametersP3} D_cr::FT end -# Positional 5-arg constructor — derives the cached thresholds and -# graupel density from `(F_rim, ρ_rim)`. The production hot path -# (`state_from_prognostic`) wraps this with input regularisation; -# tests / docs call this directly with explicit `(F_rim, ρ_rim)`. -# -# When `F_rim = 0`, `ρ_g` collapses to `NaN` and `D_gr = D_cr = Inf` -# (the "no graupel regime" sentinel). Downstream code is expected to -# gate on the unrimed branch before reading those fields. function P3State(params::CMP.ParametersP3, ρq_ice, ρn_ice, F_rim, ρ_rim) FT = eltype(ρq_ice) (; mass, ρ_i) = params @@ -80,33 +61,31 @@ ShowMethods.field_units(::P3State) = (; """ state_from_prognostic(params, ρq_ice, ρn_ice, ρq_rim, ρb_rim) -Construct a [`P3State`](@ref) from the volumetric prognostic ice -variables directly, computing the (clamped, regularised) rime mass -fraction and rime density. **The production / upstream-model entry -point.** Positional and broadcast-friendly. +Construct a [`P3State`](@ref) from the volumetric prognostic ice variables directly, +computing the (clamped, regularised) rime mass fraction and rime density. -The regularised ratios come from `Utilities.rime_mass_fraction` and -`Utilities.rime_density`, which smoothly go to zero when their -denominators are near machine precision — avoiding the hard discontinuity -at `q_ice = ϵ` / `b_rim = ϵ`. The upper clamps `F_rim < 1 − ε` and +The regularised ratios come from [`UT.rime_mass_fraction`](@ref) and +[`UT.rime_density`](@ref), which smoothly go to zero when their +denominators are near machine precision, avoiding the discontinuity +at `q_ice = ϵ` / `b_rim = ϵ`. The upper clamps `F_rim < 1 - ε` and `ρ_rim ≤ 0.8·ρ_l ≈ 730 kg/m³` keep the result inside the domain of validity of the threshold formulas evaluated by the [`P3State`](@ref) constructor. !!! note "TODO — revisit the `ρ_rim ≤ 0.8·ρ_l` cap" - The closed-form graupel density `ρ_g = F_rim·ρ_rim + (1−F_rim)·ρ_d` + The closed-form graupel density `ρ_g = F_rim·ρ_rim + (1-F_rim)·ρ_d` can mathematically exceed `ρ_l` — `ρ_d` (the unrimed portion's density, [`get_ρ_d`](@ref)) is linear in `ρ_rim` with no built-in upper clamp, so feeding `ρ_rim` near `ρ_l` can produce `ρ_g > ρ_l`. That breaks the threshold ordering `D_th < D_gr < D_cr` that the P3 - partitioning assumes (`D_gr ∝ (6α_va/(π·ρ_g))^{1/(3−β_va)}` shrinks - as `ρ_g` grows; eventually `D_gr < D_th`). The 0.8 cushion keeps - `ρ_g` comfortably below `ρ_l` for the realistic `(F_rim, ρ_rim)` - regime — Macklin-type rime-density formulations rarely give - `ρ_rim > 700 kg/m³` anyway, so the cushion costs us nothing - physically. To lift the cap to `ρ_l` we'd need to (i) explicitly - bound `ρ_g` (e.g. `min(ρ_g, ρ_l)`) or rederive `ρ_d` so it's - monotone-bounded by `ρ_l`, and (ii) accept that bulk rime densities - 800–917 kg/m³ are off the calibration domain of the original P3 fit. + partitioning assumes (`D_gr ∝ (6α_va/(π·ρ_g))^{1/(3-β_va)}` shrinks + as `ρ_g` grows; eventually `D_gr < D_th`). The 0.8-factor keeps + `ρ_g` comfortably below `ρ_l` for the realistic `(F_rim, ρ_rim)` regime. + Rime Density formulations structured like Macklin (1962) rarely give + `ρ_rim > 700 kg/m³` anyway, so the upper bound is usually inert. + To lift the cap to `ρ_l` we'd need to (i) explicitly bound `ρ_g` + (e.g. `min(ρ_g, ρ_l)`) or rederive `ρ_d` so it's monotone-bounded by + `ρ_l`, and (ii) accept that bulk rime densities 800-917 kg/m³ are off + the calibration domain of the original P3 fit. # Arguments - `params`: [`CMP.ParametersP3`](@ref) @@ -118,7 +97,7 @@ of the threshold formulas evaluated by the [`P3State`](@ref) constructor. function state_from_prognostic(params::CMP.ParametersP3, ρq_ice, ρn_ice, ρq_rim, ρb_rim) FT = eltype(ρq_ice) F_rim = min(UT.rime_mass_fraction(ρq_rim, ρq_ice), one(FT) - eps(FT)) - ρ_rim = min(UT.rime_density(ρq_rim, ρb_rim), FT(0.8) * params.ρ_l) + ρ_rim = min(UT.rime_density(ρq_rim, ρb_rim), FT(0.8) * params.ρ_l) # TODO: Make this limit configurable return P3State(params, ρq_ice, ρn_ice, F_rim, ρ_rim) end @@ -147,23 +126,6 @@ Exact solution for the density of the unrimed portion of the particle as # Returns - `ρ_d`: density of the unrimed portion of the particle [kg/m³] - -# Examples - -```jldoctest -julia> import CloudMicrophysics.Parameters as CMP, - ClimaParams as CP, - CloudMicrophysics.P3Scheme as P3 - -julia> FT = Float64; - -julia> mass = CMP.MassPowerLaw(CP.create_toml_dict(FT)); - -julia> F_rim, ρ_rim = FT(0.5), FT(916.7); - -julia> ρ_d = P3.get_ρ_d(mass, F_rim, ρ_rim) -488.9120789986412 -``` """ function get_ρ_d((; β_va)::CMP.MassPowerLaw, F_rim, ρ_rim) k = (1 - F_rim)^(-1 / (3 - β_va)) @@ -248,8 +210,8 @@ get_D_cr(mass::CMP.MassPowerLaw, F_rim, ρ_g) = _get_threshold(mass, ρ_g * (1 - """ segment_boundaries(state::P3State, D_min = 0, D_max = Inf) -Return the flat 5-tuple `(D_min, D_th, D_gr, D_cr, D_max)` of P3 mass- -regime boundaries clamped into the requested integration window +Return the 5-tuple `(D_min, D_th, D_gr, D_cr, D_max)` of P3 mass-regime +boundaries clamped into the requested integration window `[D_min, D_max]`. Suitable as the `bnds` argument to [`integrate`](@ref) / [`subintervals`](@ref). diff --git a/src/P3_processes.jl b/src/P3_processes.jl index 6b5d87aa6..7f943522a 100644 --- a/src/P3_processes.jl +++ b/src/P3_processes.jl @@ -93,10 +93,35 @@ function ice_melt( return (; dNdt, dLdt) end -function collision_cross_section_ice_liquid(state, Dᵢ, Dₗ) - rᵢ_eff(Dᵢ) = √(ice_area(state, Dᵢ) / π) - return π * (rᵢ_eff(Dᵢ) + Dₗ / 2)^2 # collision cross section -end +""" + collision_cross_section_ice_liquid_coeffs(rᵢ) + collision_cross_section_ice_liquid_coeffs(state, Dᵢ) + +Monomial coefficients `(k₀, k₁, k₂)` of the ice-liquid collision cross-section as +a polynomial in the liquid diameter `Dₗ`, + +```math +σ(Dᵢ, Dₗ) = π (rᵢ + Dₗ/2)² = k₀ + k₁ Dₗ + k₂ Dₗ², +``` + +with `k₀ = π rᵢ²`, `k₁ = π rᵢ`, `k₂ = π/4`, where the ice effective radius +is `rᵢ = √(ice_area(state, Dᵢ)/π)`; see [`ice_area`](@ref). + +Used in [`collision_cross_section_ice_liquid`](@ref) +""" +@inline collision_cross_section_ice_liquid_coeffs(rᵢ::FT) where {FT} = + (π * rᵢ^2, π * rᵢ, FT(π / 4)) +@inline collision_cross_section_ice_liquid_coeffs(state, Dᵢ) = + collision_cross_section_ice_liquid_coeffs(√(ice_area(state, Dᵢ) / π)) + +""" + collision_cross_section_ice_liquid(state, Dᵢ, Dₗ) + +Ice-liquid collision cross-section [m²], `π (rᵢ(Dᵢ) + Dₗ/2)²`, evaluated by +Horner from the shared [`collision_cross_section_ice_liquid_coeffs`](@ref). +""" +collision_cross_section_ice_liquid(state, Dᵢ, Dₗ) = + evalpoly(Dₗ, collision_cross_section_ice_liquid_coeffs(state, Dᵢ)) """ volumetric_collision_rate_integrand(state, velocity_params, ρₐ) @@ -290,213 +315,69 @@ function get_liquid_integrals(n, ∂ₜV, m_liq, ρ′_rim, liq_bounds; quad = C return liquid_integrals end -# --------------------------------------------------------------------------- -# Closed-form rain inner integral (A3 / issue 003 / A4 Phase 3) -# -# The rain SB2006 PSD is pure-exponential `n_r(D)=N₀r e^{-D/Dr_mean}` -# and Chen-2022 rain `v_l(D)=Σⱼ Aⱼ D^{Bⱼ} e^{-Cⱼ D}`; the cross-section -# is a degree-2 poly in `Dₗ`. The only obstruction `|v_i(Dᵢ)−v_l(Dₗ)|` -# is removed EXACTLY by splitting the inner integral at the velocity -# crossover `D*(Dᵢ)` where `v_l(D*)=v_i(Dᵢ)`. Each piece is a finite sum -# of regularized-incomplete-gamma moments. This is exact (no physical -# approximation); the only numeric is the scalar `D*` root, found to -# ~1e-12 — far tighter than the ~1e-3 bulk-tendency budget. Validated -# vs ChebyshevGauss(1024): max relerr 5.4e-6 (N) / 2.0e-6 (M) over a -# 1944-cell grid; AD-clean under the frozen-logλ/Float64-param regime -# (`gamma_inc` 2nd-arg only). B-rim stays numerical (piecewise-rational -# `1/ρ′_rim`). See `audits/p3_audit/A3_closed_form_reductions.md` §4. - -# The shared incomplete-gamma partial-moment kernel `gamma_inc_moment` -# (the linear-space twin of `loggamma_inc_moment`, used by the -# `closed_rain_inner_NM` sign-alternating split) is defined next to -# `loggamma_inc_moment` in `P3_size_distribution.jl` so the log/linear -# pair stays co-located and in sync (P1-2). See its docstring there. - """ crossover_diameter(v_target, v_l, D_min, D_max) -`D* ∈ [D_min,D_max]` with `v_l(D*) = v_target`, by bisection. -`v_l` is the rain terminal-velocity callable, passed in from the -**canonical** CM API `CO.particle_terminal_velocity(vel.rain, ρₐ)` — -the *exact same* Chen-2022 rain `v(D)` the numerical fallback uses in -`volumetric_collision_rate_integrand` / `compute_local_rime_density` -(`v_liq = CO.particle_terminal_velocity(velocity_params.rain, ρₐ)`). -Reusing that callable here (rather than re-rolling the -`CO.Chen2022_monodisperse_pdf` summation) guarantees the closed-form -root solve cannot silently diverge from CM's Chen rain law if that -table/form ever changes. -Returns `D_min` if `v_target ≤ v_l(D_min)` / `D_max` if -`v_target ≥ v_l(D_max)` (the `|·|` never flips in the bracket). -Chen-2022 rain `v_l` is monotone except a tiny ρₐ-dependent dimple at -`D≈5–9 mm` with rel amplitude ≤1e-4 at the extreme PSD tail; measured -single-`D*` worst-case error ≤1.4e-6 (below the bulk budget), so a -single root is used. See A3 §3.4 for the general (deep-dimple / -non-Chen) multi-root fallback (defensive — not reached for Chen rain). -""" -function crossover_diameter( - v_target, v_l::F, D_min, D_max; - tol = nothing, max_iters::Int = 60, -) where {F} - # AD-generic: `v_target` (∝ v̄ᵢ, Dual via state) and the bracket - # `D_min,D_max` (Dual via the rain PSD) need not share a type — work - # in their promotion (the over-tie was the latent AD blocker, P1-3). - T = float(promote_type(typeof(v_target), typeof(D_min), typeof(D_max))) - tol_ = tol === nothing ? T(1e-12) : T(tol) - D_min, D_max = T(D_min), T(D_max) +Find the diameter `D` in `[D_min, D_max]` where `v_l(D) = v_target` +""" +function crossover_diameter(v_target, v_l::F, D_min, D_max) where {F} + FT = float(promote_type(typeof(v_target), typeof(D_min), typeof(D_max))) f(D) = v_l(D) - v_target - flo, fhi = f(D_min), f(D_max) - flo > 0 && return D_min # v_target below the bracket — no flip - fhi < 0 && return D_max # v_target above the bracket — no flip - a, b = D_min, D_max - for _ in 1:max_iters - m = (a + b) / 2 - fm = f(m) - b - a < tol_ * max(b, one(T)) && return m - if (flo < 0) == (fm < 0) - a, flo = m, fm - else - b, fhi = m, fm - end - end - return (a + b) / 2 + sol = RS.find_zero(f, + RS.BrentsMethod(FT(D_min), FT(D_max)), RS.CompactSolution(), + RS.SolutionTolerance(FT(1e-12)), 60, + ) + return sol.root end """ - closed_rain_inner_NM(Dᵢ, v_i_at_Dᵢ, v_l, rᵢ, ρ_w, ai, bi, ci, - D_min, D_max, N₀r, D̄r) - -Closed-form `(∂ₜN_col, ∂ₜM_col)` for the rain inner integral at one -outer `Dᵢ`. `rᵢ = √(ice_area(state,Dᵢ)/π)`; `v_l` is the canonical -rain terminal-velocity callable `CO.particle_terminal_velocity(vel.rain, -ρₐ)` (reused for the `D*` crossover root — see [`crossover_diameter`]); -`(ai,bi,ci) = CO.Chen2022_vel_coeffs(vel.rain, ρₐ)` are the *same* Chen -coefficients that callable is built from, kept here because the closed -form's `αⱼ=α0+cⱼ` exponent merge and the `Aⱼ`/`Bⱼ+m` partial-moment -weights need the per-term `(aⱼ,bⱼ,cⱼ)` split, for which CM exposes no -existing API (`CO.Chen2022_exponential_pdf` is the *full-domain* `[0,∞)` -1M analogue — not a drop-in here, the `D*` split needs *partial* moments). -Rain PSD `n_r(D)=N₀r e^{-D/D̄r}`. -""" -function closed_rain_inner_NM( - Dᵢ, v_i_at_Dᵢ, v_l::F, rᵢ, ρ_w, ai, bi, ci, D_min, D_max, N₀r, D̄r, -) where {F} - # AD-generic working type. In the implicit-2M+P3 Jacobian (A2) the - # `Dᵢ` quadrature node is Float64 while the state/PSD-derived inputs - # (`v̄ᵢ, rᵢ, N₀r, D̄r`, the split bounds) are `Dual`; the inputs are a - # *mix* of Float64 and Dual, so the working type is their promotion - # (no `where {T}` over-tie — that was the latent AD blocker the - # harness's central-FD evidence masked; P1-3 / P1-test gap 2). - T = float( - promote_type( - typeof(v_i_at_Dᵢ), typeof(rᵢ), typeof(ρ_w), - typeof(D_min), typeof(D_max), typeof(N₀r), typeof(D̄r), - eltype(ai), - ), + closed_rain_inner_NM( + Dᵢ, v_i_at_Dᵢ, v_l, rᵢ, ρw, ai, bi, ci, D_min, D_max, N₀r, Dr_mean, ) - # `Tp` is the *non-Dual* exponent type (the Chen-coeff float type): - # the `gamma_inc_moment` moment order `p` (hence `z = p+1`) MUST stay - # a plain float, never a `Dual` — that is exactly the A2-safe - # `SF.gamma_inc(z, α·D)` "2nd-arg-Dual only" invariant (a Dual - # *first* arg `z` hits the missing `_gamma_inc` rule). The - # cross-section / Chen / m_liq exponents are physical constants and - # are never differentiated, so this is correct, not a workaround. - Tp = float(eltype(ai)) - K0 = T(π) * rᵢ^2 - K1 = T(π) * rᵢ - K2 = T(π) / 4 - α0 = inv(D̄r) + +Closed-form `(∂ₜN_col, ∂ₜM_col)` for the rain inner integral at one outer `Dᵢ`. +""" +function closed_rain_inner_NM(Dᵢ, v_i_at_Dᵢ, v_l::F, rᵢ, ρw, ai, bi, ci, D_min, D_max, N₀r, Dr_mean) where {F} + FT = float(eltype(ai)) + λ = inv(Dr_mean) # rain PSD slope: n_r(D) ∝ e^{-λ D} Dstar = crossover_diameter(v_i_at_Dᵢ, v_l, D_min, D_max) - mfac = ρ_w * T(π) / 6 - # `Tt`/`Tpt` are passed as `::Type` ARGUMENTS, not captured from the - # enclosing scope. A computed *type value* (`T`, `Tp` above) captured by an - # inner closure is boxed by Julia as `::DataType` (it loses its `Type{FT}` - # precision), which makes `zero(Tt)`, `Tpt(0)`, … infer `::Any` and widens - # this closure's return to `Tuple{Any,Any}` — the root of the 1.10 JET - # failures and the GPU dynamic-dispatch / mpfr InvalidIRError on the - # closed-form rain path. Taking them as type arguments keeps them concrete. - function piece_NM(::Type{Tt}, ::Type{Tpt}, a, b, sgn) where {Tt, Tpt} - b > a || return (zero(Tt), zero(Tt)) - vi_part_N = - sgn * v_i_at_Dᵢ * - ( - K0 * gamma_inc_moment(a, b, Tpt(0), α0) + - K1 * gamma_inc_moment(a, b, Tpt(1), α0) + - K2 * gamma_inc_moment(a, b, Tpt(2), α0) - ) - vi_part_M = - sgn * v_i_at_Dᵢ * mfac * - ( - K0 * gamma_inc_moment(a, b, Tpt(3), α0) + - K1 * gamma_inc_moment(a, b, Tpt(4), α0) + - K2 * gamma_inc_moment(a, b, Tpt(5), α0) - ) - v_part_N = zero(Tt) - v_part_M = zero(Tt) - @inbounds for j in eachindex(ai) - αj = α0 + ci[j] - Aj = ai[j] - Bj = bi[j] - v_part_N += - -sgn * Aj * - ( - K0 * gamma_inc_moment(a, b, Bj, αj) + - K1 * gamma_inc_moment(a, b, Bj + Tpt(1), αj) + - K2 * gamma_inc_moment(a, b, Bj + Tpt(2), αj) - ) - v_part_M += - -sgn * Aj * mfac * - ( - K0 * gamma_inc_moment(a, b, Bj + Tpt(3), αj) + - K1 * gamma_inc_moment(a, b, Bj + Tpt(4), αj) + - K2 * gamma_inc_moment(a, b, Bj + Tpt(5), αj) - ) + + # Compute rain PSD incomplete moments weighted by ice-liquid collision + # cross-section `K`, and sedimentation velocity difference `|vᵢ - vₗ|` + coeffs = collision_cross_section_ice_liquid_coeffs(rᵢ) + Iᵖ(a, b, p, α) = UU.unrolled_sum( + k * gamma_inc_moment(a, b, p + i - 1, α) for (i, k) in enumerate(coeffs) + ) + function flux(a, b, p) # ≡ ∫ₐᵇ K(Dᵢ, Dₗ) ⋅ (vᵢ(Dᵢ) - vₗ(Dₗ)) ⋅ n_r(Dₗ) dDₗ + s = v_i_at_Dᵢ * Iᵖ(a, b, p, λ) # vᵢ ⋅ ∫ₐᵇ K ⋅ n_r dDₗ + s -= UU.unrolled_mapreduce(+, ai, bi, ci) do aⱼ, bⱼ, cⱼ # - ∫ₐᵇ K ⋅ vₗ ⋅ n_r dDₗ + aⱼ * Iᵖ(a, b, p + bⱼ, λ + cⱼ) end - return (vi_part_N + v_part_N, vi_part_M + v_part_M) + return s end - lower = piece_NM(T, Tp, D_min, Dstar, +one(T)) - upper = piece_NM(T, Tp, Dstar, D_max, -one(T)) - return (N₀r * (lower[1] + upper[1]), N₀r * (lower[2] + upper[2])) + crossing(p) = flux(D_min, Dstar, p) - flux(Dstar, D_max, p) # sign flip at Dstar + mfac = ρw * CO.volume_sphere_D(one(FT)) # m_liq(D) = mfac Dₗ³ + return (N₀r * crossing(FT(0)), N₀r * mfac * crossing(FT(3))) # number: D⁰, mass: D³ end """ - get_liquid_integrals_rain_closed(psd_r, vel, n_r, ρₐ, L_r, N_r, state, - ∂ₜV, m_liq, ρ′_rim, bounds_r; quad) - -Closed-form rain inner integrals: a `liquid_integrals(Dᵢ)` returning -`(∂ₜN_col, ∂ₜM_col, ∂ₜB_col)` where N and M are the exact incomplete- -gamma closed form and B-rim falls back to the numerical 1-D `integrate` -(piecewise-rational `1/ρ′_rim`; see A3 §4 — deliberately not closed). -Dispatched only for `(RainParticlePDF_SB2006, Chen2022VelType)`. - -The Chen-2022 rain `v_l` callable is sourced from the **canonical** CM -API `CO.particle_terminal_velocity(vel.rain, ρₐ)` — the *same* -construction the numerical fallback uses (`volumetric_collision_rate_ -integrand`, `compute_local_rime_density`) — so the closed path tracks -any future change to CM's Chen law instead of forking from it. The -B-rim quadrature reuses the **passed-in `n_r` closure** (identical to -the numerical path's `∂ₜV·n_r·m_liq/ρ′_rim` integrand) rather than -re-deriving `N₀r e^{-D/D̄r}` by hand. - -`α0 = inv(Dr_mean) > 0` is an invariant of the `(SB2006, Chen)` -bundle (`Dr_mean > 0` for any valid SB2006 PSD, Chen rain `cⱼ > 0` ⇒ -`αⱼ = α0+cⱼ > 0`), so `gamma_inc_moment`'s `α≤0` NaN sentinel is -unreachable here. As a belt-and-suspenders matching the numerical -fallback's robustness (and the `iszero(N₀r)` guard above — one NaN -cell would otherwise abort a whole TRMM column), a degenerate -non-finite `(N,M)` is folded to `(0,0)` rather than NaN-poisoning the -outer integral. + get_liquid_integrals_rain_closed( + psd_r::RainParticlePDF_SB2006, vel::Chen2022VelType, + n_r, ρₐ, L_r, N_r, state, ∂ₜV, m_liq, ρ′_rim, bounds_r; quad + ) + +Returns a function `liquid_integrals(Dᵢ) -> (∂ₜN_col, ∂ₜM_col, ∂ₜB_col)` +where N and M are the exact incomplete-gamma closed form and +B_rim is computed by quadrature """ function get_liquid_integrals_rain_closed( psd_r::CMP.RainParticlePDF_SB2006, vel::CMP.Chen2022VelType, n_r, ρₐ, L_r, N_r, state, ∂ₜV, m_liq, ρ′_rim, bounds_r; quad, ) FT = eltype(state) - ρ_w = psd_r.ρw + ρw = psd_r.ρw (; N₀r, Dr_mean) = CM2.pdf_rain_parameters(psd_r, L_r / ρₐ, ρₐ, N_r) ai, bi, ci = CO.Chen2022_vel_coeffs(vel.rain, ρₐ) - # Canonical Chen-2022 rain v(D) — the exact same callable the - # numerical fallback builds (`v_liq = CO.particle_terminal_velocity( - # velocity_params.rain, ρₐ)`, see volumetric_collision_rate_integrand). v_l = CO.particle_terminal_velocity(vel.rain, ρₐ) v_i = ice_particle_terminal_velocity(vel, ρₐ, state) D_min, D_max = bounds_r @@ -504,21 +385,15 @@ function get_liquid_integrals_rain_closed( if iszero(N₀r) || !(D_max > D_min) return (zero(FT), zero(FT), zero(FT)) end - vi_Dᵢ = v_i(Dᵢ) + v_i_at_Dᵢ = v_i(Dᵢ) rᵢ = sqrt(ice_area(state, Dᵢ) / FT(π)) ∂ₜN_col, ∂ₜM_col = closed_rain_inner_NM( - FT(Dᵢ), vi_Dᵢ, v_l, rᵢ, ρ_w, ai, bi, ci, + FT(Dᵢ), v_i_at_Dᵢ, v_l, rᵢ, ρw, ai, bi, ci, D_min, D_max, N₀r, Dr_mean, ) - # Degrade like the numerical path (→0) instead of propagating a - # NaN through the outer 10-tuple if a degenerate state ever made - # the unreachable `α≤0` sentinel fire (see docstring / P1-1). if !(isfinite(∂ₜN_col) && isfinite(∂ₜM_col)) return (zero(FT), zero(FT), zero(FT)) end - # B-rim: piecewise-rational `1/ρ′_rim` — keep numerical (A3 §4). - # Reuse the passed-in `n_r` closure (identical integrand to the - # numerical `get_liquid_integrals` B-rim) — no hand re-derivation. ∂ₜB_col = integrate( D -> ∂ₜV(Dᵢ, D) * n_r(D) * m_liq(D) / ρ′_rim(Dᵢ, D), bounds_r, @@ -529,8 +404,6 @@ function get_liquid_integrals_rain_closed( return liquid_integrals end -# Dispatch: closed form for the SB2006-exp × Chen-2022 bundle, else the -# existing numerical path (behavior byte-unchanged for other bundles). _rain_inner_integrals( psd_r::CMP.RainParticlePDF_SB2006, vel::CMP.Chen2022VelType, n_r, ∂ₜV, m_liq, ρ′_rim, bounds_r, ρₐ, L_r, N_r, state; quad, @@ -661,10 +534,8 @@ function ∫liquid_ice_collisions( ∂ₜM_max = compute_max_freeze_rate(aps, tps, vel, ρₐ, T, state) # ∂ₜM_max(Dᵢ) cloud_integrals = get_liquid_integrals(n_c, ∂ₜV, m_liq, ρ′_rim, bounds_c; quad) # (∂ₜN_c_col, ∂ₜM_c_col, ∂ₜB_c_col) - # Rain inner: exact closed form for the (SB2006-exp PSD, Chen-2022) - # bundle (N,M via incomplete gamma; B-rim numerical). Numerical - # fallback for any other PSD/velocity type. Cloud inner + outer axis - # unchanged (compose with GL40). See A3 §4 / A4 Phase 3 / issue 003. + # Rain inner: exact closed form for the (SB2006-exp PSD, Chen-2022) pair + # Numerical fallback for any other PSD/velocity type. rain_integrals = _rain_inner_integrals( psd_r, vel, n_r, ∂ₜV, m_liq, ρ′_rim, bounds_r, ρₐ, L_r, N_r, state; quad, diff --git a/src/P3_size_distribution.jl b/src/P3_size_distribution.jl index 3e5c005c1..f74da48d0 100644 --- a/src/P3_size_distribution.jl +++ b/src/P3_size_distribution.jl @@ -69,28 +69,13 @@ Compute `log(Iᵏ)` where `Iᵏ` is the following integral: In log-space, this is: `- z log(λ) + logΓ(z) + log(q_D₁ - q_D₂)` - !!! note "Log/linear twin" - [`gamma_inc_moment`](@ref) (just below) is the **linear-space** - sibling of this function: it returns the same incomplete-gamma - partial moment `Γ(z)(q_D₁−q_D₂)/λ^z` *without* the `log`. The - closed-form rain-inner collision split (`closed_rain_inner_NM`, - `P3_processes.jl`) needs the linear form because it accumulates - **sign-alternating** addends (the `±sgn` split pieces and the - `vᵢ_part − v_part` subtraction) that cannot be done under a `log` - / `logsumexp`. The two are kept adjacent so the pair stays in - sync (P1-2). `Common.Chen2022_exponential_pdf` is the *full-domain* - `[0,∞)` 1M analogue (complete, not incomplete, gamma) — not a - drop-in for the closed form, which needs *partial* `[D₁,D₂]` - moments because of the velocity-crossover `D*` split. - +See also [`gamma_inc_moment`](@ref) """ function loggamma_inc_moment(D₁, D₂, μ, logλ, k = 0, scale = 1) FT = eltype(logλ) D₁ < D₂ || return log(FT(0)) # return log(0) if D₁ ≥ D₂ z = k + μ + 1 - # NOTE: We use `LogExpFunctions.xexpy(D, logλ)` to compute `λD = D * exp(logλ)`. - # When `logλ` is large, `exp(logλ) = Inf`, so the naive product `D * exp(logλ)` - # yields `0 * Inf = NaN` when `D = 0`. `xexpy` correctly returns `0` in that case. + # `λ⋅D ≡ xexpy(D, logλ) ≡ D * exp(logλ)` (numerically stable) (_, q_D₁) = SF.gamma_inc(z, LogExpFunctions.xexpy(D₁, logλ)) (_, q_D₂) = SF.gamma_inc(z, LogExpFunctions.xexpy(D₂, logλ)) return -z * logλ + SF.loggamma(z) + log(q_D₁ - q_D₂) + log(FT(scale)) @@ -99,47 +84,20 @@ end """ gamma_inc_moment(D₁, D₂, p, α) -`∫_{D₁}^{D₂} D^p e^{-α D} dD = α^{-(p+1)} Γ(p+1) [Q(p+1,αD₁) − Q(p+1,αD₂)]` +`∫_{D₁}^{D₂} D^p e^{-α D} dD = α^{-(p+1)} Γ(p+1) [Q(p+1,αD₁) - Q(p+1,αD₂)]` with `Q` the regularized upper incomplete gamma. -The **linear-space twin** of [`loggamma_inc_moment`](@ref) (same -`Γ(z)(q_{D₁}−q_{D₂})/λ^z`, `z = p+1` reduction, without the `log` and -without `loggamma_inc_moment`'s `logλ`/`xexpy` overflow guard, which is -unneeded here — the closed-form slopes `α0=1/D̄r`, `αⱼ=α0+cⱼ` are plain -positive scalars well inside range). Used by the closed-form rain-inner -collision split (`closed_rain_inner_NM`, `P3_processes.jl`), whose -**sign-alternating** addends cannot be accumulated under a `log` / -`logsumexp` — see the "Log/linear twin" note on -[`loggamma_inc_moment`](@ref). Keep the two adjacent and in sync (P1-2). - -`Common.Chen2022_exponential_pdf` is the *full-domain* `[0,∞)` 1M -analogue (complete gamma); it is **not** a drop-in here because the -velocity-crossover `D*` split needs *partial* `[D₁,D₂]` (incomplete) -moments. - -Returns `0` if `D₂ ≤ D₁`. Returns `NaN` if `α ≤ 0` (a mis-set slope); -for the dispatched `(SB2006, Chen)` bundle `α > 0` is an invariant -(`D̄r > 0`, Chen rain `cⱼ > 0`) so that branch is unreachable, and -`get_liquid_integrals_rain_closed` folds any non-finite result to `0` -as a defensive belt-and-suspenders (P1-1). `z = p+1 > 0` is likewise -guaranteed (Chen rain `Bⱼ > 0`, `bounds_r` lower > 0) so `SF.gamma(z)` -is always finite (P2-2). `p+1` is a Float64 constant — the `Dual` -enters only via `gamma_inc`'s **2nd** argument (AD-safe, A2). +Returns `0` if `D₂ ≤ D₁`, and `NaN` if `α ≤ 0` + +See also [`loggamma_inc_moment`](@ref) """ @inline function gamma_inc_moment(D₁, D₂, p, α) - # AD-generic working type: promotes Float64 (e.g. a Float64 D₁/D₂ - # quadrature endpoint) and a `Dual` slope `α` (Dual via state/PSD - # params) to the Dual so the partial moment differentiates. `p` is - # a plain Float64/Int exponent (z = p+1 stays a non-Dual constant — - # the A2-safe `gamma_inc` 2nd-arg-only pattern; never a Dual 1st arg). - T = float(promote_type(typeof(D₁), typeof(D₂), typeof(α))) - D₂ > D₁ || return zero(T) - α > 0 || return T(NaN) - z = p + 1 # > 0: Chen rain Bⱼ>0 & bounds_r lower>0 ⇒ SF.gamma(z) finite + FT = float(promote_type(typeof(D₁), typeof(D₂), typeof(α))) + D₂ > D₁ || return zero(FT) + α > 0 || return FT(NaN) + z = p + 1 (_, q1) = SF.gamma_inc(z, α * D₁) (_, q2) = SF.gamma_inc(z, α * D₂) - # Float64 * Dual / Dual auto-promotes to the Dual `T` (the `zero(T)` - # / `T(NaN)` early returns are why `T` is computed explicitly). return SF.gamma(z) * (q1 - q2) / α^z end @@ -275,33 +233,7 @@ where `m(D)` is the mass of a particle at diameter `D` (see [`ice_mass`](@ref)). # Arguments - `state`: The [`P3State`](@ref) -- `logλ_guess`: Optional warm-start seed. When `logλ_guess` is finite and - strictly inside `[logλ_min, logλ_max]`, one extra evaluation of the - residual at the guess is used to **halve the Brent bracket** by sign: - whichever endpoint sits on the same side of the root as `f(logλ_guess)` - is replaced by the guess, and Brent then converges from the tighter - bracket. If the guess is `nothing` / non-finite / out-of-bracket, or - evaluates to a non-finite residual, the solver falls back to the full - bracket with no behaviour change. Default: `nothing`. - - **Caveat — monotonicity depends on the `μ(λ)` parameterization.** - - - Under [`CMP.SlopeConstant`](@ref) (μ fixed, independent of λ), - `log(L/N)(logλ)` is a strictly decreasing function of `logλ` — - larger slope ⇒ smaller mean mass. Brent finds the unique root on - any valid bracket and the sign-based warm-start narrowing is - exact. - - Under [`CMP.SlopePowerLaw`](@ref) (μ = clamp(a·λ^b - c, 0, μ_max) - — piecewise flat, rising, flat), the residual is *not* globally - monotonic and the same target `L/N` can have several roots (see - `docs/src/plots/P3SlopeParameterizations.jl` for a three-root - example). In that regime the halved bracket still contains at - least one root, but it may not be the root closest to the guess; - the warm-start is heuristic rather than exact. For smooth cell - evolution (step-to-step `logλ` continuity) this is usually fine - because the guess and the current root stay on the same branch, - but callers that need all roots should use - `get_distribution_logλ_all_solutions`. +- `logλ_guess`: Optional initial guess - `logλ_min`: The minimum value of the search bounds [log(1/m)], default is `log(1e1)` - `logλ_max`: The maximum value of the search bounds [log(1/m)], default is `log(1e7)` """ @@ -314,26 +246,11 @@ function get_distribution_logλ(state, logλ_guess = nothing, logλ_min = 2, log target_log_LdN = log(ρq_ice) - log(ρn_ice) shape_problem(logλ) = logLdivN(state, logλ) - target_log_LdN - - # Cold bracket. Brent converges to *some* root inside a valid bracket. - # Whether that root is unique over the full bracket depends on the - # `μ(λ)` parameterization: monotone for `SlopeConstant`, possibly - # multi-valued for `SlopePowerLaw` (see docstring caveat). If either - # endpoint residual is non-finite or the endpoints share a sign, the - # target `L/N` lies outside the representable PSD range — saturate at - # the nearer endpoint rather than failing. lo, hi = FT(logλ_min), FT(logλ_max) f_lo, f_hi = shape_problem(lo), shape_problem(hi) if !isfinite(f_lo) || !isfinite(f_hi) || f_lo * f_hi > 0 return abs(f_lo) ≤ abs(f_hi) ? lo : hi end - - # Optional hot-start: one probe at the prior step's `logλ` halves the - # bracket by sign, and Brent takes over from there. Any further - # narrowing is Brent's job — its interpolating iteration is more - # efficient per f-eval than any fixed-offset probe we could add here. - # Non-monotonic states may have the halved bracket land on a different - # root than the guess; see the docstring caveat. (lo, f_lo, hi, f_hi) = _narrow_bracket(shape_problem, lo, f_lo, hi, f_hi, logλ_guess) @@ -361,13 +278,9 @@ function get_distribution_logλ_from_prognostic( return get_distribution_logλ(state, args...) end -# One-probe narrowing of a valid bracket `[lo, hi]` using a probe point `p`. -# No-op if `p` is not usable (nothing / non-finite / outside the bracket) or -# if `f(p)` is non-finite. Otherwise, replace whichever endpoint sits on the -# same side of the root as `p` — a single sign check on `f_lo * f_p`. @inline _narrow_bracket(_sp, lo, f_lo, hi, f_hi, ::Nothing) = (lo, f_lo, hi, f_hi) @inline function _narrow_bracket(shape_problem, lo, f_lo, hi, f_hi, p::Real) - p_ = oftype(lo, p) # match FT of the bracket for type stability + p_ = oftype(lo, p) (isfinite(p_) && lo < p_ < hi) || return (lo, f_lo, hi, f_hi) f_p = shape_problem(p_) isfinite(f_p) || return (lo, f_lo, hi, f_hi) diff --git a/src/P3_terminal_velocity.jl b/src/P3_terminal_velocity.jl index 791f80f7e..9e120d6c6 100644 --- a/src/P3_terminal_velocity.jl +++ b/src/P3_terminal_velocity.jl @@ -5,13 +5,6 @@ Returns a single-argument function `v_term(D)` that gives the Chen 2022 terminal velocity of an ice particle of maximum dimension `D`. -The size-independent coefficient work (`Chen2022_vel_coeffs` and -monodisperse-PDF construction for both small- and large-ice regimes) is -done once at call time; the returned closure only does the per-`D` -evaluation and — if `use_aspect_ratio = true` — the aspect-ratio -correction `cbrt(ϕᵢ(state, D))`. `ϕᵢ(state, D)` is O(1) because `state` -caches the regime thresholds. - # Arguments - `velocity_params`: A [`CMP.Chen2022VelType`](@ref) - `ρₐ`: Air density [kg/m³] @@ -26,19 +19,9 @@ caches the regime thresholds. FT = typeof(ρₐ) (; small_ice, large_ice) = velocity_params D_cutoff = small_ice.cutoff - ρᵢ = FT(916.7) + ρᵢ = FT(916.7) # TODO: Use parameter v_term_small = CO.particle_terminal_velocity(small_ice, ρₐ, ρᵢ) v_term_large = CO.particle_terminal_velocity(large_ice, ρₐ, ρᵢ) - - # Return ONE closure (not one of two closures chosen by `use_aspect_ratio`). - # A `use_aspect_ratio ? closureA : closureB` return makes the closure's - # *type* depend on a runtime value and adds a closure-nesting layer; when - # that closure is composed into an integrand (`n(D) * v_term(D)`) Julia 1.10 - # inference exceeds its closure-recursion limit and falls back to `::Any`, - # which poisons `integrate` (runtime dispatch / boxing on CPU, InvalidIRError - # on GPU). Folding the aspect-ratio choice into a scalar branch *inside* the - # single closure keeps the return type concrete (`use_aspect_ratio` is - # constant-propagated through this `@inline`). v_term(D) = let vₜ = D <= D_cutoff ? v_term_small(D) : v_term_large(D) use_aspect_ratio ? cbrt(ϕᵢ(state, D)) * vₜ : vₜ diff --git a/src/parameters/IceNucleation.jl b/src/parameters/IceNucleation.jl index 30cc6bacf..ca070dc15 100644 --- a/src/parameters/IceNucleation.jl +++ b/src/parameters/IceNucleation.jl @@ -127,7 +127,7 @@ Returns the volumetric freezing rate in SI units [m⁻³ s⁻¹]. The stored `het_B` is in [cm⁻³ s⁻¹]; the factor 10⁶ converts cm³ → m³. """ @kwdef struct RainFreezing{FT} <: ParametersType - "empirical parameter [°C⁻¹]" + "empirical parameter [°K⁻¹]" het_a::FT "water-type dependent parameter [cm⁻³ s⁻¹]" het_B::FT @@ -136,7 +136,7 @@ end function RainFreezing(td::CP.ParamDict) name_map = (; :BarklieGokhale1959_a_parameter => :het_a, - :BarklieGokhale1959_B_parameter => :het_B, + :BarklieGokhale1959_B_parameter => :het_B, # TODO: Fix units in CP, then here ) parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics") return RainFreezing(; parameters...) @@ -203,31 +203,7 @@ end # F23 INP-activation memory models # --------------------------------------------------------------------------- -export AbstractINPDepletion, NIceProxyDepletion - -""" - AbstractINPDepletion - -Abstract type for the model of how F23 INP-activation budgets are -"depleted" within an air parcel. Choose a concrete subtype to control -the value subtracted from the F23 INPC target in the F23 deposition and -immersion-cap rates: - -``` -∂ₜn_frz = max(0, INPC(T)/ρ - n_active) / τ_act -``` - -where `n_active` is supplied by the host. The depletion model also -carries `τ_act`, the F23 activation relaxation timescale, so the host -doesn't need to wire that knob separately. The model selects how the -host sources `n_active`: - -- [`NIceProxyDepletion`](@ref): use existing `n_ice` (zero memory). - -(A prognostic activation-memory model with a finite decay timescale lived -here previously and is deferred to a follow-up PR.) -""" -abstract type AbstractINPDepletion end +export NIceProxyDepletion """ NIceProxyDepletion{FT} @@ -246,7 +222,7 @@ clean air below it and the F23 channel artificially shuts off. # Fields $(DocStringExtensions.FIELDS) """ -struct NIceProxyDepletion{FT} <: AbstractINPDepletion +struct NIceProxyDepletion{FT} "F23 activation relaxation timescale [s] (default `300`)" τ_act::FT end diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl index a61e9015b..7477357ca 100644 --- a/test/gpu_tests.jl +++ b/test/gpu_tests.jl @@ -1200,24 +1200,12 @@ function test_gpu(FT) if VERSION < v"1.12" # The collision path exceeds Julia <= 1.11's inference depth, so this # kernel does not compile there (InvalidIRError: dynamic dispatch). - # 1.12 inference resolves it. TODO: ungate once GPU CI is >= 1.12. + # 1.12 inference resolves it. TODO: fix once GPU CI is >= 1.12. TT.@test_broken VERSION >= v"1.12" else kernel!( - mp_2m_p3, - tps, - output, - ρ, - T, - q_tot, - q_lcl, - n_lcl, - q_rai, - n_rai, - q_ice, - n_ice, - q_rim, - b_rim; + mp_2m_p3, tps, output, ρ, T, q_tot, + q_lcl, n_lcl, q_rai, n_rai, q_ice, n_ice, q_rim, b_rim; ndrange, ) TT.@test allequal(Array(output)) diff --git a/test/heterogeneous_ice_nucleation_tests.jl b/test/heterogeneous_ice_nucleation_tests.jl index 910a2e510..ecba515b8 100644 --- a/test/heterogeneous_ice_nucleation_tests.jl +++ b/test/heterogeneous_ice_nucleation_tests.jl @@ -250,25 +250,25 @@ function test_heterogeneous_ice_nucleation(FT) τ = FT(300) # Above T_freeze: zero - r_warm = CMI_het.f23_immersion_limit_rate( + r_warm = CMI_het.immersion_limit_rate( ip_frostenberg, T_freeze + FT(0.1), ρ; τ, ) TT.@test r_warm.∂ₜn_frz == FT(0) # Cold T: rate is positive and matches INPC/(ρ·τ) T_cold = T_freeze - FT(20) - r_cold = CMI_het.f23_immersion_limit_rate(ip_frostenberg, T_cold, ρ; τ) + r_cold = CMI_het.immersion_limit_rate(ip_frostenberg, T_cold, ρ; τ) INPC_expected = exp(CMI_het.INP_concentration_mean(ip_frostenberg, T_cold)) / ρ / τ TT.@test r_cold.∂ₜn_frz ≈ INPC_expected rtol = sqrt(eps(FT)) # Colder ⇒ larger rate - r_colder = CMI_het.f23_immersion_limit_rate( + r_colder = CMI_het.immersion_limit_rate( ip_frostenberg, T_freeze - FT(30), ρ; τ, ) TT.@test r_colder.∂ₜn_frz > r_cold.∂ₜn_frz # log_inpc_shift > 0 ⇒ larger rate - r_shifted = CMI_het.f23_immersion_limit_rate( + r_shifted = CMI_het.immersion_limit_rate( ip_frostenberg, T_cold, ρ; τ, inpc_log_shift = FT(1), ) TT.@test r_shifted.∂ₜn_frz ≈ r_cold.∂ₜn_frz * exp(FT(1)) rtol = sqrt(eps(FT)) @@ -295,7 +295,7 @@ function test_heterogeneous_ice_nucleation(FT) # Strict defaults are picked up automatically when not overridden. T_test = T_freeze - FT(20) # below the -15 °C strict gate q_sat_test = TDI.saturation_vapor_specific_content_over_ice(tps, T_test, ρ) - r_default = CMI_het.f23_deposition_rate( + r_default = CMI_het.deposition_rate( ip_frostenberg, tps, T_test, ρ, 2 * q_sat_test, FT(0), FT(0), n_ice; m_nuc, τ_act, ) @@ -303,7 +303,7 @@ function test_heterogeneous_ice_nucleation(FT) # And they're closed at T = -10 °C (above the -15 °C gate) T_warm_test = T_freeze - FT(10) q_sat_warm_test = TDI.saturation_vapor_specific_content_over_ice(tps, T_warm_test, ρ) - r_default_closed = CMI_het.f23_deposition_rate( + r_default_closed = CMI_het.deposition_rate( ip_frostenberg, tps, T_warm_test, ρ, 2 * q_sat_warm_test, FT(0), FT(0), n_ice; m_nuc, τ_act, ) @@ -316,7 +316,7 @@ function test_heterogeneous_ice_nucleation(FT) T_cold = T_freeze - FT(20) q_sat_ice = TDI.saturation_vapor_specific_content_over_ice(tps, T_cold, ρ) q_vap_super = 2 * q_sat_ice # S_i ≈ 1.0 - r_active = CMI_het.f23_deposition_rate( + r_active = CMI_het.deposition_rate( ip_frostenberg, tps, T_cold, ρ, q_vap_super, FT(0), FT(0), n_ice; m_nuc, T_thresh = T_thresh_default, S_i_thresh = S_i_thresh_default, τ_act, @@ -329,7 +329,7 @@ function test_heterogeneous_ice_nucleation(FT) # n_ice = INPC ⇒ depleted to zero ⇒ both n and q rates vanish INPC_at_T = exp(CMI_het.INP_concentration_mean(ip_frostenberg, T_cold)) / ρ - r_depleted = CMI_het.f23_deposition_rate( + r_depleted = CMI_het.deposition_rate( ip_frostenberg, tps, T_cold, ρ, q_vap_super, FT(0), FT(0), INPC_at_T; m_nuc, T_thresh = T_thresh_default, S_i_thresh = S_i_thresh_default, τ_act, @@ -343,7 +343,7 @@ function test_heterogeneous_ice_nucleation(FT) T_warm = T_freeze - FT(10) q_sat_warm = TDI.saturation_vapor_specific_content_over_ice(tps, T_warm, ρ) q_vap_super_warm = 2 * q_sat_warm - r_T_gate = CMI_het.f23_deposition_rate( + r_T_gate = CMI_het.deposition_rate( ip_frostenberg, tps, T_warm, ρ, q_vap_super_warm, FT(0), FT(0), n_ice; m_nuc, T_thresh = T_thresh_default, S_i_thresh = S_i_thresh_default, τ_act, @@ -353,7 +353,7 @@ function test_heterogeneous_ice_nucleation(FT) # Subsaturated wrt ice ⇒ both rates zero (S_i gate closed, and # the vapor cap independently zeros ∂ₜq_frz). - r_sub = CMI_het.f23_deposition_rate( + r_sub = CMI_het.deposition_rate( ip_frostenberg, tps, T_cold, ρ, FT(0.5) * q_sat_ice, FT(0), FT(0), n_ice; m_nuc, T_thresh = T_thresh_default, S_i_thresh = S_i_thresh_default, τ_act, @@ -362,7 +362,7 @@ function test_heterogeneous_ice_nucleation(FT) TT.@test r_sub.∂ₜq_frz == FT(0) # Tunable thresholds: warming the T_thresh opens the warm-side gate - r_relaxed = CMI_het.f23_deposition_rate( + r_relaxed = CMI_het.deposition_rate( ip_frostenberg, tps, T_warm, ρ, q_vap_super_warm, FT(0), FT(0), n_ice; m_nuc, T_thresh = T_freeze - FT(5), S_i_thresh = S_i_thresh_default, τ_act, @@ -371,7 +371,7 @@ function test_heterogeneous_ice_nucleation(FT) TT.@test r_relaxed.∂ₜq_frz > FT(0) # Permissive thresholds (used by BMT) ⇒ always passes the gate - r_permissive = CMI_het.f23_deposition_rate( + r_permissive = CMI_het.deposition_rate( ip_frostenberg, tps, T_warm, ρ, q_vap_super_warm, FT(0), FT(0), n_ice; m_nuc, T_thresh = FT(2000), S_i_thresh = FT(-2), τ_act, ) @@ -385,7 +385,7 @@ function test_heterogeneous_ice_nucleation(FT) q_sat_super = TDI.saturation_vapor_specific_content_over_ice(tps, T_super, ρ) q_vap_tiny_excess = q_sat_super * FT(1.001) # S_i = 0.001 q_excess_tiny = q_vap_tiny_excess - q_sat_super - r_capped = CMI_het.f23_deposition_rate( + r_capped = CMI_het.deposition_rate( ip_frostenberg, tps, T_super, ρ, q_vap_tiny_excess, FT(0), FT(0), FT(0); m_nuc = FT(1e3), # absurd m_nuc forces vapor cap to bind T_thresh = FT(2000), S_i_thresh = FT(-2), τ_act, diff --git a/test/p3_tests.jl b/test/p3_tests.jl index 372419933..326dbe1b3 100644 --- a/test/p3_tests.jl +++ b/test/p3_tests.jl @@ -38,12 +38,6 @@ function test_p3_state_creation(FT) # Test thresholds for rimed state (; D_th, D_gr, D_cr) = state_rimed @test D_th < D_gr < D_cr - - # Note: `P3State` no longer validates `F_rim ∈ [0, 1)` or - # `ρ_rim ∈ [0, ρ_l]` at construction. Production input regularisation - # lives in `state_from_prognostic` (clamps via `min(...)`); test code - # is expected to pass valid `(F_rim, ρ_rim)` directly. See the - # `P3State` docstring. end end @@ -61,10 +55,6 @@ function test_thresholds_solver(FT) ρ_rim_good = (FT(200), FT(400), FT(800)) # representative ρ_rim values F_rim_good = (FT(0.5), FT(0.8), FT(0.95)) # representative F_rim values - # Note: `P3State` no longer asserts `F_rim ∈ [0, 1)` or - # `ρ_rim ∈ [0, ρ_l]`. Domain enforcement happens upstream in - # `state_from_prognostic` via clamps. Tests pass valid inputs. - # Test if the P3 scheme solution satisifies the conditions # from eqs. 14-17 in Morrison and Milbrandt 2015 function get_ρ_d_paper((; α_va, β_va)::CMP.MassPowerLaw; D_cr, D_gr) @@ -658,9 +648,6 @@ function test_p3_melting(FT) T_vwarm = FT(273.15 + 0.1) rate = P3.ice_melt(vel, aps, tps, T_vwarm, ρₐ, state, logλ) - - # Uncapped melt rate (dt/availability limiting removed): the raw rate - # exceeds the available L/N. Reference values are output from the code. if FT == Float64 ref_vwarm_dNdt = FT(1.7186681756049155e6) ref_vwarm_dLdt = FT(8.593340878024577e-4) @@ -912,25 +899,23 @@ function test_p3_ice_self_collection(FT) end end -# Closed-form rain-inner liquid-ice collision integral (A3 / issue 003). -# The exact incomplete-gamma reduction of the rain inner {N,M} for the -# (RainParticlePDF_SB2006, Chen2022VelType) bundle, validated vs a -# high-order numerical reference; plus the dispatch fallback contract. function test_p3_closed_form_rain_inner(FT) - @testset "P3 closed-form rain inner (N,M)" begin + # The closed form is the exact analytic reduction of the rain inner integral. + # Testset checks: + # - N and M matches an adaptive (QuadGK) reference numerical quadrature + # - correctness of the rime volume quadrature + # - smoke test for collision cross section + @testset "P3 closed-form rain inner (N, M, B) + cross-section coeffs" begin params = CMP.ParametersP3(FT) vel = CMP.Chen2022VelType(FT) - sb = CMP.SB2006(FT) - psd_r = sb.pdf_r + psd_r = CMP.SB2006(FT).pdf_r ρₐ = FT(1) ρ_w = psd_r.ρw p = eltype(params)(1e-5) m_liq(D) = ρ_w * FT(π) / 6 * D^3 - # rtol ≫ the measured ~5.4e-6 (N) / 2.0e-6 (M) closed-vs-CG(1024) - # margin (and the FT=Float32 round-off floor), so the test is a - # robust regression guard, not a brittle exact-digits check. - rtol = FT == Float64 ? FT(1e-3) : FT(2e-2) - ref_n = FT == Float64 ? 1024 : 256 + rtol = FT == Float64 ? FT(1e-10) : FT(1e-3) + qrtol = FT == Float64 ? FT(1e-12) : FT(1e-7) + v_l = CO.particle_terminal_velocity(vel.rain, ρₐ) for (L_ice, N_ice, F_rim, ρ_rim) in ( (1e-3, 1e6, 0.5, 500), (1e-2, 1e8, 0.95, 800), @@ -938,105 +923,50 @@ function test_p3_closed_form_rain_inner(FT) ), (L_r, N_r) in ((1e-6, 1e4), (1e-4, 1e3), (2e-3, 5e2)) - state = P3.P3State( - params, FT(L_ice), FT(N_ice), FT(F_rim), FT(ρ_rim), - ) + state = P3.P3State(params, FT(L_ice), FT(N_ice), FT(F_rim), FT(ρ_rim)) n_r = DT.size_distribution(psd_r, FT(L_r) / ρₐ, ρₐ, FT(N_r)) ∂ₜV = P3.volumetric_collision_rate_integrand(vel, ρₐ, state) ρ′_rim = P3.compute_local_rime_density(vel, ρₐ, FT(270), state) - bnds = CM2.get_size_distribution_bounds( - psd_r, FT(L_r) / ρₐ, ρₐ, FT(N_r), p, - ) - bnds[2] > bnds[1] || continue + D_min, D_max = bnds = CM2.get_size_distribution_bounds(psd_r, FT(L_r) / ρₐ, ρₐ, FT(N_r), p) + D_max > D_min || continue rc = P3.get_liquid_integrals_rain_closed( psd_r, vel, n_r, ρₐ, FT(L_r), FT(N_r), state, ∂ₜV, m_liq, ρ′_rim, bnds; quad = P3.ChebyshevGauss(40), ) - rn = P3.get_liquid_integrals( - n_r, ∂ₜV, m_liq, ρ′_rim, bnds; - quad = P3.ChebyshevGauss(ref_n), + rn = P3.get_liquid_integrals( # numerical fallback + n_r, ∂ₜV, m_liq, ρ′_rim, bnds; quad = P3.ChebyshevGauss(40), ) + v_i = P3.ice_particle_terminal_velocity(vel, ρₐ, state) for Dᵢ in FT.(10 .^ range(-5, -2; length = 5)) - Nc, Mc, _ = rc(Dᵢ) - Nr, Mr, _ = rn(Dᵢ) - @test isapprox(Nc, Nr; rtol) - @test isapprox(Mc, Mr; rtol) - end - end - end + vi = v_i(Dᵢ) + Dstar = P3.crossover_diameter(vi, v_l, D_min, D_max) - # P1-test gap 1: lock the closed-vs-numerical-fallback contract on - # the *typed* bundle itself. The B-rim output now reuses the same - # `n_r` closure the numerical path uses (P1-5), so on the typed - # `(SB2006,Chen)` bundle the closed path's B-rim must equal the - # numerical path's B-rim *exactly* at matched quadrature, and N/M - # must agree to the closed-form accuracy margin. This proves the - # P0-1 canonical-Chen-API swap is numerically equivalent (same - # physical law) and that B-rim is a verbatim dedup, not a re-roll. - @testset "P3 closed vs numerical fallback (typed bundle)" begin - params = CMP.ParametersP3(FT) - vel = CMP.Chen2022VelType(FT) - psd_r = CMP.SB2006(FT).pdf_r - ρₐ = FT(1) - ρ_w = psd_r.ρw - p = eltype(params)(1e-5) - m_liq(D) = ρ_w * FT(π) / 6 * D^3 - rtolNM = FT == Float64 ? FT(1e-3) : FT(2e-2) - ref_n = FT == Float64 ? 1024 : 256 - for (L_ice, N_ice, F_rim, ρ_rim) in - ((1e-3, 1e6, 0.5, 500), (1e-2, 1e8, 0.95, 800)), - (L_r, N_r) in ((1e-4, 1e3), (2e-3, 5e2)) - - state = - P3.P3State(params, FT(L_ice), FT(N_ice), FT(F_rim), FT(ρ_rim)) - n_r = DT.size_distribution(psd_r, FT(L_r) / ρₐ, ρₐ, FT(N_r)) - ∂ₜV = P3.volumetric_collision_rate_integrand(vel, ρₐ, state) - ρ′_rim = P3.compute_local_rime_density(vel, ρₐ, FT(270), state) - bnds = CM2.get_size_distribution_bounds( - psd_r, FT(L_r) / ρₐ, ρₐ, FT(N_r), p, - ) - bnds[2] > bnds[1] || continue - # B-rim: closed path reuses `n_r` → must match the numerical - # path's B-rim *exactly* at the SAME quadrature order. - rc = P3.get_liquid_integrals_rain_closed( - psd_r, vel, n_r, ρₐ, FT(L_r), FT(N_r), state, ∂ₜV, - m_liq, ρ′_rim, bnds; quad = P3.ChebyshevGauss(40), - ) - rn_match = P3.get_liquid_integrals( - n_r, ∂ₜV, m_liq, ρ′_rim, bnds; - quad = P3.ChebyshevGauss(40), - ) - rn_ref = P3.get_liquid_integrals( - n_r, ∂ₜV, m_liq, ρ′_rim, bnds; - quad = P3.ChebyshevGauss(ref_n), - ) - for Dᵢ in FT.(10 .^ range(-5, -2; length = 4)) + # N, M: closed form vs adaptive reference Nc, Mc, Bc = rc(Dᵢ) - _, _, Bm = rn_match(Dᵢ) - Nr, Mr, _ = rn_ref(Dᵢ) - # B-rim is the same integrand & quadrature ⇒ identical. - @test isapprox(Bc, Bm; rtol = sqrt(eps(FT))) - # N,M: closed form vs high-order numerical reference. - @test isapprox(Nc, Nr; rtol = rtolNM) - @test isapprox(Mc, Mr; rtol = rtolNM) + σ(D) = P3.collision_cross_section_ice_liquid(state, Dᵢ, D) + gN(D) = σ(D) * abs(vi - v_l(D)) * n_r(D) + Nref = QGK.quadgk(gN, D_min, Dstar, D_max; rtol = qrtol)[1] + Mref = QGK.quadgk(D -> gN(D) * m_liq(D), D_min, Dstar, D_max; rtol = qrtol)[1] + @test isapprox(Nc, Nref; rtol) + @test isapprox(Mc, Mref; rtol) + + # B: closed quadrature vs numerical fallback at the same order + @test isapprox(Bc, rn(Dᵢ)[3]; rtol = sqrt(eps(FT))) + + # smoke test: collision cross section has form: `π(rᵢ + Dₗ/2)²`, with rᵢ derived from `P3.ice_area` + rᵢ = sqrt(P3.ice_area(state, Dᵢ) / FT(π)) + K = P3.collision_cross_section_ice_liquid_coeffs(state, Dᵢ) + @test K == P3.collision_cross_section_ice_liquid_coeffs(rᵢ) + @test evalpoly(Dstar, K) ≈ FT(π) * (rᵢ + Dstar / 2)^2 + @test evalpoly(D_max, K) ≈ FT(π) * (rᵢ + D_max / 2)^2 end end end - # P1-test gap 2: in-suite ForwardDiff AD smoke. The headline - # justification of the closed form (vs the numerical `abs` path) is - # AD-cleanliness under the A2 regime: in the implicit 2M+P3 - # Jacobian the `Dᵢ` quadrature node is a Float64 constant and the - # `Dual` enters only through the *state-dependent* physical inputs - # (`v̄ᵢ`, `rᵢ`, `D̄r`, `N₀r`) — `gamma_inc(z, α·D)` with `z`=Float64 - # const, Dual in the 2nd arg only (A2-safe). The D* velocity-split - # removes the `|v̄ᵢ−v_l|` kink, so `closed_rain_inner_NM` is C¹ in - # each of these inputs. We now `import ForwardDiff` (added to - # `test/Project.toml`) and check `FD.derivative` is finite and - # matches a central finite difference for each differentiated - # input — promoting the previously harness-only AD evidence - # (`A3_rain_inner_ad_smoothness.jl`) into the suite. - @testset "P3 closed-form ForwardDiff AD smoke (v̄ᵢ,rᵢ,D̄r,N₀r)" begin + # AD smoke test for the closed form. + # find D where ice/liquid sedimentation velocities are equal, + # separate the integrals, then check that closed form is correct. + @testset "P3 closed-form ForwardDiff AD smoke (v_i,r_i,Dr,N₀r)" begin vel = CMP.Chen2022VelType(FT) psd_r = CMP.SB2006(FT).pdf_r ρₐ = FT(1) @@ -1046,33 +976,19 @@ function test_p3_closed_form_rain_inner(FT) L_r, N_r = FT(1e-4), FT(1e3) (; N₀r, Dr_mean) = CM2.pdf_rain_parameters(psd_r, L_r / ρₐ, ρₐ, N_r) - Dᵢ0 = FT(3e-3) # Float64 quadrature node (real path) + Dᵢ0 = FT(3e-3) rᵢ0 = FT(1e-3) vi0 = FT(3.0) D_min, D_max = FT(1e-5), FT(5e-3) - # The headline AD-cleanliness claim — finite derivative, no - # `gamma_inc`-1st-arg-Dual MethodError, mixed Float64/Dual - # accepted — is asserted for BOTH FT. The *tight* FD-agreement - # is checked only for Float64: several M-moment derivatives are - # O(1e-5) values formed from cancelling O(1) incomplete-gamma - # terms, and a central FD of the bisection root, both of which - # are at the Float32 round-off floor (≈1e-7 rel) — not a defect - # of the closed form (Float64 agrees to ≤1e-11). The Float64 - # numerical check is the meaningful AD↔FD regression guard; - # Float32 only guards "AD runs cleanly and is finite". - check_fd = FT == Float64 rtolAD = FT(1e-4) - # N (out index 1) and M (out index 2) vs each differentiated - # input; `closed_rain_inner_NM(Dᵢ, v̄ᵢ, v_l, rᵢ, ρ_w, ai,bi,ci, - # D_min, D_max, N₀r, D̄r)`. cases = ( - ("v̄ᵢ", vi0, + ("v_i", vi0, x -> P3.closed_rain_inner_NM(Dᵢ0, x, v_l, rᵢ0, ρ_w, ai, bi, ci, D_min, D_max, N₀r, Dr_mean)), - ("rᵢ", rᵢ0, + ("r_i", rᵢ0, x -> P3.closed_rain_inner_NM(Dᵢ0, vi0, v_l, x, ρ_w, ai, bi, ci, D_min, D_max, N₀r, Dr_mean)), - ("D̄r", Dr_mean, + ("Dr", Dr_mean, x -> P3.closed_rain_inner_NM(Dᵢ0, vi0, v_l, rᵢ0, ρ_w, ai, bi, ci, D_min, D_max, N₀r, x)), ("N₀r", N₀r, @@ -1082,28 +998,17 @@ function test_p3_closed_form_rain_inner(FT) for (_, x0, g) in cases, idx in (1, 2) f(x) = g(x)[idx] d_ad = FD.derivative(f, x0) - @test isfinite(d_ad) # AD-clean (both FT) - if check_fd + @test isfinite(d_ad) # check that AD works + if FT == Float64 h = max(abs(x0) * FT(1e-5), FT(1e-12)) d_fd = (f(x0 + h) - f(x0 - h)) / (2h) - @test isapprox(d_ad, d_fd; rtol = rtolAD, - atol = 100 * eps(FT) * - max(abs(d_ad), abs(d_fd), one(FT))) + @test isapprox(d_ad, d_fd; + rtol = rtolAD, + atol = 100 * eps(FT) * max(abs(d_ad), abs(d_fd), one(FT)), + ) end end - # D* root differentiability. The fixed-iteration bisection is - # (by construction) AD-insensitive to its `v_target` — its - # branch comparisons discard the Dual, so FD-through-bisection - # returns ∂D*/∂v̄ᵢ = 0. That is the *correct* "frozen-D*-per-step" - # behaviour (A2: D* is frozen exactly like logλ). It is SAFE for - # the closed form because the split integrand is **continuous at - # D*** (v̄ᵢ−v_l(D*)=0 there) ⇒ the Leibniz boundary term - # `integrand(D*)·∂D*` vanishes regardless of ∂D* — which is why - # the `v̄ᵢ` cases above already match FD to ~1e-12 with a frozen - # root. Confirm both: (a) FD through the bisection is 0/finite - # (frozen-root), and (b) the elementary IFT the differentiated - # path *would* use, ∂D*/∂v̄ᵢ = 1/v_l′(D*), is correct (a central - # FD of the root vs the AD-clean elementary `v_l′`). + # Check that the crossover diameter is differentiable v_tgt = (v_l(D_min) + v_l(D_max)) / 2 dDstar_bisect = FD.derivative( vt -> P3.crossover_diameter(vt, v_l, D_min, D_max), v_tgt, @@ -1112,11 +1017,8 @@ function test_p3_closed_form_rain_inner(FT) Dstar = P3.crossover_diameter(v_tgt, v_l, D_min, D_max) vlp = FD.derivative(v_l, Dstar) # v_l′ elementary, AD-clean @test isfinite(vlp) && vlp > 0 - # IFT correctness via a central FD of the root — Float64 only: - # the bisection root is only accurate to ~eps(FT), so a FD over - # it is at the Float32 round-off floor (≈7% rel) while Float64 - # confirms ∂D*/∂v̄ᵢ = 1/v_l′(D*) to <1e-3. - if check_fd + # Check that it matches a numerical derivative + if FT == Float64 hv = abs(v_tgt) * FT(1e-6) dDstar_fd = ( @@ -1128,17 +1030,13 @@ function test_p3_closed_form_rain_inner(FT) end end - # P1-test gap 3: D*-at-bracket-end edges. Drive `crossover_diameter` - # past the rain velocity band so it must return D_min (v̄ᵢ below the - # whole band ⇒ lower piece empty) and D_max (v̄ᵢ above ⇒ upper piece - # empty). Result must stay finite and the absent-crossover case must - # collapse one split piece to zero (exercised via the public closed - # path with synthetic `v̄ᵢ` targets through `closed_rain_inner_NM`). - @testset "P3 closed-form D*-at-bracket-end (v̄ᵢ outside band)" begin + # Check edge cases (e.g. ice velocity > liquid velocity for all sizes) + # still returns the correct result + @testset "P3 closed-form D*-at-bracket-end (v_i outside band)" begin vel = CMP.Chen2022VelType(FT) v_l = CO.particle_terminal_velocity(vel.rain, FT(1)) D_min, D_max = FT(1e-5), FT(5e-3) - # v̄ᵢ far below v_l(D_min) ⇒ D* = D_min ; far above ⇒ D* = D_max. + # v_i far below v_l(D_min) ⇒ D* = D_min ; far above ⇒ D* = D_max. @test P3.crossover_diameter(FT(-1), v_l, D_min, D_max) == D_min @test P3.crossover_diameter(FT(1e6), v_l, D_min, D_max) == D_max # An interior target lands strictly inside the bracket. @@ -1146,13 +1044,13 @@ function test_p3_closed_form_rain_inner(FT) Dstar = P3.crossover_diameter(v_mid, v_l, D_min, D_max) @test D_min <= Dstar <= D_max @test isapprox(v_l(Dstar), v_mid; rtol = FT(1e-4)) - # Full closed path with v̄ᵢ outside the band stays finite and the + # Full closed path with v_i outside the band stays finite and the # absent-crossover side collapses (one piece zero). - rᵢ = FT(2e-4) + r_i = FT(2e-4) ai, bi, ci = CO.Chen2022_vel_coeffs(vel.rain, FT(1)) - for v̄ᵢ in (FT(-1), FT(1e6)) + for v_i in (FT(-1), FT(1e6)) N, M = P3.closed_rain_inner_NM( - FT(1e-3), v̄ᵢ, v_l, rᵢ, FT(1000), ai, bi, ci, + FT(1e-3), v_i, v_l, r_i, FT(1000), ai, bi, ci, D_min, D_max, FT(1e7), FT(5e-4), ) @test isfinite(N) && isfinite(M) @@ -1160,13 +1058,6 @@ function test_p3_closed_form_rain_inner(FT) end @testset "P3 closed-form dispatch fallback (non-Chen / non-SB2006)" begin - # A non-(RainParticlePDF_SB2006, Chen2022VelType) bundle must - # select the numerical `get_liquid_integrals` path (byte-unchanged - # behavior for other bundles). We assert method-table dispatch: - # the typed first/second args resolve to the closed-form method; - # any other first/second args resolve to the `::Any,::Any` - # numerical-fallback method. The trailing args are placeholders - # (only the first two participate in dispatch selection). params = CMP.ParametersP3(FT) vel = CMP.Chen2022VelType(FT) psd_r = CMP.SB2006(FT).pdf_r @@ -1175,10 +1066,8 @@ function test_p3_closed_form_rain_inner(FT) identity, (a, b) -> a, identity, identity, identity, FT(1), FT(1e-4), FT(1e3), state, ) - m_closed = - which(P3._rain_inner_integrals, typeof((psd_r, vel, rest...))) - m_fallback = - which(P3._rain_inner_integrals, typeof((1.0, 2.0, rest...))) + m_closed = which(P3._rain_inner_integrals, typeof((psd_r, vel, rest...))) + m_fallback = which(P3._rain_inner_integrals, typeof((1.0, 2.0, rest...))) # closed-form method: first two params are the typed bundle @test m_closed.sig.parameters[2] <: CMP.RainParticlePDF_SB2006 @test m_closed.sig.parameters[3] <: CMP.Chen2022VelType diff --git a/test/performance_tests.jl b/test/performance_tests.jl index 41f399756..31b4b3fe4 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -176,7 +176,7 @@ function benchmark_test(FT) P3.P3State, P3.P3State, (params_P3, L_ice, N_ice, F_rim, ρ_rim), - 200, # ns; loose enough for slower CI hardware (≈120 ns observed) + 220, ) bench_press(FT, P3.get_distribution_logλ, (state,), 30_000) # The weighted-velocity integrals build a nested terminal-velocity closure @@ -194,7 +194,7 @@ function benchmark_test(FT) @NamedTuple{dNdt::FT, dLdt::FT}, P3.het_ice_nucleation, (kaolinite, tps, q_liq, N_liq, RH_2, T_air_2, ρ_air), - 200, + 220, ) bench_press( @NamedTuple{dNdt::FT, dLdt::FT}, @@ -238,7 +238,7 @@ function benchmark_test(FT) micro_mock = (; q_tot = FT(0.00145), q_lcl = FT(0), q_icl = FT(0), q_rai = FT(0), q_sno = FT(0)) thermo_mock = (; ρ = FT(0.8), T = FT(263)) bench_press(FT, CMN.conv_q_vap_to_q_lcl, - (CMP.CloudLiquidFormation(CP.create_toml_dict(FT)), mp_mock, tps, micro_mock, thermo_mock), 140) + (CMP.CloudLiquidFormation(CP.create_toml_dict(FT)), mp_mock, tps, micro_mock, thermo_mock), 160) @info "0-Moment Scheme" bench_press(FT, CM0.remove_precipitation, (p0m, q_liq, q_ice), 12)