Skip to content

P3 quadrature: residual collision/B_rim accuracy floor and wet-growth onset root solve #761

Description

@haakon-e

Follow-up to #741, which set the P3 default quadrature to Gauss-Legendre order 6 and placed the known non-smooth features of the size-distribution integrands (mass-regime thresholds, the terminal-velocity regime break, the ice-liquid fall-speed crossing, the distribution-tail decay scale, and the wet-growth onset) on subinterval boundaries. Two accuracy limitations remain after that change, both order-independent, so raising the quadrature order does not reduce them.

1. Collision / B_rim accuracy floor (~1.4e-3)

Per #741: without the fall-speed-crossing breakpoint, the collision-source error stays at approximately 1.4e-3 for every quadrature order up to 32. With the breakpoint, the error-family that carries this component (collision_efficiency in test/p3_quadrature_error_study.jl) still stops shrinking well before the transport family does (GL(6) max 1.2e-2 down to GL(8) max 1.1e-3, versus the transport family's continued convergence).

Root cause: the fall-speed crossing diameter is located by crossover_diameter (src/P3_processes.jl:360), which calls RS.find_zero with RS.BrentsMethod and the FixedIterations{FT} tolerance object (src/P3_size_distribution.jl:250). FixedIterations's convergence predicate always returns false, so the solve always runs exactly maxiters iterations (8 for Float32, 10 for Float64) with no residual check and no use of the solver's converged flag. That budget was calibrated empirically for a different root problem, get_distribution_logλ's shape-parameter solve, which the surrounding comment (src/P3_size_distribution.jl:301) explicitly documents as "an EMPIRICAL, curvature-dependent result, NOT a guaranteed tolerance", close to linear over its bracket. crossover_diameter and the onset solve below reuse the same fixed budget on structurally different functions without re-deriving or re-checking that budget for their curvature.

The resulting crossing-diameter error propagates through two paths that both consume the same Dstar:

  • get_liquid_integrals_rain_closed (src/P3_processes.jl:418) uses Dstar as the split point of the closed-form incomplete-gamma difference in closed_rain_inner_NM (crossing(p) = flux(D_min, Dstar, p) - flux(Dstar, D_max, p), line 401). An inexact Dstar biases this split even though the two halves are each evaluated in exact closed form.
  • The same Dstar is placed as a subinterval boundary for the B_rim quadrature (src/P3_processes.jl:443), so an inexact Dstar leaves the true |v_i - v_l| derivative discontinuity slightly off the subinterval boundary, degrading Gauss-Legendre's geometric convergence on that subinterval regardless of order.
  • crossing_integral_bounds (src/P3_processes.jl:348), used by the general (non-closed-form) liquid path, has the same dependence on crossover_diameter.

Because the error source is the location precision of Dstar, not the quadrature discretization, none of these components improve by raising the Gauss-Legendre order.

2. Wet-growth onset root solve (~0.7% onset placement, ~7e-3 wet-growth indicator at order 12)

wet_growth_onset_diameter (src/P3_processes.jl:631) locates up to two onset crossings with a 16-point log-spaced scan (n_scan = 16, line 676) for sign changes, followed by the same FixedIterations/BrentsMethod refine (lines 678-686) used above. Per #741, the closed-form balance locates the onset within about 0.7% of the exact value across tested regimes, and the wet-growth indicator component of the bulk collision source holds at approximately 7e-3 at order 12, "set by the onset placement", while the freeze/shed mass split converges much further (4e-5) since it responds to the onset location only through the smooth min() partition.

Two mechanisms limit this component, independent of quadrature order:

  • The fixed-iteration Brent refine caps the onset diameter's placement precision the same way as crossover_diameter above.
  • The 16-point log-spaced scan can miss a wet-growth window narrower than one scan interval (two sign changes inside a single bracket are invisible to the pairwise sign-change test), or resolve a window whose true width is comparable to Δl = (lhi - llo) / n_scan.

Proposed directions

  • Give crossover_diameter and the onset refine_crossing a stopping tolerance derived for their own curvature (a residual- or bracket-width-based RS.AbstractTolerance, or an empirically re-derived fixed iteration count following the same approach used for get_distribution_logλ), rather than reusing the logλ-calibrated budget unchanged.
  • Quantify the sensitivity of closed_rain_inner_NM's N/M split to Dstar placement error directly, to determine whether tightening Dstar alone removes the ~1.4e-3 floor or whether the B_rim quadrature's off-boundary derivative discontinuity remains the dominant term.
  • For the onset window, increase n_scan or use an adaptive/bisecting scan when the log-diameter range is wide, and add a regression check for the fraction of test states whose wet-growth window is narrower than one scan interval.
  • Add unit-level accuracy tests for crossover_diameter and wet_growth_onset_diameter against an independent high-precision reference solve, rather than observing their effect only indirectly through the bulk error study.

Measurement

test/p3_quadrature_error_study.jl (run_quadrature_error_study()) is the harness that produced the numbers above and should be rerun to confirm any improvement; see its docstring for usage.

Reference: #741.

Metadata

Metadata

Assignees

No one assigned

    Labels

    P3Predicted particle properties

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions