Moisture Model Checks

[AMLattanzi]

Summary

This PR expands the logical checks for moisture models. Key areas of the code check if a moisture model that has cold components was utilized. New moisture models have been added but these logical checks were not updated. This PR corrects that.

[asalmgren]

Code review

The diff is correct and can merge as-is. It also does a bit more than the summary claims: converting RadiationSimple from an enumerated model list to moisture_indices.qi fixes a latent bug — the old list set m_ice = true for Kessler, whose RhoQ3_comp is rain, not ice (Source/DataStructs/ERF_DataStruct.H:343-349). Kessler runs were folding rain water into the longwave optical depth.

Verified

• SolverChoice::moisture_indices.qi is an absolute conserved-array component (RhoQ3_comp etc.), assigned in ERF_DataStruct.H:288-359, so cons_arr(i,j,k,qi_comp) is indexed correctly. m_ice = (m_qi_comp >= 0) guards the read, so qi_comp == -1 is never dereferenced.
• The newly ice-enabled schemes are all in bounds: WSM6 (n_qstate_moist_size = 6), WDM6 (9), SuperDroplets (6). RhoQ3_comp is allocated and initialized in each (ERF_SuperDropletsMoistInit.cpp:345, ERF_UpdateWSM6.cpp, ERF_UpdateWDM6.cpp).
• Every MoistureType enumerator is covered by the moisture_indices if-chain, and MoistureComponentIndices() = default exists (ERF_DataStruct.H:216), so MoistureType::None is safe.
• The ERF_Tagging.cpp addition is right: erf_dermaxreflectivity (ERF_Derive.cpp:580) hardcodes RhoQ4/5/6 as rain/snow/graupel, and WSM6/WDM6 both use exactly that layout. The warm and no-ice schemes are correctly still excluded — they don't allocate RhoQ5/RhoQ6. The new list also brings tagging into agreement with the pre-existing gates at ERF_Plotfile.cpp:552 and :565.
• Regression impact: only Exec/CanonicalTests/DYCOMS2RF01/inputs_DYCOMS2RF01 sets erf.radiation_model = "Simple" (with SuperDroplets), and it is not in Tests/CTestList.cmake. No gold files affected.

Medium — the same stale check survives in WriteSubvolume, and there it's out of bounds

Source/IO/ERF_WriteSubvolume.cpp:73

This is the identical pattern the PR is fixing, but it was missed. The else fall-through is commented "For moisture_type SAM and Morrison we have all six variables" — and that branch is now reached by MoistNoCondensation, which has only two moist state variables.

Failure scenario: moisture_model = MoistNoCondensation with erf.subvol_vars = rhoQ3 (or rhoQ4..rhoQ6). The name passes the filter, and the unguarded

mf.ParallelCopy(vars_new[lev][Vars::cons], i, ...);

at line 249 then reads source component 6 from a MultiFab with ncomp = NDRY + NSCALARS + 2 = 6 — out of bounds.

The same else also permits rhoQ7..rhoQ11 for every model except Morrison (11) and WDM6 (9), so moisture_model = SAM + rhoQ7 reads component 10 of a 10-component array. The derived-name chain at line 133 has the identical hole: qi / qsnow / qgraup are accepted for MoistNoCondensation.

Driving this off solverChoice.moisture_indices, exactly as the PR does in RadiationSimple, would close it.

Low — the tagging check is still an enumerated list, and still drifts

Source/Refinement/ERF_Tagging.cpp:302

The new list is correct today but remains the wrong shape: it enumerates models where it means to test data layout. SuperDroplets has qr/qs/qg at RhoQ4/5/6 (ERF_DataStruct.H:350-358) and would compute max reflectivity correctly, yet a run with moisture_model = SuperDroplets and a max_reflectivity refinement indicator hits the Abort at line 317.

Testing mi.qr >= 0 && mi.qs >= 0 && mi.qg >= 0 instead would be self-maintaining, and would keep this in sync with the two copies in ERF_Plotfile.cpp that currently have to be edited in lockstep every time a scheme is added.

Low — RadiationSimple is only half converted to the index map

Source/PhysicsInterfaces/Radiation/Simple/ERF_RadiationSimple.cpp:85

qi now comes from moisture_indices.qi, but qv and qc on the two lines above still hardcode RhoQ1_comp / RhoQ2_comp and are gated on the coarse m_moist flag. Harmless today — every scheme maps qv→RhoQ1 and qc→RhoQ2, and every moist scheme has at least two moist variables — but it leaves in place the exact trap the PR just removed, for the next scheme whose layout differs. Caching m_qv_comp / m_qc_comp alongside m_qi_comp would make the whole kernel index-map-driven.

Pre-existing, not a finding against this PR

The max_reflectivity tagging branch iterates with mfi.tilebox() under TileNoZ(), while erf_dermaxreflectivity asserts spans_full_column — that assert fires if the grid is decomposed in z. This predates the PR and applies equally to Morrison and SAM, though widening the allowed model set does expose it to more runs.

[AMLattanzi]

@asalmgren Updated by added helper routines to the MoistureIndices struct so that we can test on the indices rather than the models everywhere. This requires us to just add the moisture model indices in datastruct rather than all over the code as we move forward.

[asalmgren]

Review of the moisture-model check refactor. The direction here is right: replacing the hardcoded moisture_type == X || ... chains with queries against SolverChoice::moisture_indices is the correct move, and it fixes two genuine pre-existing bugs along the way — the old setSubVolVariables derived loop iterated derived_subvol_names.size() while indexing derived_names[i] (it only worked because the first five entries happened to line up), and the old cons filter would have let rhoQ10/rhoQ11 through for WDM6 and for MoistureType::None, which is a real out-of-bounds ParallelCopy.

