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

[haakon-e]

Set the default P3 quadrature rule to Gauss-Legendre of order 6, and place every known feature of the P3 size-distribution integrands on a quadrature subinterval boundary so that a low-order rule suffices.

Stacked on #757 (mechanical refactors and small fixes).

---

Each P3 size-distribution integrand has a small set of known non-smooth features: the mass-regime thresholds, the terminal-velocity regime break, the ice-liquid fall-speed crossing, the exponential decay of the distribution tail, and the wet-growth onset. Placing each feature on a subinterval boundary leaves the integrand smooth within every subinterval, where Gauss-Legendre converges geometrically. Order 6 then meets (a handwavily reasonable) target accuracy (see below).

API

• The quadrature rule lives on P3IceParams. build_quadrature and the stored quadrature_order field are removed; the parameter constructors take a quadrature_order (or quad) keyword.
• integrate and the P3 process functions require quad.
• VolumetricCollisionRate stores the ice and liquid terminal-velocity closures as fields, so the collision integrals read the velocities and their structure from the integrand itself: the crossing bounds, the closed-form rain integrals (coefficients read from the closure), and the wet-growth onset balance.
• velocity_breakpoints(v_term) returns the diameters where a terminal-velocity closure changes functional form. velocity_integral_bounds consumes the closure, so the integral bounds no longer assume the Chen 2022 cutoff field.

Breakpoints

• Terminal-velocity regime break (Chen small/large ice): a subinterval boundary (velocity_integral_bounds) for the velocity-weighted integrals, ice_melt (through its ventilation factor), the outer collision integral, and ice self-collection.
• Ice-liquid fall-speed crossing v_l(D) = v_i(Dᵢ): a subinterval boundary of the inner collision integrals (crossing_integral_bounds). The rain path computes the crossover diameter once per outer node and shares it across the closed-form N and M components and the B_rim quadrature. With the crossing on a boundary, the collision-efficiency family (collision sources and ice self-collection) converges with order rather than saturating at a fixed error: worst-case relative error falls from 1.2e-2 at order 6 to 5.0e-6 at order 32.
• Distribution tail at three decay lengths (3/λ): the upper bound is the (1 - p) quantile, so without this breakpoint the last subinterval spans most of the log(1/p) ≈ 11.5 decay lengths and sets the low-order error in large-mean-size states. With the breakpoint, ice self-collection in the hail-core states reaches 7.1e-5 relative error at order 12 (1.1e-2 at order 6).
• Wet-growth onset (where collected liquid mass balances the Musil freeze limit): a step in the outer integrand, through the freeze/shed min() partition and the wet-growth indicator. wet_growth_onset_diameter locates up to two balance points (the window closes again at large sizes) from a closed-form balance refined by bracketed root-finding, and places them on outer subinterval boundaries. At order 12 the shedding component reaches ~1e-7 relative error; the wet-growth-indicator component sits near 2e-4 and is nearly order-independent (2.5e-4 at order 6, 2.3e-4 at order 12), set by the onset placement rather than the node count.
• Ice self-collection integrates over the upper triangle D₁ ≤ D₂, which places the D₁ = D₂ derivative discontinuity of the relative fall speed on a boundary and halves the integrand evaluations. A test cross-checks the triangle against the full-square double integral with the ½ pair-counting factor.

Accuracy and the default order

I wrote a script that let's me test the accuracy for different orders (since I suspect I'll revisit this in the future).
Over each P3 integral, the collision source vector, and the full 2M+P3 bulk tendency vector, across 17 "regime" states including hail cores. Baseline ("truth") is a GL(128) reference (converged). The main column runs the same script on main against main's own GL(128) reference (both without breakpoints). The bulk-tendency error separates into transport components (sedimentation velocities, melt) and components proportional to the collision efficiency (collision sources, ice self-collection, rime):

	main (GL(16) default)	GL(5)	GL(6) (default)	GL(7)	GL(8)
transport p95 / max	—	1.7e-3 / 2.0e-2	7.1e-4 / 4.5e-3	3.2e-4 / 1.2e-3	4.4e-5 / 3.9e-4
E-family p95 / max	—	3.9e-2 / 6.9e-2	1.1e-2 / 1.2e-2	2.0e-3 / 7.8e-3	6.1e-4 / 1.1e-3
bulk p95 / max	1.4e-3 / 6.2e-2	2.8e-2 / 6.9e-2	8.9e-3 / 1.4e-2	1.5e-3 / 7.8e-3	2.9e-4 / 9.0e-4

