Skip to content

Make the Kessler scheme conserve water and θˡⁱ on any vertical grid - #1030

Open
glwagner wants to merge 5 commits into
mainfrom
fix/kessler-conservation-validation
Open

glwagner wants to merge 5 commits into
mainfrom
fix/kessler-conservation-validation

Conversation

@glwagner

@glwagner glwagner commented Sep 25, 2026 •

Copy link
Copy Markdown
Member

This PR fixes three conservation issues in Breeze’s DCMIP2016 Kessler implementation. Rain sedimentation now uses finite-volume cell thicknesses and includes the top cell in the CFL bound. It transports rain partial density directly, leaving vapor and cloud water unchanged and making the reported surface precipitation flux consistent with the water removed from the column.

Condensation and evaporation now preserve the liquid–ice potential temperature used by the model, with a temperature response and saturation-adjustment slope consistent with that formulation. Sedimentation retains the existing fixed-temperature convention, and accretion uses the post-sedimentation rain. Unsupported atmospheric use with StaticEnergyFormulation now produces an explicit error.

glwagner and others added 3 commits September 25, 2026 18:26
Three coupled defects of the DCMIP2016 Kessler kernel, found by the SCM
calibration audits, are fixed together because they share the same
bookkeeping:

1. Sedimentation geometry: the upwind rain-flux divergence divided by the
   distance between cell centers (and a half cell at the top) instead of the
   finite-volume cell thickness, so the column rain budget did not close on
   stretched grids (every junction cell) and the top cell lost rain at twice
   the correct rate. Fluxes ρqʳ𝕎ʳ are now divided by Δzᶜᶜᶜ, with zero inflow
   at the top and the CFL bound on every cell's thickness.

2. Density basis: on the anelastic core the kernel moved ρᵣ rʳ 𝕎ʳ in mixing
   ratio space at fixed total density and converted back with the changed rᵗ,
   rescaling vapor and cloud in every cell rain passed through and reporting
   a surface flux that was not the rain removed. Sedimentation now updates the
   prognostic partial density ρqʳ directly; the Kessler processes act on
   dry-air mixing ratios with ρᵈ = ρᵣ - ρᵗ (anelastic) or the prognostic ρᵈ
   (compressible), so the prognostic water budget closes to the substep-mean
   surface flux exactly.

3. Phase-change thermodynamics: T was incremented by ℒˡᵣ Δrˡ / cᵖᵈ and θˡⁱ
   re-derived from it, a percent-level inconsistency with the formulation's
   own θˡⁱ ↔ T relation that acted as a net θˡⁱ sink. Sedimentation is now
   applied at fixed T and phase change at fixed θˡⁱ, with the DCMIP single
   Newton step linearized by ∂T/∂rˡ at fixed invariant
   (`phase_change_temperature_slope`), matching `SaturationAdjustment`. The
   static-energy formulation is rejected with a clear error; the parcel path
   uses the same slope for its own invariant.

Tests: a new `dcmip2016_kessler_conservation.jl` (budgets on uniform,
junction and geometric grids, top cell, vapor untouched by sedimentation,
all-process water closure, θˡⁱ conservation for condensation, cloud and rain
evaporation across 1000-470 hPa, finite-difference check of the slope,
Float32 and Float64), and the column reference implementation in
`dcmip2016_kessler.jl` rewritten to the documented algorithm on a stretched
grid, with a compressible water-budget check.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…Kessler tests

- The constructor docstring now states that negative water partial densities are
  clipped to zero on entry (unchanged behaviour), that clipping creates the clipped
  mass, and that the closure holds for non-negative inputs only. A new test pins the
  water created by clipping to exactly the clipped mass.
- The conservation tests read the model temperature after `set!`; `update_state!`
  diagnoses it from the scheme's diagnostic mass fractions, which it refreshes only
  after reading them, so the helper now runs two passes before recording the
  pre-step temperature (this is what made the latent-heating-sign checks fail on CI).
- The reference-comparison test places rain in the bottom cell as well, so the
  surface flux it asserts on is actually exercised within one 10 s step.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…ce; retune the clipping test

- `test/dcmip2016_kessler.jl`: the DCMIP2016 Fortran translation that the previous
  "Julia vs Fortran" test used is restored verbatim as `dcmip2016_fortran_kessler!`
  (documented as the legacy algorithm, not the current one), together with its
  per-cell adjustment block `legacy_kessler_adjustment`. Two new test sets compare the
  current kernel with it only where the physics is unchanged: the process formulas
  (autoconversion, accretion, single-step saturation adjustment, rain evaporation)
  agree to 1e-14 when the current step is given the legacy latent-heating slope
  ℒ/cᵖᵈ, and the sedimentation of dry-air-carried rain on a uniform grid away from
  the top cell (compressible core, dry column, evaporation off) gives identical rain,
  θˡⁱ and surface flux. Whole-column equality is not claimed where the geometry,
  density basis or temperature update were deliberately corrected; those are
  covered by `kessler_column_reference!` and the conservation tests.
- `test/dcmip2016_kessler_conservation.jl`: the negative-input clipping test uses a
  dry raining column and a tolerance scaled by the inventory round-off, so the
  created mass is resolvable at Float32 (the previous relative tolerance on the
  created mass itself failed by 1e-4 on CI).

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
@glwagner
glwagner marked this pull request as ready for review September 25, 2026 20:16
@glwagner glwagner added the build all examples 🏗️ PRs which should build all examples label Sep 25, 2026
@codecov

codecov Bot commented Sep 25, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!

glwagner and others added 2 commits September 25, 2026 20:32
Documenter cannot resolve `@ref` links to `LiquidIcePotentialTemperatureState` and
`StaticEnergyState`, which have no docstrings, so the docs build failed; refer to them
as plain code instead.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…the clipped-relative bound

- `test/dcmip2016_kessler.jl`: a new test set measures, on the former "Physical
  fidelity" column (anelastic, uniform 100 m, all processes, one 10 s step), how far
  the corrected kernel departs from the legacy Fortran translation where the physics
  was deliberately corrected, and bounds it: vapor < 2e-3, cloud < 3e-2, rain < 1e-1
  of the species maxima; θˡⁱ by less than 2 % of the step's latent heating (measured
  0.56 %) and by more than 1e-3 K (so the deliberate departure cannot silently
  vanish); the legacy convention's water residual is nonzero while the current one is
  at round-off.
- `test/dcmip2016_kessler_conservation.jl`: the clipping test keeps the dry raining
  fixture but returns to the bound relative to the created mass
  (`≈ clipped rtol=budget_rtol(FT)`), which the dry fixture satisfies with margin
  (independent measurement: Float32 error 9.3e-9 against a 1e-6 bound); the injected
  negatives are doubled and the signal-to-round-off margin
  `clipped / (eps(FT) W₀ Nz) ≥ 1000` is asserted explicitly.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

build all examples 🏗️ PRs which should build all examples

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant