feat(P3): Gauss-Legendre quadrature, integrand breakpoints -> GL(6) default

docs/src/API.md

@@ -221,6 +221,7 @@ These methods integrate over the particle size distribution.
 P3Scheme.D_m
 P3Scheme.ice_particle_terminal_velocity
 P3Scheme.ice_terminal_velocity_number_weighted
+P3Scheme.ice_terminal_velocity_number_weighted_from_prognostic
 P3Scheme.ice_terminal_velocity_mass_weighted
 ```
 
@@ -247,10 +248,19 @@ P3Scheme.bulk_liquid_ice_collision_sources
 Supporting methods:
 
 ```@docs
+P3Scheme.VolumetricCollisionRate
 P3Scheme.volumetric_collision_rate_integrand
 P3Scheme.compute_max_freeze_rate
 P3Scheme.compute_local_rime_density
 P3Scheme.get_liquid_integrals
+P3Scheme.crossing_integral_bounds
+P3Scheme.crossover_diameter
+P3Scheme.closed_rain_inner_NM
+P3Scheme.collision_cross_section_ice_liquid
+P3Scheme.collision_cross_section_ice_liquid_coeffs
+P3Scheme.collision_cross_section_ice_ice
+P3Scheme.ice_self_collection
+P3Scheme.wet_growth_onset_diameter
 P3Scheme.∫liquid_ice_collisions
 ```
 
@@ -262,8 +272,9 @@ P3Scheme.integrate
 P3Scheme.subintervals
 P3Scheme.ChebyshevGauss
 Quadrature.GaussLegendre
-Quadrature.build_quadrature
 P3Scheme.integral_bounds
+P3Scheme.velocity_integral_bounds
+P3Scheme.velocity_breakpoints
 ```
 
 # Aerosol model
