Support reference-state thermodynamics in anelastic microphysics

[pressel]

Summary

This PR adds consistent pressure and temperature handling for Eulerian microphysics when ERF is run with anelastic dynamics.

For anelastic levels, the Kessler, SAM, and Morrison model families, along with WSM6, now use ERF's hydrostatic reference state:

$$\theta = \frac{\rho\theta}{\rho}, \qquad T = \theta \pi_0, \qquad p = p_0$$

Here, ($p_0$) is the hydrostatic reference pressure and ($\pi_0$) is the corresponding reference Exner function.

This avoids reconstructing microphysics pressure from the compressible equation of state on anelastic levels while preserving the existing compressible behavior.

Changes

Reference-state thermodynamics

• Add the reference Exner function, ($\pi_0$), to the ERF base-state thermodynamic contract.
• Initialize reference pressure and Exner consistently through the shared hydrostatic base-state path.
• Treat the stored reference Exner function as authoritative when diagnosing anelastic temperature.
• Reconstruct ($\pi_0$) when reading legacy checkpoints that contain ($p_0$) but predate the new base-state component.
• Preserve base-state consistency through initialization, AMR level construction, ghost filling, and restart handling.

Microphysics integration

Add a shared thermodynamic diagnosis for Eulerian microphysics and use it for the following anelastic model selections:

• Kessler
• Kessler_NoRain
• SAM
• SAM_NoIce
• SAM_NoPrecip_NoIce
• Morrison
• Morrison_NoIce
• WSM6

For these schemes, anelastic microphysics uses:

$$T = \theta \pi_0$$

and:

$$p = p_0$$

Scheme-specific pressure conventions remain unchanged:

• Kessler and SAM store pressure internally in mbar/hPa.
• Morrison and WSM6 retain pressure in Pa.
• ERF performs these conversions internally; users do not need to provide scheme-specific pressure units.
• WSM6 continues to reconstruct conserved potential temperature from its updated temperature and held pressure after the microphysics source update.

Compressible configurations continue to use the existing local equation-of-state diagnosis for pressure and temperature.

SuperDroplets compatibility

SuperDroplets is not currently supported when any AMR level uses anelastic dynamics.

ERF now rejects this configuration early because the current SuperDroplets thermodynamic path reconstructs pressure through the compressible equation of state.

Users can instead:

• run SuperDroplets with compressible dynamics; or
• select a supported Eulerian moisture model for anelastic simulations.

The existing particle-enabled build requirement for SuperDroplets remains unchanged.

Cloud Chamber behavior

• Strengthen the Cloud Chamber initialization path for the prescribed reference thermodynamic state.
• Preserve independently prescribed Cloud Chamber density.
• Use the configured reference temperature and pressure consistently.
• Retain the analytic hydrostatic pressure profile.
• Improve validation of incompatible Cloud Chamber configurations.
• Strengthen deterministic pressure and configuration checks.

Restart compatibility

• Add support for reconstructing the reference Exner function when reading older checkpoints that contain the previous base-state component count.
• Preserve current checkpoint behavior when (\pi_0) is already present.
• Avoid changing unrelated checkpoint formats or restart semantics.

Documentation

Update the Sphinx microphysics documentation to explain:

• how anelastic dynamics are enabled;
• which Eulerian microphysics models use the hydrostatic reference state;
• how pressure and temperature are supplied to those schemes;
• that users do not need to provide an additional microphysics pressure field;
• that users do not need to perform pressure-unit conversions;
• that compressible microphysics behavior is unchanged;
• how legacy checkpoints reconstruct a missing reference Exner function;
• why SuperDroplets is not supported with anelastic dynamics.

The SuperDroplets compatibility warning and particle build requirement are retained in the dedicated Super-Droplet documentation section.

User impact

Users enable anelastic dynamics through the existing input:

erf.anelastic = 1

The option may be specified once for all AMR levels or as one value per level.

No additional microphysics pressure field is required. ERF supplies the hydrostatic reference pressure and corresponding temperature to supported Eulerian schemes and handles internal pressure-unit conversions automatically.

Existing compressible configurations are unchanged.

Testing

This PR adds or strengthens coverage for:

• shared anelastic microphysics thermodynamic diagnosis;
• production Array4 and MultiFab wiring;
• Kessler reference-state copy-in;
• SAM reference-state copy-in and pressure handling;
• Morrison reference-state copy-in;
• WSM6 reference-state copy-in and conserved-state copy-out;
• hydrostatic reference pressure and Exner consistency;
• internal and ghost-cell base-state handling;
• legacy checkpoint reconstruction of missing (\pi_0);
• SatAdj temperature and pressure diagnostics;
• CUDA-portable SatAdj test initialization;
• Cloud Chamber thermodynamic checks;
• Cloud Chamber invalid-configuration failures;
• SuperDroplets and anelastic compatibility checks.

The ABL_MOST_IMP_DIFF_TKE regression reference was regenerated after integrating current upstream changes. The resulting differences are at roundoff scale, and no regression comparison tolerances were loosened.

Notes for reviewers

The stored reference Exner function is authoritative for anelastic temperature diagnosis. Microphysics does not recompute it from reference pressure inside each scheme kernel.

The compressible thermodynamic path is intentionally unchanged.

The updated TKE gold file reflects the current combined branch after upstream integration. The update does not loosen tolerances or intentionally change the physical configuration.

[asalmgren]

PR #3564 Review — Support reference-state thermodynamics in anelastic microphysics

Diff reviewed: git diff $(git merge-base pr3564 development) pr3564 — 31 files, +1495/−153.

Summary

The core change is correct. T = θ·π is exact for ERF's moist EOS (getTgivenRandRTh
reduces to θ·π with θ = ρθ/ρ dry), so T = θ·π₀ with p = p₀ is the right anelastic
diagnosis, and it matches the pre-existing SatAdj anelastic path. The erf_enforce_hse
refactor is bit-identical (getExnergivenP now takes the register value instead of the
just-stored memory value — exact for IEEE doubles). Both base_state[lev] and cons carry
ComputeGhostCells+1 ghosts, so the growntilebox() copy-in loops do not read out of
bounds. Scheme coverage is complete for the current MoistureType enum (SuperDroplets now
rejected, SatAdj already handled, MoistNoCondensation pressure-independent, WDM6 not
user-selectable).

Findings

1. Gold file regeneration will break CI — medium/high

Tests/ERFGoldFiles/ABL_MOST_IMP_DIFF_TKE/Level_0/Cell_H:14

Nothing in this diff can change a dry ABL+MOST+TKE run: the only production path it touches
for such a run is erf_enforce_hse, and that refactor is bit-identical. Yet the gold data
was regenerated. Worse, development already updated this same gold in #3900, and the PR's
values differ from development's by up to 1.75e-6 relative (comp 3:
9.52576489081652267e-01 vs 9.52578159161245441e-01), far above the "roundoff scale"
claimed in the PR body. Merging as-is will conflict with / regress #3900 and make
ABL_MOST_IMP_DIFF_TKE fail.

Fix: drop the gold change and rebase.

2. Copy-out breaks the "stored π₀ is authoritative" invariant the PR introduces — medium

Source/Microphysics/WSM6/ERF_WSM6.H:106, Source/Microphysics/SAM/ERF_SAMUtils.H:178

Copy-in computes T = θ·π₀_stored; copy-out computes
θ = getThgivenTandP(T, p₀) = T·(p_0/p₀)^rdOcp, i.e. it recomputes Exner from p₀
instead of dividing by the stored π₀.

Concrete scenario: any base state where π₀ != getExnergivenP(p₀) — e.g. the
interpolated/BC-filled base state at AMR level boundaries and domain-exterior ghosts, or a
future user/reader-prescribed base state — makes each microphysics call apply
θ → θ·π₀_stored/π₀_eos, a systematic θ jump per step even with zero microphysical
tendency. The PR's own new scalar test
(ERF_GTestMicrophysicsUtilsScalar.cpp,
EXPECT_NE(low.temperature, getTgivenPandTh(p0_low, theta, RdoCp))) asserts precisely that
these two differ. Morrison has the same asymmetry (ERF_AdvanceMorrison.cpp:257 builds
pii = pow(pres/p0, rdcp)).

Fix: use the stored π₀ on both sides, or drop the "authoritative π₀" contract and
derive π₀ from p₀ everywhere as SatAdj already does.

3. Unintended compressible behavior change in SAM + new dead state — medium