Every bulk error above 0.45% is in a collision-efficiency component, which sets E = 1. The package's own single-moment scheme assigns collision efficiencies from 0.1 to 1.0 for the analogous liquid-frozen and rain-frozen collisions (cloud_liquid_snow_collision_efficiency = 0.1, cloud_ice_rain_collision_efficiency = 1.0), so the E = 1 choice is an upper bound uncertain to a factor of order ten. That uncertainty exceeds the ~1% quadrature error on those rates, so resolving them below ~1% adds no physical accuracy while E is fixed. The transport components, which do not include the E = 1 assumption, converge to at most 0.45% (worst-case 4.5e-3 at order 6). Order 5 reaches 2.0e-2 in the transport family and is rejected. Order 7 matches main's bulk p95 and is far below main on the worst-case bulk error; order 8 is below main on both. Both remain one-keyword changes (quadrature_order). Parameterizing the collision efficiency is the top follow-up, and this table should be revisited when it lands.

Runtime

Then testing the runtime with these changes in Float64 on Julia 1.12.

cell	main (GL(16))	this PR (GL(6) + breakpoints)	ratio
warm (ice-free)	0.48 µs	0.81 µs	1.7
mixed-phase	5011 µs	2063 µs	0.41
hail core	4824 µs	2143 µs	0.44

The warm cell performs no size-distribution integrals; its sub-microsecond time differs between the two by about 0.3 µs, reproducibly across runs and negligible against the ice-path cost. The he/p3-refactors base (the PR base) gives the same GL(16) baseline as main (5011 / 4824 µs here versus 5017 / 4869 µs for main in a separate run, for the two ice cells), so the speedup is attributable to the order reduction, the closed-form rain integrals, and the triangle self-collection, net of the onset cost. The onset location costs about a quarter of the collision call (a log-spaced scan, RootSolvers Brent refinements of the closed-form balance, and two additional outer subintervals). Removing the step also removes a state-dependent discontinuity from the tendencies that the implicit solver's Jacobian must differentiate.

The absolute GL(16) baseline depends on the Julia version: on Julia 1.12 the 16-node collision double integral runs at ~5000 µs, roughly 1.8× the figure a lower Julia version reports, while the GL(6) time is stable across versions. The ratio above is therefore the Julia 1.12 result (on lower versions the ratio is closer to 0.75). The measurement ran on the clima cluster.

---

Tests:

• test/p3_quadrature_error_study.jl (run_quadrature_error_study(); usage in the developer guide) regenerates the accuracy table when the quadrature machinery or the process integrands change.
• test/bulk_tendencies_quadrature_tests.jl locks orders 6 through 12 as a regression test against a GL(200) reference, with graupel and hail-core states in the state set.
• The P3 and 2M integral tests use Gauss-Legendre in place of Chebyshev-Gauss, matching the default: scheme-path calls use GL(6), and the standalone reference integrals use the smallest order meeting their tolerance. On these integrands Gauss-Legendre is 1-2 orders of magnitude more accurate than Chebyshev-Gauss at matched order and cost; the endpoint behavior that originally motivated Chebyshev-Gauss was the unresolved crossing and onset structure, which the breakpoints remove.

[haakon-e]

This change is part of the following stack:

• feat(P3): Gauss-Legendre quadrature, integrand breakpoints -> GL(6) default #741 ◀
  • feat(Rosenbrock): unified RosenbrockAverage substepping framework + 2M+P3 exact-AD mode #726
    • feat(Rosenbrock): 1M donor / coupled / exact modes on the unified framework #731
      • feat(Rosenbrock): post-solve per-process verbose diagnostics #732
        • fix(Rosenbrock): saturation adjustment on max(S_ice, S_liq) #737
          • perf(BMT): recover the primal tendency from the Rosenbrock Jacobian pass #742
            • feat(BMT): analytic ManualJacobian for the 2M+P3 Rosenbrock substep #743
              • feat(BMT): coupled Tier-3 donor linearization for the 2M+P3 ManualJacobian #744
  • refactor(Common,P3): consolidate the Chen 2022 velocity curves #762

_{Change managed by git-spice.}

[codecov]

Codecov Report

❌ Patch coverage is 96.66667% with 4 lines in your changes missing coverage. Please review.
✅ Project coverage is 92.81%. Comparing base (cd6581c) to head (355b541).
⚠️ Report is 3 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #741      +/-   ##
==========================================
- Coverage   92.92%   92.81%   -0.11%     
==========================================
  Files          56       56              
  Lines        2896     2952      +56     
==========================================
+ Hits         2691     2740      +49     
- Misses        205      212       +7     

Components	Coverage Δ
src	93.59% <96.66%> (-0.13%)	⬇️
ext	69.47% <ø> (ø)

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[haakon-e]

just rewriting without branches here an a few other places.

[trontrytel]

Claude is asking if ForwardDiff will propagate the non-finite rate through FT(rate) before ifelse discards it. Will Dual arithmetic with Inf work for the derivative?