@@ -406,6 +417,8 @@ Utilities.fac
 
 ```@docs
 Common
+Common.Chen2022VelocityCurve
+Common.particle_terminal_velocity
 Common.G_func_liquid
 Common.G_func_ice
 Common.logistic_function

docs/src/DevelopersGuide.md

@@ -206,3 +206,25 @@ Add a comment of what you are testing and use `@test` to create your test.
 The GPU tests are ran twice: for `Float64` and `Float32`.
 Similar as with performance tests, some trial and error is needed
   to find good tolerances for both options.
+
+### P3 quadrature error study
+
+The P3 scheme evaluates its size-distribution integrals with a Gauss-Legendre
+rule whose default order is recorded in `src/parameters/Microphysics2MParams.jl`.
+The study that supports the default lives in `test/p3_quadrature_error_study.jl`:
+it measures the relative error of every quadrature-consuming P3 quantity against
+a high-order reference, across a set of physically plausible column states, and
+attributes it to the transport components (sedimentation velocities, melt), the
+components proportional to the collision efficiency, and the full bulk tendency
+vector.
+
+Run it from the repository root with
+
+```
+julia --project=test -e 'include("test/p3_quadrature_error_study.jl");
+                         run_quadrature_error_study()'
+```
+
+Rerun the study when changing the quadrature rule, the integral bounds or
+breakpoints, or the P3 process integrands, and revisit the default order and the
+tolerances in `test/bulk_tendencies_quadrature_tests.jl` accordingly.

docs/src/plots/P3Melting.jl

@@ -9,6 +9,7 @@ FT = Float64
 # parameters
 params = CMP.ParametersP3(FT; slope_law = :constant)
 vel = CMP.Chen2022VelType(FT)
+quad = CMP.Microphysics2MParams(FT; with_ice = true).ice.quad  # the package default rule
 aps = CMP.AirProperties(FT)
 tps = TDI.PS(FT)
 
@@ -34,31 +35,31 @@ Fᵣ = FT(0.8)
 label1 = "ρₐ=$ρₐ kg/m³, Fᵣ=$Fᵣ, ρᵣ=$(Int(ρᵣ)) kg/m³"
 state = P3.P3State(params, Lᵢ, Nᵢ, Fᵣ, ρᵣ)
 logλ = P3.get_distribution_logλ(state)
-melt1 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ)
+melt1 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ; quad)
 dLdt1 = getfield.(melt1, :dLdt)
 dNdt1 = getfield.(melt1, :dNdt)
 
 Fᵣ = FT(0.2)
 label2 = "ρₐ=$ρₐ kg/m³, Fᵣ=$Fᵣ, ρᵣ=$(Int(ρᵣ)) kg/m³"
 state = P3.P3State(params, Lᵢ, Nᵢ, Fᵣ, ρᵣ)
 logλ = P3.get_distribution_logλ(state)
-melt2 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ)
+melt2 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ; quad)
 dLdt2 = getfield.(melt2, :dLdt)
 dNdt2 = getfield.(melt2, :dNdt)
 
 ρᵣ = FT(200)
 label3 = "ρₐ=$ρₐ kg/m³, Fᵣ=$Fᵣ, ρᵣ=$(Int(ρᵣ)) kg/m³"
 state = P3.P3State(params, Lᵢ, Nᵢ, Fᵣ, ρᵣ)
 logλ = P3.get_distribution_logλ(state)
-melt3 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ)
+melt3 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ; quad)
 dLdt3 = getfield.(melt3, :dLdt)
 dNdt3 = getfield.(melt3, :dNdt)
 
 ρₐ = FT(0.5)
 label4 = "ρₐ=$ρₐ kg/m³, Fᵣ=$Fᵣ, ρᵣ=$(Int(ρᵣ)) kg/m³"
 state = P3.P3State(params, Lᵢ, Nᵢ, Fᵣ, ρᵣ)
 logλ = P3.get_distribution_logλ(state)
-melt4 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ)
+melt4 = P3.ice_melt.(vel, aps, tps, params.T_freeze .+ ΔT_range, ρₐ, state, logλ; quad)
 dLdt4 = getfield.(melt4, :dLdt)
 dNdt4 = getfield.(melt4, :dNdt)
 

docs/src/plots/P3TerminalVelocityPlots.jl

@@ -27,6 +27,7 @@ function get_values(
     D_m = zeros(x_resolution, y_resolution)
     D_m_regimes = zeros(x_resolution, y_resolution)
     ϕᵢ = zeros(x_resolution, y_resolution)
+    quad = CMP.Microphysics2MParams(FT; with_ice = true).ice.quad  # the package default rule
 
     for i in 1:x_resolution
         for j in 1:y_resolution
@@ -35,8 +36,8 @@ function get_values(
             state = P3.P3State(params, L, N, F_rim, ρ_rim)
             state_noar = P3.P3State(params_noar, L, N, F_rim, ρ_rim)
             logλ = P3.get_distribution_logλ(state)
-            V_m[i, j] = P3.ice_terminal_velocity_mass_weighted(Chen2022, ρ_a, state_noar, logλ)
-            V_m_ϕ[i, j] = P3.ice_terminal_velocity_mass_weighted(Chen2022, ρ_a, state, logλ)
+            V_m[i, j] = P3.ice_terminal_velocity_mass_weighted(Chen2022, ρ_a, state_noar, logλ; quad)
+            V_m_ϕ[i, j] = P3.ice_terminal_velocity_mass_weighted(Chen2022, ρ_a, state, logλ; quad)
             D_m[i, j] = P3.D_m(state, logλ)
             D_m_regimes[i, j] = D_m[i, j]
             ϕᵢ[i, j] = P3.ϕᵢ(state, D_m[i, j])

src/P3_integral_properties.jl

@@ -41,7 +41,26 @@ Compute the integration bounds for the P3 size distribution,
 
     # Only integrate up to the maximum diameter, `D_max`, including intermediate thresholds
     # If `F_rim` is very close to 1, `D_cr` may be greater than `D_max`, in which case it is disregarded.
-    return segment_boundaries(state, D_min, D_max)
+    bnds = segment_boundaries(state, D_min, D_max)
+    # `D_max` sits `log(1/p)` decay lengths into the size-distribution tail; a
+    # breakpoint at the decay scale keeps each subinterval resolvable at low order
+    D_e = clamp(3 / λ, D_min, D_max)
+    return Tuple(SA.sort(SA.SVector(bnds..., D_e)))
+end
+
+"""
+    velocity_integral_bounds(state::P3State, logλ, v_term; p, moment_order = 0)
+
+Compute the integration bounds for a velocity-weighted P3 integral: the
+mass-regime [`integral_bounds`](@ref) with the [`velocity_breakpoints`](@ref)
+of the terminal-velocity closure `v_term` clamped into `[D_min, D_max]` and
+re-sorted, so each breakpoint coincides with a subinterval boundary. Returns a
+fixed-length tuple.
+"""
+function velocity_integral_bounds(state::P3State{FT}, logλ, v_term::V; p, moment_order = 0) where {FT, V}
+    bnds = integral_bounds(state, logλ; p, moment_order)
+    breaks = map(D -> clamp(FT(D), first(bnds), last(bnds)), velocity_breakpoints(v_term))
+    return Tuple(SA.sort(SA.SVector(bnds..., breaks...)))
 end
 
 """