IGR integration

[StevenXue98]

Description

This PR integrates IGR with more of the shared RHS/Riemann infrastructure while preserving the existing IGR numerical behavior. This is an intermediate integration milestone, not the final IGR architecture. It intentionally preserves the current IGR reconstruction/Sigma/viscous numerics while moving RHS accumulation, RK update behavior, state layout, and LF Riemann dispatch toward shared MFC subroutines.

Main changes:

• Store the final volume fraction for IGR so the IGR state layout matches the normal model-equations layout more closely.
• Route IGR face states through the shared s_riemann_solver path using the LF Riemann solver.
• Use shared RHS accumulation helpers for IGR:
  • s_compute_advection_source_term
  • s_compute_additional_physics_rhs
• Split the old fused s_igr_riemann_solver behavior into smaller IGR-specific routines called from the shared RHS path:
  • s_igr_compute_jac_rhs
  • s_igr_reconstruct_cell_boundary_values
  • s_reconstruct_cell_boundary_values_visc_deriv
  • s_igr_correct_lf_fluxes
  • s_igr_sigma
• Reconstruct stored IGR velocity-gradient arrays into the shared LF-facing derivative arrays, so the LF solver constructs viscous source fluxes for IGR.
• Preserve legacy IGR behavior where needed:
  • skip the normal LF MPP limiter for IGR, since it was not part of the old IGR path
  • keep IGR’s fixed-coefficient viscous-derivative reconstruction so the shared LF viscous-source construction preserves legacy IGR numerics
• Remove the special padded rhs_vf allocation previously needed by IGR-specific direct RHS updates.

IGR-specific face reconstruction, Sigma correction, and viscous source-flux construction are still present in m_igr.fpp. Subroutine s_igr_riemann_solver is also still present in m_igr.fpp, which explains the increased line count.

Deferred follow-up work:

• Evaluate whether IGR can switch from reconstruct-conservative-then-convert to the normal convert-to-primitive-then-reconstruct order, i.e. using s_reconstruct_cell_boundary_values. This means switching from the current cons -> reconstruction -> prim order for IGR to the cons -> prim -> reconstruction order used in the standard path, which will likely cause numerical differences. Some preliminary experiments showed differences as large as ~1e-3 on some cases.
• evaluate whether IGR viscous derivatives can eventually be routed through s_get_viscous without changing numerics
• if numerically acceptable, route IGR face reconstruction and viscous derivatives through the standard path
• revisit compatibility with MHD and other currently prohibited physics

Type of change

• Refactor

Testing

Local testing:

./mfc.sh test --add-new-variables --only IGR
./mfc.sh test --add-new-variables --from F926448C --to F926448C
./mfc.sh test --add-new-variables --from 18A71AB2 --to 18A71AB2
./mfc.sh test --add-new-variables --from 43B5FEBD --to 43B5FEBD
./mfc.sh test --add-new-variables --from E49EF7B6 --to E49EF7B6
./mfc.sh test --only IGR

Additional local comparison/debug checks used during development:

python3 local_check_igr_full_alpha.py --default-igr
python3 local_igr_diagnostics/scripts/compare_igr_packs.py --default-igr

The IGR golden files were updated with --add-new-variables to include the newly stored final volume fraction. Existing golden entries were preserved, and the added final-alpha fields satisfy the expected volume-fraction sum checks, which is verified by the local Python scripts. The updated IGR tests match the previous conserved-variable outputs to machine-level roundoff, aside from the newly added final volume-fraction fields.

Checklist

Check these like this [x] to indicate which of the below applies.

• I added or updated tests for new behavior
• I updated documentation if user-facing behavior changed

See the developer guide for full coding standards.

GPU changes (expand if you modified src/simulation/)

• GPU results match CPU results
• Tested on NVIDIA GPU or AMD GPU
• Verified CPU only

AI code reviews

Reviews are not retriggered automatically. To request a review, comment on the PR:

• @claude full review — Claude full review (also triggers on PR open/reopen/ready)
• Or add label claude-full-review — Claude full review via label

[codecov]

Codecov Report

❌ Patch coverage is 82.44681% with 165 lines in your changes missing coverage. Please review.
✅ Project coverage is 57.89%. Comparing base (f4ab1aa) to head (267de82).
⚠️ Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
src/simulation/m_rhs.fpp	66.31%	30 Missing and 34 partials ⚠️
src/simulation/m_igr.fpp	91.60%	0 Missing and 57 partials ⚠️
src/simulation/m_time_steppers.fpp	31.11%	15 Missing and 16 partials ⚠️
src/simulation/m_riemann_state.fpp	11.11%	6 Missing and 2 partials ⚠️
src/simulation/m_riemann_solver_lf.fpp	0.00%	0 Missing and 3 partials ⚠️
src/simulation/m_riemann_solvers.fpp	0.00%	1 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1615      +/-   ##
==========================================
- Coverage   60.43%   57.89%   -2.55%     
==========================================
  Files          83       83              
  Lines       19871    20570     +699     
  Branches     2956     3009      +53     
==========================================
- Hits        12010    11908     -102     
- Misses       5860     6681     +821     
+ Partials     2001     1981      -20     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[StevenXue98]

I ran some local performance comparisons against the legacy IGR path.

Representative CPU results:

18A71AB2, 50 steps, 100 runs:

• runtime: 3.48s refactor vs 2.74s legacy (~1.27x)
• peak RSS: 71.6 MB refactor vs 57.4 MB legacy (~1.25x)

