Support reference pressure in anelastic microphysics

[pressel]

Summary

This PR adds consistent reference-pressure thermodynamics for supported Eulerian microphysics schemes when ERF is run with anelastic dynamics.

On an anelastic AMR level, ERF uses

$$\theta = \frac{\rho\theta}{\rho}, \qquad p_{\mathrm{micro}} = p_{\mathrm{base}},$$

and diagnoses temperature as

$$T = \theta \left( \frac{p_{\mathrm{base}}}{p_{\mathrm{ref}}} \right)^{R_d/c_p}.$$

Here, (p_{\mathrm{base}}) is ERF's hydrostatic base-state pressure and (p_{\mathrm{ref}}) is the fixed reference pressure used in the definition of potential temperature.

The anelastic pressure-temperature conversion uses ERF's configured erf.c_p.

Compressible AMR levels retain their existing local equation-of-state behavior.

Supported Eulerian microphysics

The reference-pressure treatment applies to:

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

SatAdj already uses the hydrostatic base-state pressure in anelastic configurations.

MoistNoCondensation does not have a pressure-dependent condensation source requiring this treatment.

Pressure conventions

ERF handles the scheme-specific pressure units internally:

• Kessler and SAM use hPa/mbar internally.
• Morrison, WSM6, and WDM6 use Pa.

Users do not provide an additional microphysics pressure field or perform these pressure-unit conversions themselves.

Compressible non-regression

The existing compressible thermodynamic path is intentionally preserved.

In particular, the historical WSM6 and WDM6 temperature-to-potential-temperature conversion behavior is retained in compressible configurations, including when erf.c_p differs from the default value.

For anelastic levels, those conversions instead use the configured (R_d/c_p) together with the held hydrostatic base-state pressure.

SAM

The existing compressible SAM coefficient construction is unchanged:

1. horizontally average density, potential temperature, and water-vapor mixing ratio;
2. diagnose pressure and temperature from those averaged state variables.

The anelastic path instead uses thermodynamic fields diagnosed from the hydrostatic base-state pressure.

Mixed AMR hierarchies

erf.anelastic may be specified per AMR level.

The microphysics thermodynamic path is selected independently on each level:

• anelastic levels use the hydrostatic base-state pressure;
• compressible levels use the existing local equation of state.

SuperDroplets

SuperDroplets is not currently supported if any AMR level is anelastic.

The current SuperDroplets implementation does not consume ERF's hydrostatic base-state pressure. Instead, it diagnoses pressure and temperature through the compressible equation of state and uses those thermodynamic fields in saturation calculations, phase change, particle mass change, terminal-velocity calculations, and coalescence.

That pressure diagnosis is not valid for anelastic dynamics, so ERF now detects and rejects this configuration during input validation rather than allowing physically inconsistent microphysics.

Users should either:

• run SuperDroplets with compressible dynamics; or
• choose one of the supported Eulerian microphysics schemes for an anelastic simulation.

This PR does not modify SuperDroplets particle physics.

Scope

This rebuilt PR intentionally does not:

• introduce pi0 as a new independent microphysics input;
• treat stored pi0 as authoritative for anelastic microphysics;
• reconstruct supposedly missing pi0 from legacy checkpoints;
• alter checkpoint/restart formats;
• refactor hydrostatic base-state initialization;
• modify Cloud Chamber implementation;
• regenerate TKE regression gold files;
• modify SuperDroplets physics.

Documentation

The Sphinx documentation now describes:

• the anelastic microphysics thermodynamic contract;
• the relationship among potential temperature, hydrostatic base-state pressure, and temperature;
• supported Eulerian schemes;
• scheme pressure units;
• mixed compressible/anelastic AMR hierarchies;
• the current SuperDroplets incompatibility.

Validation

The rebuilt implementation was validated with:

• focused anelastic microphysics wiring tests: 14/14;
• full serial unit suite: 461/461;
• parallel unit suite: 16/16;
• Cloud Chamber and ABL regression suite: 10/10;
• a Spack/MPI AppleClang build;
• git diff --check;
• strict Sphinx builds against both clean development and the feature branch.

The Sphinx feature build introduces no new documentation warnings relative to clean development.

Focused regression coverage includes:

• Kessler reference-pressure diagnosis;
• Morrison reference-pressure diagnosis;
• SAM anelastic and compressible heterogeneous-plane behavior;
• WSM6 anelastic pressure handling;
• WSM6 custom-c_p compressible non-regression;
• WDM6 anelastic temperature-to-potential-temperature conversion;
• WDM6 custom-c_p compressible non-regression;
• SuperDroplets/anelastic configuration rejection.

GitHub CI results should be interpreted from the current PR checks rather than from this description.

[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.