feat(P3): Gauss-Legendre quadrature parameter, velocity-aware bounds, and P3 property refactors

Once every integrand kink is a quadrature subinterval boundary, the P3
size-distribution integrals are smooth, so a low-order Gauss-Legendre rule
reaches the accuracy that previously needed a high-order Chebyshev-Gauss rule.

Quadrature:
- Default to Gauss-Legendre at every order (arbitrary n >= 1); at matched order
  it is ~20x more accurate than Chebyshev-Gauss on the smooth P3 integrands,
  verified against an adaptive QuadGK reference.
- Pass the quadrature rule directly on the parameters (P3IceParams.quad); remove
  build_quadrature and the quadrature_order field/kwarg, so the rule is a single
  explicit object rather than an order-to-rule indirection.
- Lower the default order to 12 (~1e-5 worst-case relative error on the
  domain-reduced integrals, the accuracy previously reached at order 16).
- Make quad a required argument of integrate and the P3 process functions,
  removing the hidden ChebyshevGauss(100) default; the rule is threaded from the
  parameters.
- Resolve the Chen small/large-ice velocity D_cutoff kink as a subinterval
  boundary (velocity_integral_bounds), so velocity-weighted integrals converge
  geometrically instead of plateauing at the kink.

Ice self-collection:
- Integrate the double integral over the upper triangle D_1 <= D_2, which folds
  the D_1 = D_2 diagonal and every regime breakpoint onto quadrature bounds and
  halves the integrand evaluations.

Particle properties:
- regime_value: one branchless selector for the five P3 mass/area regimes,
  replacing the duplicated regime conditionals in ice_area, ice_mass_coeffs, and
  ϕ_material_density.
- ϕᵢ uses the per-regime material density (oblate spheroid model), with aᵢ^{3/2}
  written as aᵢ * sqrt(aᵢ), a hardware sqrt and a multiply in place of a
  runtime-exponent power call.
- Select the ice terminal-velocity aspect-ratio correction by type
  (UseAspectRatio / NoAspectRatio) so the branch resolves at compile time.

Robustness:
- Floor the air density to non-negative in the Chen velocity coefficients; a
  slightly-negative sub-domain density otherwise throws in their fractional
  powers.

Author: haakon-e

docs/src/API.md

@@ -189,6 +189,7 @@ P3Scheme.ice_mass
 P3Scheme.ice_density
 P3Scheme.∂ice_mass_∂D
 P3Scheme.ice_area
+P3Scheme.ϕ_material_density
 P3Scheme.ϕᵢ
 ```
 
@@ -261,7 +262,6 @@ P3Scheme.integrate
 P3Scheme.subintervals
 P3Scheme.ChebyshevGauss
 Quadrature.GaussLegendre
-Quadrature.build_quadrature
 P3Scheme.integral_bounds
 ```
 

docs/src/P3Scheme.md

