model_eqns = 1 (Gamma/Pi-infinity) is unusable in simulation: no gamma-law path in conversion or Riemann layers

[sbryngelson]

Summary

model_eqns = 1 (the Γ/Π∞ model) is documented and accepted by the case validator, and pre_process / post_process implement it, but simulation has no working path for it. A gamma-law case run through simulation computes rho = gamma = pi_inf = 0 in the conversion and Riemann layers.

This is pre-existing on master. Filing after an audit of the whole branch; nothing here is introduced by an open PR.

Why it fails

For model_eqns = 1 the state vector carries gamma and pi_inf as transported variables (m_global_parameters_common.fpp:115):

eqn_idx%cont%beg = 1;  eqn_idx%cont%end = 1        ! single density
eqn_idx%mom%beg  = 2;  eqn_idx%mom%end  = 1 + num_vels
eqn_idx%E        = mom%end + 1
eqn_idx%adv%beg  = E + 1                            ! = eqn_idx%gamma
eqn_idx%adv%end  = adv%beg + 1                      ! = eqn_idx%pi_inf
sys_size         = adv%end

Correct conversion is therefore a direct read (s_convert_mixture_to_mixture_variables), not the volume-fraction sum used by the other models.

case_validator.py:355 prohibits setting num_fluids when model_eqns == 1, so a valid case never supplies it. pre_process and post_process compensate:

• src/pre_process/m_global_parameters.fpp:441 — if (model_eqns == model_eqns_gamma_law) num_fluids = 1
• src/post_process/m_global_parameters.fpp:326 — same

simulation has no equivalent, so num_fluids stays at dflt_int = -100. Consequences in simulation:

1. s_compute_species_fraction declares alpha_rho_K(num_fluids) / alpha_K(num_fluids) — zero-size arrays.
2. Every do i = 1, num_fluids accumulation loop is empty, so rho_K, gamma_K, pi_inf_K, qv_K stay 0.
3. s_accumulate_mixture_properties (m_riemann_state.fpp:1015, called from m_riemann_solver_hll.fpp:202-203 and m_riemann_solver_lf.fpp:197-198) returns zeros for both Riemann states.

Forcing num_fluids = 1 alone is necessary but not sufficient. With it, s_compute_species_fraction takes its num_fluids == 1 branch and sets alpha_K(1) = q_vf(eqn_idx%adv%beg) — which for gamma-law is the transported gamma, not a volume fraction. The kernel then computes gamma_K = alpha_K(1)*gammas(1) and pi_inf_K = alpha_K(1)*pi_infs(1), i.e. the transported gamma scaled by fluid_pp-derived constants.

What is missing in simulation

• num_fluids = 1 for gamma-law — the one-line gap that pre_process and post_process both close.
• A gamma-law branch in the conversion layer. s_convert_conservative_to_primitive_variables (m_variables_conversion.fpp:~473) guards only model_eqns /= model_eqns_4eq, so gamma-law reaches the species kernel. s_convert_species_to_mixture_variables_kernel has no gamma-law case, and it is a GPU_ROUTINE, so the device path needs one too.
• A gamma-law branch in s_convert_primitive_to_flux_variables, which also calls the species kernel unconditionally.
• A gamma-law branch in the Riemann states. s_accumulate_mixture_properties must read gamma/pi_inf from the state vector instead of summing over fluids.
• Any test or example. No case in tests/ or examples/ sets model_eqns = 1, which is why this has never surfaced.

Partial support that already exists

The endpoints were clearly intended to work, which is what makes the gap easy to miss:

• src/simulation/m_global_parameters.fpp:778 — skips Re_idx setup for gamma-law
• src/simulation/m_cbc.fpp:836 — has an explicit gamma-law case
• src/pre_process/m_check_patches.fpp:476-480 — validates patch rho/gamma/pi_inf for gamma-law
• src/post_process/{m_start_up,m_data_output,m_derived_variables}.fpp — write gamma/pi_inf output
• src/common/m_variables_conversion.fpp:62 — the host dispatcher already routes gamma-law correctly

So pre_process can build a gamma-law initial condition and post_process can render it; only the solver in between is missing.

Suggested resolution

Either:

(a) Implement it — add num_fluids = 1 in simulation, plus gamma-law cases in the conversion kernel (host and device), s_convert_primitive_to_flux_variables, and s_accumulate_mixture_properties; add a 1D regression case.

(b) Reject it — have case_validator.py prohibit model_eqns = 1 for simulation, and mark it pre/post-only in docs/documentation/case.md:492,567-568, which currently presents it as a supported model on equal footing with 2 and 3.

Option (b) is much cheaper and matches the current reality; option (a) is only worth it if the Γ/Π∞ solver is actually wanted.

Note for reviewers of #1679

That PR changes the conversion guard from #ifdef MFC_SIMULATION to #ifdef MFC_GPU, so a CPU simulation build now reaches the correct host dispatcher for gamma-law while a GPU build still takes the species kernel. It does not introduce the bug, but it does make it asymmetric between CPU and GPU, which is worth resolving deliberately.

[sbryngelson]

What it would take to make model_eqns = 1 usable in simulation

Audited the full branch. The obstacle is structural, not just the missing num_fluids line: gamma-law has one fluid but two advected entries, so nothing that loops do i = 1, num_fluids over eqn_idx%adv can express it.

eqn_idx%adv%beg = eqn_idx%gamma     ! transported gamma
eqn_idx%adv%end = eqn_idx%pi_inf    ! transported pi_inf
num_fluids      = 1

