CliMA / CliMA/CloudMicrophysics.jl
P3 quadrature: residual collision/B_rim accuracy floor and wet-growth onset root solve
- Dominant language
- Julia
- Stars
- 55
- Forks
- 11
- Avg merge
- 1d 5h
- Merged PRs (30d)
- 8
Description
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.
Contributor guide
No contributing guide indexed for this repository
Research direction
Start with crossover_diameter and wet_growth_onset_diameter in src/P3_processes.jl, then read FixedIterations in src/P3_size_distribution.jl and its use in get_distribution_logλ. Run test/p3_quadrature_error_study.jl and inspect the referenced error families before choosing how to validate the root solves. Done means the collision/B_rim and wet-growth accuracy floors improve, with unit-level checks against an independent high-precision reference.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- julia
- Domain
- performance, testing
- Issue type
- Bug
- Difficulty
- 4/5
- Estimated time
- 3-5 days
- Activity status
- Quiet
- Clarity
- Mostly clear
- Newbie friendliness
- 48/100