Source/Microphysics/SAM/ERF_InitSAM.cpp:178-211

Compute_Coefficients switched from tabs1d(k) = T(⟨ρ⟩, ⟨ρ⟩⟨θ⟩, ⟨q_v⟩) to
tabs1d(k) = ⟨T⟩ (and pres1d likewise). These differ at O(perturbation variance), not
roundoff, and tabs1d feeds sam_compute_coefficient_row →
accrsi/coefice/evaps*/evapg*/evapr*. So compressible SAM results change, contradicting
"Compressible configurations continue to use the existing local equation-of-state
diagnosis"; there is no SAM regression test in CTestList.cmake to catch it
(SquallLine_2D is Kessler).

Separately, pres1d is now written but read nowhere in the tree — the plane-average of
pres and its host/device copies are pure overhead.

4. Dead code — advertised legacy-checkpoint feature does not exist — medium

Source/IO/ERF_Checkpoint.cpp:751, Source/IO/ERF_BaseStateRestart.H:19

BaseState::pi0_comp == 2 (it has existed since #1908 — this PR does not add it), while
ncomp_base_to_read is initialized to 3 at ERF_Checkpoint.cpp:768 and only ever read
from the header as 3/4/5. So ncomp_base_to_read <= BaseState::pi0_comp is never true and
the whole reconstruction block, the new header field, and the "legacy checkpoints
reconstruct a missing π₀" note in Microphysics.rst are unreachable. The PR's own unit test
has to pass legacy_ncomp = BaseState::pi0_comp — a value no real checkpoint can carry — to
exercise it.

Fix: either remove the block (and the docs note) or fix the intended predicate.

5. Latched mode flag consumed by a different call — low/medium

Source/Microphysics/WSM6/ERF_InitWSM6.cpp:48

m_use_anelastic_reference_pressure is set in Copy_State_to_Micro and read much later in
Copy_Micro_to_State (ERF_UpdateWSM6.cpp:9,20). The still-present single-arg entry points
(WSM6::Copy_State_to_Micro(cons) at ERF_InitWSM6.cpp:39,
WSM6::Update_Micro_Vars(cons) at ERF_WSM6.H:196, both overrides reachable through the
NullMoist interface) forward nullptr and silently reset the flag to false. Any future
caller landing between copy-in and copy-out on an anelastic level would make copy-out take
the compressible getThgivenRandT branch and write a wrong RhoTheta, with no diagnostic.
Kessler/SAM/Morrison avoid this by threading the mode through the call; WSM6 should too.

Relatedly, the guard at ERF_UpdateWSM6.cpp:9 can never fire: m_rdOcp defaults to
RdoCp and Define sets R_d/c_p, both > 0.

6. Test oracle hardcodes deck values it does not read — low

Source/Prob/ERF_CloudChamber.H:25

reference_base_temperature = 292.0 and prescribed_reference_density() = p_0/(R_d*292)
duplicate prob.T_0 and use the global p_0 rather than prob.p_inf, even though
Exec/ERF_Prob.cpp:36 derives rho_0 = p_inf/(R_d*T_0). It happens to be bit-exact for all
five current decks (p_inf = 100000.0, T_0 = 292.0), but any deck that changes
prob.T_0/prob.p_inf, or supplies prob.rho_0 directly, makes CloudChamberCheck fail at
Tests/CloudChamberCheck.cpp:417 with a "model" error that is really a stale test constant.

Fix: read the values from the inputs file, or assert them.

7. Unguarded decode of a sentinel — low

Tests/CloudChamberCheck.cpp:505

If the tie-match at lines 488-501 finds no cell (excess > 0 && |excess - max| <= 8·eps·max
can in principle miss), worst_code stays numeric_limits<Long>::max() and
static_cast<int>(worst_code / plane) is an out-of-range narrowing; the second reduction
then returns -DBL_MAX for every field and the failure message prints garbage instead of the
offending cell.

Fix: guard with if (worst_code == numeric_limits<Long>::max()) before decoding.

8. Dead member — low

Source/Microphysics/Kessler/ERF_Kessler.H:123,254

m_rdOcp is assigned in Define but never read anywhere in Kessler (unlike SAM/Morrison/
WSM6, which all use theirs). It is also declared without an initializer, so it would be
indeterminate if Define were ever skipped.