Every conversion and Riemann path derives gamma/pi_inf by summing alpha_K(i)*gammas(i) over num_fluids. For gamma-law those two values must be read from the state vector, not derived. So the fix is to give each of those sites an explicit gamma-law branch rather than to generalize the loops.

Five changes, all additive.

1. src/simulation/m_global_parameters.fpp — force num_fluids = 1

In s_initialize_global_parameters_module, before s_initialize_eqn_idx, matching pre_process/m_global_parameters.fpp:441 and post_process/m_global_parameters.fpp:326:

! Gamma/Pi_inf: single fluid (matches pre_process and post_process)
if (model_eqns == model_eqns_gamma_law) num_fluids = 1

Necessary but not sufficient on its own — without the changes below, alpha_K(1) still holds the transported gamma and the kernel computes gamma_K = gamma*gammas(1).

2. src/common/m_variables_conversion.fpp — gamma-law branch in cons -> prim

In s_convert_conservative_to_primitive_variables, replace the outer guard:

                    if (model_eqns == model_eqns_gamma_law) then
                        rho_K = qK_cons_vf(1)%sf(j, k, l)
                        gamma_K = qK_cons_vf(eqn_idx%gamma)%sf(j, k, l)
                        pi_inf_K = qK_cons_vf(eqn_idx%pi_inf)%sf(j, k, l)
                        qv_K = 0._wp
                    else if (model_eqns /= model_eqns_4eq) then
                        ! ... existing MFC_GPU / host branches unchanged ...
                    end if

These are plain reads, valid on host and device alike, so this branch also removes the CPU/GPU divergence noted at the end of the issue: both builds would take the same path for gamma-law.

3. src/common/m_variables_conversion.fpp — gamma-law branch in prim -> flux

In s_convert_primitive_to_flux_variables, ahead of the hypoelasticity dispatch (note this array is 4-D, indices last):

                    if (model_eqns == model_eqns_gamma_law) then
                        rho_K = qK_prim_vf(j, k, l, 1)
                        gamma_K = qK_prim_vf(j, k, l, eqn_idx%gamma)
                        pi_inf_K = qK_prim_vf(j, k, l, eqn_idx%pi_inf)
                        qv_K = 0._wp
                    else if (hypoelasticity) then
                        ! ... existing ...

4. The three Riemann solvers — read gamma/pi_inf instead of accumulating

m_riemann_solver_hll.fpp, m_riemann_solver_hllc.fpp, m_riemann_solver_lf.fpp all call s_accumulate_mixture_properties(num_fluids, ...) (m_riemann_state.fpp:1015), which sums over fluids. Guard those calls. qL_prim_rsx_vf / qR_prim_rsx_vf are single reshaped buffers reused for all three directions, and ${SF('')}$ / ${SF(' + 1')}$ are the existing index expressions from the enclosing #:for:

                            if (model_eqns == model_eqns_gamma_law) then
                                rho_L = alpha_rho_L(1)
                                rho_R = alpha_rho_R(1)
                                gamma_L = qL_prim_rsx_vf(${SF('')}$, eqn_idx%gamma)
                                gamma_R = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%gamma)
                                pi_inf_L = qL_prim_rsx_vf(${SF('')}$, eqn_idx%pi_inf)
                                pi_inf_R = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%pi_inf)
                                qv_L = 0._wp
                                qv_R = 0._wp
                            else
                                call s_accumulate_mixture_properties(num_fluids, alpha_rho_L, alpha_L, rho_L, gamma_L, pi_inf_L, &
                                                                     & qv_L)
                                call s_accumulate_mixture_properties(num_fluids, alpha_rho_R, alpha_R, rho_R, gamma_R, pi_inf_R, &
                                                                     & qv_R)
                            end if

alpha_rho_L(1) is already correct: the preceding loop runs do i = 1, eqn_idx%cont%end, and cont%end == 1 for gamma-law.

The mpp_lim clamp block just above needs no guard — case_validator.py:359 already prohibits mpp_lim with model_eqns = 1.

5. Regression coverage

There is currently no test or example with model_eqns = 1, which is why this was never caught. A 1D shock tube in toolchain/mfc/test/cases.py would cover conversion, reconstruction, the Riemann solver, and CBC in one case. Worth adding an examples/1D_gammalaw_shocktube/ too, since pre_process and post_process already support the model end to end.

Needs no change

• s_convert_primitive_to_conservative_variables — copies do i = eqn_idx%adv%beg, eqn_idx%adv%end verbatim, which is already right for gamma-law (this is why pre_process works).
• m_cbc.fpp:836 — already has an explicit gamma-law case.
• m_global_parameters.fpp:778 — already skips Re_idx setup for gamma-law.
• s_compute_speed_of_sound — takes gamma/pi_inf as arguments, so it is correct once they are.
• s_compute_species_fraction — not reached under the new branches.
• m_riemann_solver_hlld.fpp — MHD only; does not call s_accumulate_mixture_properties.

Scope and risk

Roughly 40 added lines across 5 files, all inside new model_eqns == model_eqns_gamma_law branches, so every existing model keeps its exact code path — the golden files should not move. The real cost is validation: gamma-law has no reference results in the repo today, so the new test's golden file has to be justified physically (e.g. against an analytic Sod solution) rather than simply recorded.

If the Γ/Π∞ solver is not actually wanted, option (b) in the issue body — reject model_eqns = 1 for simulation in the validator and mark it pre/post-only in docs/documentation/case.md:492,567-568 — remains far cheaper and is the honest description of the current state.