18A71AB2, 500 steps, 10 runs:

• runtime: 29.64s refactor vs 22.61s legacy (~1.31x)
• peak RSS: 71.6 MB refactor vs 57.3 MB legacy (~1.25x)

18A71AB2, larger grid, ~100^3, 20 steps, 3 runs:

• runtime: 58.54s refactor vs 55.31s legacy (~1.06x)
• peak RSS: 848.2 MB refactor vs 276.5 MB legacy (~3.07x)

RSS stands for Maximum resident set size, which is measured by the command:
/usr/bin/time -v -o "${time_file}" mpirun -np 1 "${exe}"

[github-actions]

Claude Code Review

Head SHA: 267de82

Files changed:

• 13
• src/common/include/3dHardcodedIC.fpp
• src/common/m_global_parameters_common.fpp
• src/common/m_variables_conversion.fpp
• src/pre_process/m_assign_variables.fpp
• src/simulation/m_igr.fpp
• src/simulation/m_rhs.fpp
• src/simulation/m_riemann_solver_lf.fpp
• src/simulation/m_riemann_solvers.fpp
• src/simulation/m_riemann_state.fpp
• src/simulation/m_start_up.fpp
  (+2 more: src/simulation/m_time_steppers.fpp, toolchain/mfc/case_validator.py, toolchain/mfc/test/cases.py)

Findings:

• src/simulation/m_igr.fpp — new s_igr_reconstruct_cell_boundary_values stores the L/R face states with the offset convention backwards, so the shared LF Riemann solver reads mismatched left/right states. Throughout the new subroutine (x/y/z blocks, e.g. qL_rsx_vf(j + 1, k, l, i) = alpha_rho_L(i) / qR_rsx_vf(j, k, l, i) = alpha_rho_R(i)), the reconstructed L value for the face computed at stencil index j is written to slot j + 1, while R is written unshifted to slot j. This is backwards relative to the codebase's established convention: src/simulation/m_riemann_solver_lf.fpp's s_lf_riemann_solver (unchanged code, confirmed via its SF/STENCIL_VAR macro expansion) always reads qL_prim_rsx_vf(${SF('')}$, ...) (unshifted) for the left state and qR_prim_rsx_vf(${SF(' + 1')}$, ...) (shifted +1) for the right state — i.e. L should be stored unshifted and R shifted +1, the exact opposite of what the new reconstruction routine does. The newly added s_igr_correct_lf_fluxes (same file) also reads with the solver's convention (u_L = qL_vf(${SF('')}$,...), u_R = qR_vf(${SF(' + 1')}$,...)), so it is likewise fed mismatched data by the reconstruction routine. Only the also-new s_igr_get_viscous happens to read with the same (backwards) offset the reconstruction wrote, which is internally consistent but disagrees with the solver/flux-correction routines that consume the very same qL_rsx_vf/qR_rsx_vf arrays.
  • Failure scenario: for any IGR run (which per case_validator.py's new riemann_solver != 5 check must now go through this shared LF path), each call to s_riemann_solver/s_lf_riemann_solver and s_igr_correct_lf_fluxes combines the left state reconstructed for face j-1 with the right state reconstructed for face j+1 instead of the matching pair for face j. This silently corrupts the flux/advection-source computation at every cell face for every IGR simulation (density, momentum, energy, and volume-fraction fluxes alike), producing wrong results without any crash or diagnostic — the class of bug this project's regression tests exist to catch, but a systemic index-convention swap here would shift results uniformly rather than in a way that is easy to spot from a corner case.

[StevenXue98]

Regarding the Claude code review comment, it was correct that the new subroutine s_igr_reconstruct_cell_boundary_values is storing in this fashion:

qL_rsx_vf(j + 1, k, l, i) = alpha_rho_L(i)
qR_rsx_vf(j, k, l, i) = alpha_rho_R(i)

This is the same convention used by s_weno and s_muscl; for example MUSCL writes vL_rs_vf_x(j) as the left cell boundary and vR_rs_vf_x(j) as the right cell boundary. I believe the reason this works is that s(_igr)_reconstruct_cell_boundary_values and s_lf_riemann_solver interpret L/R differently. In s(_igr)_reconstruct_cell_boundary_values, L/R face values are reconstructed relative to cell centers, but in s_lf_riemann_solver, L/R face values are read relative to cell interfaces. So the left cell face value of cell j + 1 is the right state at cell interface j + 1/2, and the right cell face value of cell j is the left state at cell interface j + 1/2. This is why we do

call s_riemann_solver(qR_rsx_vf, dqR_prim_dx_n(id)%vf, dqR_prim_dy_n(id)%vf, dqR_prim_dz_n(id)%vf, qR_prim(id)%vf, &
                                  & qL_rsx_vf, dqL_prim_dx_n(id)%vf, dqL_prim_dy_n(id)%vf, dqL_prim_dz_n(id)%vf, qL_prim(id)%vf, &
                                  & q_prim_qp%vf, flux_n(id)%vf, flux_src_n(id)%vf, flux_gsrc_n(id)%vf, id, irx, iry, irz)

in s_compute_rhs. As we store the face values at stencil j to slot j + 1 for the left arrays and stencil j to slot j for the right arrays, they are used in s_lf_riemann_solver as (shifted +1) for the right state and (unshifted) for the left state, which is consistent with the codebase's established convention.