@@ -354,8 +354,42 @@ An alternative parameterization for $μ$ is a constant value:
 We use the Chen et al. [Chen2022](@cite) velocity parameterization,
   see [here](https://clima.github.io/CloudMicrophysics.jl/dev/TerminalVelocity.html#Chen-et.-al.-2022) for details.
 
+### Aspect ratio
+
+We model the ice particle as an oblate spheroid of maximum dimension $D$, projected
+area $a_i = $ `ice_area(state, D)`, and mass $m_i = $ `ice_mass(state, D)`. Writing
+the equatorial diameter as $a = 2\sqrt{a_i/π}$ and the polar diameter as $c = ϕ\,a$,
+and equating the spheroid mass $ρ\,(π/6)\,a^2 c$ to $m_i$ gives the oblate aspect
+ratio ($κ = 1/3$)
+
+```math
+ϕ = \frac{3 \sqrt{π}\, m_i}{4\, ρ\, a_i^{3/2}},
+```
+
+where $ρ$ is the particle's material density [`ϕ_material_density`](@ref) — the
+density of the solid the particle is made of ($ρ_i$, or $ρ_g$ for graupel) — not the
+size-dependent effective density $m_i / V_{sphere}(D)$ [`ice_density`](@ref). Using
+the effective density would cancel the particle mass and pin $ϕ ≡ 1$, disabling the
+$\sqrt[3]{ϕ}$ fall-speed correction in [`ice_particle_terminal_velocity`](@ref).
+
+Within each mass regime the material density is constant, so $ϕ$ tracks $m_i / a_i^{3/2}$.
+In the spherical regimes ($D < D_{th}$, and graupel $D_{gr} ≤ D < D_{cr}$) the mass and
+area are exactly spherical, so $ϕ = 1$. In the dense-nonspherical regime $ϕ < 1$ (oblate)
+and decreases with $D$, slowing large unrimed ice.
+
+!!! note "Residual $ϕ > 1$ band just above $D_{th}$"
+    $ϕ$ can slightly exceed 1 (peak $≈ 1.2$) in a narrow size band immediately above
+    $D_{th}$. This is not a sign error: it is caused by the area discontinuity at
+    $D_{th}$, where the projected area drops from the spherical law $(π/4) D^2$ to the
+    nonspherical law $γ D^σ$ while the mass stays continuous. The smaller area inflates
+    $m_i / a_i^{3/2}$ just past the threshold. $\sqrt[3]{ϕ}$ remains finite and bounded
+    there ($\sqrt[3]{1.2} ≈ 1.06$), so no clamp is applied. Reconciling the mass and area
+    laws at $D_{th}$ (so the area is continuous too) is a separate area power-law item,
+    left untouched here.
+
 The figure below shows the implied aspect ratio for different particle size regimes.
-Note that $ϕ = 1$ corresponds to spherical particles (small spherical ice ($D < D_{th}$) and graupel ($D_{gr} < D < D_{cr}$)).
+Note that $ϕ = 1$ corresponds to spherical particles (small spherical ice ($D < D_{th}$) and graupel ($D_{gr} < D < D_{cr}$)),
+and that $ϕ$ exceeds 1 in a narrow band just above $D_{th}$ (the area discontinuity noted above; a follow-on item).
 ```@example
 include("plots/P3AspectRatioPlot.jl")
 

docs/src/plots/P3AspectRatioPlot.jl

@@ -10,7 +10,7 @@ axis_theme = Theme(
         xminorticksvisible = true,
         xminorticks = IntervalsBetween(5),
         xticks = [0.01, 0.1, 1, 10],
-        limits = ((0.01, 10.0), (0, 1.05)),
+        limits = ((0.01, 10.0), (0, 1.3)),
         xgridvisible = false,
         ygridvisible = false,
         xlabel = "D (mm)",

src/Common.jl

@@ -289,6 +289,7 @@ See [Chen2022](@cite) for more details.
 """
 @inline function Chen2022_vel_coeffs(coeffs::CMP.Chen2022VelTypeRain, ρₐ)
     (; ρ0, a, a3_pow, b, b_ρ, c) = coeffs
+    ρₐ = max(ρₐ, zero(ρₐ))
     # Table B1
     q = exp(ρ0 * ρₐ)
     ai = (a[1] * q, a[2] * q, a[3] * q * ρₐ^a3_pow)
@@ -303,6 +304,7 @@ end
 @inline function Chen2022_vel_coeffs(coeffs::CMP.Chen2022VelTypeSmallIce, ρₐ, ρᵢ)
     FT = eltype(coeffs)
     (; A, B, C, E, F, G) = coeffs
+    ρₐ = max(ρₐ, zero(ρₐ))
     # Table B3 - cache sqrt for reuse
     log_ρᵢ = log(ρᵢ)
     sqrt_ρᵢ = sqrt(ρᵢ)
@@ -325,6 +327,7 @@ end
 @inline function Chen2022_vel_coeffs(coeffs::CMP.Chen2022VelTypeLargeIce, ρₐ, ρᵢ)
     FT = eltype(coeffs)
     (; A, B, C, E, F, G, H) = coeffs
+    ρₐ = max(ρₐ, zero(ρₐ))
     # Table B5 - cache sqrt for reuse
     log_ρᵢ = log(ρᵢ)
     sqrt_ρᵢ = sqrt(ρᵢ)

src/P3_integral_properties.jl

@@ -42,6 +42,20 @@ function integral_bounds(state::P3State{FT}, logλ; p, moment_order = 0) where {
     return segment_boundaries(state, D_min, D_max)
 end
 
+"""
+    velocity_integral_bounds(state::P3State, logλ, D_cutoff; p, moment_order = 0)
+
+Compute the integration bounds for a velocity-weighted P3 integral: the
+mass-regime [`integral_bounds`](@ref) with the Chen 2022 small/large-ice velocity
+breakpoint `D_cutoff` clamped into `[D_min, D_max]` and re-sorted, so the `v(D)`
+kink lands on a subinterval boundary. Returns a fixed 6-tuple.
+"""
+function velocity_integral_bounds(state::P3State{FT}, logλ, D_cutoff; p, moment_order = 0) where {FT}
+    bnds = integral_bounds(state, logλ; p, moment_order)
+    D_c = clamp(FT(D_cutoff), first(bnds), last(bnds))
+    return Tuple(SA.sort(SA.SVector(bnds..., D_c)))
+end
+
 """
     D_m(state::P3State, logλ)
 

src/P3_particle_properties.jl

@@ -304,6 +304,33 @@ function weighted_average(f_a, a, b)
     return f_a * a + (1 - f_a) * b
 end
 
+"""
+    regime_value(state::P3State, D; small, unrimed, dense_rimed, graupel, partially_rimed)
+
+Select the value for the P3 mass/area regime that the maximum dimension `D` falls in:
+ - `small`:             small spherical ice     (`D < D_th`),
+ - `unrimed`:           large unrimed ice       (`F_rim = 0` ∧ `D_th ≤ D`),
+ - `dense_rimed`:       dense rimed ice         (`D_th ≤ D < D_gr`),
+ - `graupel`:           graupel (rimed)         (`D_gr ≤ D < D_cr`),
+ - `partially_rimed`:   partially rimed ice     (`D_cr ≤ D`).
+
+All five values are evaluated and must share a type.
+"""
+@inline function regime_value(
+    state::P3State, D;
+    small, unrimed, dense_rimed, graupel, partially_rimed,
+)
+    (; F_rim, D_th, D_gr, D_cr) = state
+    #! format: off
+    return ifelse(D < D_th,      small,             # small spherical ice
+           ifelse(iszero(F_rim), unrimed,           # large nonspherical unrimed ice
+           ifelse(D < D_gr,      dense_rimed,       # dense nonspherical rimed ice
+           ifelse(D < D_cr,      graupel,           # graupel (rimed)
+                                 partially_rimed,   # partially rimed ice
+    ))))
+    #! format: on
+end
+
 """
     ice_mass_coeffs(state::P3State, D)
 
@@ -318,21 +345,16 @@ Return the coefficients for the ice mass power law at diameter `D`.
 """
 function ice_mass_coeffs(state::P3State, D)
     FT = promote_type(eltype(state), eltype(D))
-    (; params, F_rim, ρ_g, D_th, D_gr, D_cr) = state
+    (; params, F_rim, ρ_g) = state
     (; ρ_i) = params
     (; α_va, β_va) = params.mass
-    a, b = if D < D_th       # small spherical ice
-        (ρ_i * π / 6, 3)
-    elseif iszero(F_rim)     # large nonspherical unrimed ice
-        (α_va, β_va)
-    elseif D_th ≤ D < D_gr   # dense nonspherical rimed ice
-        (α_va, β_va)
-    elseif D_gr ≤ D < D_cr   # graupel (rimed)
-        (ρ_g * π / 6, 3)
-    else # D_cr ≤ D          # partially rimed ice
-        (α_va / (1 - F_rim), β_va)
-    end
-    return FT(a), FT(b)
+    m_sphere(ρ) = (FT(ρ * π / 6), FT(3))   # spherical mass ∝ ρ (π/6) D³
+    m_power(a) = (FT(a), FT(β_va))         # power-law mass ∝ a D^β_va
+    return regime_value(state, D;
+        small = m_sphere(ρ_i), unrimed = m_power(α_va),
+        dense_rimed = m_power(α_va), graupel = m_sphere(ρ_g),
+        partially_rimed = m_power(α_va / (1 - F_rim)),
+    )
 end
 
 """
@@ -361,7 +383,6 @@ Return the density of a particle with diameter `D`
 # Notes:
  The density of nonspherical particles is assumed to be the particle mass divided
  by the volume of a sphere with the same D [MorrisonMilbrandt2015](@cite).
- Needed for aspect ratio calculation, so we assume zero liquid fraction.
 """
 ice_density(state::P3State, D) = ice_mass(state, D) / CO.volume_sphere_D(D)
 
@@ -395,46 +416,63 @@ Return the cross-sectional area of a particle based on where it falls in the
  - `D`: maximum particle dimension [m]
 """
 function ice_area(state::P3State, D)
-    (; params, F_rim, D_th, D_gr, D_cr) = state
+    (; params, F_rim) = state
     (; γ, σ) = params.area
-    spherical_area(D) = D^2 * π / 4
-    nonspherical_area(D) = γ * D^σ
-    return if D < D_th       # small spherical ice
-        spherical_area(D)
-    elseif iszero(F_rim)     # large nonspherical unrimed ice
-        nonspherical_area(D)
-    elseif D_th ≤ D < D_gr   # dense nonspherical rimed ice
-        nonspherical_area(D)
-    elseif D_gr ≤ D < D_cr   # graupel (rimed)
-        spherical_area(D)
-    else # D_cr ≤ D          # partially rimed ice
-        weighted_average(F_rim, spherical_area(D), nonspherical_area(D))
-    end
+    spherical = D^2 * π / 4
+    nonspherical = γ * D^σ
+    return regime_value(state, D;
+        small = spherical, unrimed = nonspherical, dense_rimed = nonspherical,
+        graupel = spherical,
+        partially_rimed = weighted_average(F_rim, spherical, nonspherical),
+    )
+end
+
+"""
+    ϕ_material_density(state::P3State, D)
+
+Return the material density of the ice particle's solid at diameter `D`, used in
+the aspect-ratio closure [`ϕᵢ`](@ref). This is `ρ_i` in every mass regime except
+graupel (`D_gr ≤ D < D_cr`), where it is the graupel density `ρ_g`.
+
+This is the density of the actual solid material, not the size-dependent effective
+density `mᵢ / Vsphere(D)` returned by [`ice_density`](@ref); see the
+[P3 documentation](@ref "Aspect ratio") for the distinction.
+"""
+function ϕ_material_density(state::P3State, D)
+    (; params, ρ_g) = state
+    (; ρ_i) = params
+    return regime_value(state, D;
+        small = ρ_i, unrimed = ρ_i, dense_rimed = ρ_i, graupel = ρ_g, partially_rimed = ρ_i,
+    )
 end
 
 """
     ϕᵢ(state::P3State, D)
 
-Returns the aspect ratio (ϕ) for an ice particle with diameter `D`
+Return the oblate aspect ratio `ϕ = 3√π mᵢ / (4 ρ aᵢ^{3/2})` (`κ = 1/3`) for an
+ice particle of maximum dimension `D`, with mass `mᵢ = ice_mass(state, D)`,
+projected area `aᵢ = ice_area(state, D)`, and material density
+`ρ = ϕ_material_density(state, D)`. Assumes zero liquid fraction.
 
 # Arguments
  - `state`: The [`P3State`](@ref)
  - `D`: maximum dimension of ice particle [m]
 
-# Notes
- The density of nonspherical particles is assumed to be equal to the particle mass
- divided by the volume of a spherical particle with the same D_max [MorrisonMilbrandt2015](@cite).
- Assuming zero liquid fraction and oblate shape.
+See also [`ϕ_material_density`](@ref) and the
+[aspect-ratio section of the P3 documentation](@ref "Aspect ratio") for the
+spheroid derivation and the residual `ϕ > 1` band above `D_th`.
 """
 function ϕᵢ(state::P3State, D)
     FT = eltype(D)
     mᵢ = ice_mass(state, D)
     aᵢ = ice_area(state, D)
-    ρᵢ = mᵢ / CO.volume_sphere_D(D)
+    ρ = ϕ_material_density(state, D)
 
-    # TODO - prolate or oblate?
-    ϕ_ob = min(1, 3 * sqrt(FT(π)) * mᵢ / (4 * ρᵢ * aᵢ^FT(1.5))) # κ =  1/3
-    #ϕ_pr = max(1, 16 * ρᵢ^2 * aᵢ^3 / (9 * FT(π) * mᵢ^2))       # κ = -1/6
+    # Oblate aspect ratio (κ = 1/3); no clamp (see P3 documentation).
+    # `aᵢ^{3/2}` written as `aᵢ * sqrt(aᵢ)`: a hardware sqrt and a multiply
+    # instead of a runtime-exponent power call, evaluated per quadrature node.
+    ϕ_ob = 3 * sqrt(FT(π)) * mᵢ / (4 * ρ * aᵢ * sqrt(aᵢ))
+    #ϕ_pr = 16 * ρ^2 * aᵢ^3 / (9 * FT(π) * mᵢ^2)  # prolate, κ = -1/6
 
     return ifelse(D == 0, zero(ϕ_ob), ϕ_ob)
 end

src/P3_processes.jl

@@ -57,14 +57,14 @@ end
  - `logλ`: the log of the slope parameter [log(1/m)]
 
 # Keyword arguments
- - `quad`: quadrature rule, default is `ChebyshevGauss(100)`
+ - `quad`: quadrature rule (a `Quadrature.QuadratureRule`)
 
 Returns the melting rate of ice (QIMLT in Morrison and Mildbrandt (2015)).
 """
 function ice_melt(
     velocity_params, aps::CMP.AirProperties, tps::TDI.PS,
     Tₐ, ρₐ, state::P3State, logλ;
-    quad = ChebyshevGauss(100),
+    quad,
 )
     # Note: process not dependent on `F_liq`
     # (we want ice core shape params)
@@ -292,7 +292,7 @@ Returns a function `liquid_integrals(Dᵢ)` that computes the liquid particle in
 - `liq_bounds`: integration bounds for liquid particles
 
 # Keyword arguments
-- `quad`: quadrature rule, default is `ChebyshevGauss(100)`
+- `quad`: quadrature rule (a `Quadrature.QuadratureRule`)
 
 # Notes
 The function `liquid_integrals(Dᵢ)` returns a tuple `(∂ₜN_col, ∂ₜM_col, ∂ₜB_col)`
@@ -301,7 +301,7 @@ The function `liquid_integrals(Dᵢ)` returns a tuple `(∂ₜN_col, ∂ₜM_col
 - `∂ₜM_col`: mass collision rate [kg/s]
 - `∂ₜB_col`: rime volume collision rate [m³/s]
 """
-function get_liquid_integrals(n, ∂ₜV, m_liq, ρ′_rim, liq_bounds; quad = ChebyshevGauss(100))
+function get_liquid_integrals(n, ∂ₜV, m_liq, ρ′_rim, liq_bounds; quad)
     function liquid_integrals(Dᵢ)
         integrand =
             D -> SA.SVector(
@@ -440,12 +440,12 @@ Computes the bulk collision rate integrands between ice and liquid particles.
 - `ice_bounds`: integration bounds for ice particles, from [`integral_bounds`](@ref)
 
 # Keyword arguments
-- `quad`: quadrature rule, default is `ChebyshevGauss(100)`
+- `quad`: quadrature rule (a `Quadrature.QuadratureRule`)
 
 # Returns
 A tuple of 8 integrands, see [`∫liquid_ice_collisions`](@ref) for details.
 """
-function ∫liquid_ice_collisions(n_i, ∂ₜM_max, cloud_integrals, rain_integrals, ice_bounds; quad = ChebyshevGauss(100))
+function ∫liquid_ice_collisions(n_i, ∂ₜM_max, cloud_integrals, rain_integrals, ice_bounds; quad)
     function liquid_ice_collisions_integrands(Dᵢ)
         # Inner integrals over liquid particle diameters
         ∂ₜN_c_col, ∂ₜM_c_col, ∂ₜB_c_col = cloud_integrals(Dᵢ)
@@ -593,7 +593,7 @@ A `NamedTuple` of `(; ∂ₜq_c, ∂ₜq_r, ∂ₜN_c, ∂ₜN_r, ∂ₜL_rim, 
 function bulk_liquid_ice_collision_sources(
     state, logλ,
     psd_c, psd_r, L_c, N_c, L_r, N_r,
-    aps, tps, vel, ρₐ, T; quad = ChebyshevGauss(100),
+    aps, tps, vel, ρₐ, T; quad,
 )
     FT = promote_type(eltype(state), UT.promote_typeof(L_c, N_c, L_r, N_r, ρₐ, T))
     (; τ_wet, ρ_i) = state.params
@@ -687,28 +687,34 @@ while leaving mass, rime mass, and rime volume unchanged.
 - `ρₐ`: air density [kg/m³]
 
 # Keyword arguments
-- `quad`: quadrature rule, default is `ChebyshevGauss(100)`
+- `quad`: quadrature rule (a `Quadrature.QuadratureRule`)
 
 # Returns
 A `NamedTuple` of `(; dNdt)`, where:
 1. `dNdt`: ice number concentration tendency due to self-collection [1/m³/s] (always positive or zero, represents a loss rate)
 """
-function ice_self_collection(state, logλ, vel, ρₐ; quad = ChebyshevGauss(100))
+function ice_self_collection(state, logλ, vel, ρₐ; quad)
     n_i = DT.size_distribution(state, logλ)
     ∂ₜV = volumetric_ice_ice_collision_rate_integrand(vel, ρₐ, state)
 
     p = eps(one(ρₐ))
-    ice_bounds = integral_bounds(state, logλ; p)
+    ice_bounds = velocity_integral_bounds(state, logλ, vel.small_ice.cutoff; p)
+    D_min, D_max = ice_bounds[1], ice_bounds[end]
 
     function inner_integral(D_1)
+        # Inner integral over D_2 ∈ [D_1, D_max] (the upper triangle). Its
+        # subinterval boundaries are the P3 regime breakpoints restricted to that
+        # window: clamping each breakpoint up to D_1 drops those below it to
+        # zero-width (no-op) subintervals.
         integrand = D_2 -> ∂ₜV(D_1, D_2) * n_i(D_2)
-        rate_at_D1 = integrate(integrand, ice_bounds, quad)
-        return rate_at_D1 * n_i(D_1)
+        D_lo = clamp(D_1, D_min, D_max)
+        inner_bounds = map(D -> max(D, D_lo), ice_bounds)
+        return integrate(integrand, inner_bounds, quad) * n_i(D_1)
     end
 
-    total_rate = integrate(inner_integral, ice_bounds, quad)
-
-    # The 0.5 factor accounts for double-counting in self-collection
-    dNdt = (1 // 2) * total_rate
+    # The collision integrand is symmetric under `D_1 ↔ D_2`, so the triangle
+    # integral times two equals the full-square integral; that factor of two
+    # cancels the `1/2` self-collection double-counting factor.
+    dNdt = integrate(inner_integral, ice_bounds, quad)
     return (; dNdt)
 end