Four comments inline. The two I'd want resolved before merge are the has_comp() filter in ERF_WriteSubvolume.cpp (silently drops components that are allocated and written) and the newly-reachable mucape/precipitable subvolume paths (they can abort the run on an ordinary z-decomposed grid).

[asalmgren]

This asks the wrong question. has_comp() tests whether the index map names a component, but the guard we need is whether the conserved state allocates it — and those two differ.

Morrison_NoIce is mapped with only qv/qc/qr (ERF_DataStruct.H:403-411), so qi, qs, qg, nc, nr, ns, ng are all -1. But it runs on the Morrison class, which has n_qstate_moist_size = 11, and ERF_UpdateMorrison.cpp:49-59 writes RhoQ3_comp, RhoQ5_comp, RhoQ6_comp, RhoQ7_comp (nc) and RhoQ9_comp (nr) unconditionally every step. So erf.subvol.vars = rhoQ7 now silently produces no output where it previously did. The same applies to the SuperDroplets non-moist species, which live at RhoQ1_comp + m_qstate_moist_size + 2k (i.e. RhoQ7_comp onward) and are all -1 in the index map.

The 3D-plotfile path answers this identical question with an allocation-based test — cons_comp < caps.conserved_state_size (ERF_PlotfileSelection.H:315-316). That is equally safe against the out-of-bounds ParallelCopy this filter is guarding against, and it keeps the two output paths consistent.

[AMLattanzi]

This is an interesting design choice that we need to make and be consistent about. Allocating all the variables a moisture model could have allows the class to be reused with limited functionality and no remapping of the indices. However, writing of data that is allocated but not updated seems like an odd choice to me. The moisture model may recopy a state variable back to the dycore but if it did not modify said variable do we want to allow users to write the essentially dead data? We should discuss before collapsing one pathway onto the other.

[asalmgren]

Adding "precipitable" and "mucape" here makes the calculate_derived calls at ERF_WriteSubvolume.cpp:288-289 reachable for the first time — previously those names could never enter subvol_var_names, so the calls were dead.

Both kernels carry AMREX_ALWAYS_ASSERT_WITH_MESSAGE(spans_full_column(bx,geomdata), ...) (ERF_Derive.cpp:735 and :785), which is active in optimized builds. WriteSubvolume builds dmf on grids[lev_for_sub] and iterates with TileNoZ(), which prevents tiling in z but does not undo a z-decomposition of the grids themselves — and ERF decomposes in z whenever max_grid_size[0][2] < domain.length(2) (ERF_MakeNewLevel.cpp:39-40), plus always for lev > 0 and on restart.

So erf.subvol.vars = mucape with a routine amr.max_grid_size = 32 and nz > 32 will abort at the first subvolume write, and nothing at name-selection time warns about it. The 3D plotfile has the same latent constraint, so this may be an accepted project-wide rule — but it becomes newly reachable here.

[asalmgren]

Worth calling out that this is an answer change for two groups of models, and no test or gold file was touched.

Previously m_ice was true for Kessler and the kernel read RhoQ3_comp as "qi" — but for Kessler RhoQ3 is rain (ERF_DataStruct.H:343-348), so the old optical depth xk*rho*(qc+qi)*dz was effectively including rain water. Kessler maps qi = -1, so Kessler + erf.radiation_model = Simple now loses that term entirely. In the other direction, WSM6/WDM6/SuperDroplets gain an ice term they never had; the only in-tree Simple-radiation input, Exec/CanonicalTests/DYCOMS2RF01/inputs_DYCOMS2RF01:75,80, is SuperDroplets, which does maintain a real q_i in RhoQ3 (ERF_SuperDropletsMoistUtils.cpp:306).

The new mapping is the physically correct one — this is not a request to change it, just a note that it silently moves numbers and should probably be flagged in the PR description (and covered by a regression test) rather than landing unannounced.

[asalmgren]

Two related points on this helper.

The "use directly as a filter" contract is weaker than the comment suggests. The doc says names backed by the qmoist arrays return true so the predicate can filter a mixed name list, but availability of those arrays is itself scheme-dependent, and this doesn't check it. The existing tested predicate does: rel_humidity requires moisture.relative_humidity_diagnostic && (SatAdj || qmoist_size >= 6) and condensation_rate requires qmoist_size >= 4 (ERF_PlotfileSelection.H:376-378). Here both return true unconditionally, including for MoistureType::None. The unknown-name fall-through to return true has the same shape. That's also why the caller still needs the ad-hoc ok_to_add &= (mi.has_moisture() || derived_names[i] != "mucape") at ERF_WriteSubvolume.cpp:94 — so two different mechanisms now gate moisture availability inside one loop, while the sibling "precipitable" is gated inside has_derived_var.

This duplicates logic that already exists and is unit-tested. erf_plotfile::plot3d_fixed_variable_available / plot3d_moist_q_component_available / plot3d_precipitating_q_range (ERF_PlotfileSelection.H:309-378) answer exactly "is this moist variable available for this scheme", are covered by Tests/Unit/IO/ERF_GTestPlotfileSelection.cpp, and are already what setPlotVariables calls. The two copies disagree today — "qp" here is qr||qs||qg, while the selection module requires moisture.rain plus an available precipitating q-range — so the same inputs file can select qp for the plotfile but not the subvolume. They will drift further as schemes are added.

Routing setSubVolVariables through the existing predicates would collapse this duplication and also resolve the has_comp divergence I flagged in ERF_WriteSubvolume.cpp.