sbryngelson · Oct 3, 2026 · Oct 3, 2026 · Oct 3, 2026 · Oct 3, 2026 · Oct 5, 2026
@@ -20,6 +20,7 @@ project(MFC LANGUAGES C CXX Fortran)
 option(MFC_MPI           "Build with MPI"                                     ON)
 option(MFC_OpenACC       "Build with OpenACC"                                OFF)
 option(MFC_OpenMP        "Build with OpenMP"                                 OFF)
+option(MFC_FIXED_BOUNDS  "Compile-time bounds for per-thread GPU arrays"     ON)
 option(MFC_GCov          "Build with GCov"                                   OFF)
 option(MFC_Unified       "Build with unified CPU & GPU memory (GH-200 only)" OFF)
 option(MFC_Fastmath      "Build with -gpu=fastmath on NV GPUs"               OFF)

@@ -101,6 +101,12 @@ macro(HANDLE_SOURCES target useCommon)
         list(APPEND ${target}_incs ${common_incs})
     endif()
 
+    # Fixed per-thread array bounds (see shared_parallel_macros.fpp): only simulation is offloaded.
+    set(_fixed_bounds False)
+    if (MFC_FIXED_BOUNDS AND (MFC_OpenACC OR MFC_OpenMP) AND "${target}" STREQUAL "simulation")
+        set(_fixed_bounds True)
+    endif()
+
     # /path/to/*.fpp (used by <target>) -> <build>/fypp/<target>/*.f90
     file(MAKE_DIRECTORY "${CMAKE_BINARY_DIR}/fypp/${target}")
     foreach(fpp ${${target}_FPPs})
@@ -118,6 +124,7 @@ macro(HANDLE_SOURCES target useCommon)
                                  -D MFC_${CMAKE_Fortran_COMPILER_ID}
                                  -D MFC_${${target}_UPPER}
                                  -D MFC_COMPILER="${CMAKE_Fortran_COMPILER_ID}"
+                                 -D MFC_FIXED_BOUNDS=${_fixed_bounds}
 				                 -D MFC_CASE_OPTIMIZATION=False
                                  -D chemistry=False
                                  --line-numbering

@@ -9,7 +9,7 @@ file(GLOB_RECURSE _mfc_gen_inputs
     "${CMAKE_CURRENT_SOURCE_DIR}/toolchain/mfc/params/*.py"
 )
 
-# Enumerate the 18 generated .fpp files explicitly so ninja can track them as
+# Enumerate the 21 generated .fpp files explicitly so ninja can track them as
 # build-time outputs and so HANDLE_SOURCES does not need a configure-time GLOB
 # of ${CMAKE_BINARY_DIR}/include/<target>/ (which fails when the dir is empty).
 set(_mfc_gen_inc "${CMAKE_BINARY_DIR}/include")
@@ -18,6 +18,7 @@ set(_mfc_gen_files_pre_process
     "${_mfc_gen_inc}/pre_process/generated_decls.fpp"
     "${_mfc_gen_inc}/pre_process/generated_constants.fpp"
     "${_mfc_gen_inc}/pre_process/generated_eos.fpp"
+    "${_mfc_gen_inc}/pre_process/generated_eqn_idx.fpp"
     "${_mfc_gen_inc}/pre_process/generated_bcast.fpp"
     "${_mfc_gen_inc}/pre_process/generated_case_opt_decls.fpp"
 )
@@ -26,6 +27,7 @@ set(_mfc_gen_files_simulation
     "${_mfc_gen_inc}/simulation/generated_decls.fpp"
     "${_mfc_gen_inc}/simulation/generated_constants.fpp"
     "${_mfc_gen_inc}/simulation/generated_eos.fpp"
+    "${_mfc_gen_inc}/simulation/generated_eqn_idx.fpp"
     "${_mfc_gen_inc}/simulation/generated_bcast.fpp"
     "${_mfc_gen_inc}/simulation/generated_case_opt_decls.fpp"
 )
@@ -34,6 +36,7 @@ set(_mfc_gen_files_post_process
     "${_mfc_gen_inc}/post_process/generated_decls.fpp"
     "${_mfc_gen_inc}/post_process/generated_constants.fpp"
     "${_mfc_gen_inc}/post_process/generated_eos.fpp"
+    "${_mfc_gen_inc}/post_process/generated_eqn_idx.fpp"
     "${_mfc_gen_inc}/post_process/generated_bcast.fpp"
     "${_mfc_gen_inc}/post_process/generated_case_opt_decls.fpp"
 )

@@ -87,7 +87,7 @@ between them is narrow.
 mfc.sh (env bootstrap, venv, module loading, lock)
   └─ toolchain/mfc/build.py (config slugs, cmake invocation)
        └─ CMakeLists.txt + cmake/{GPU,Fypp,ParamsCodegen,MFCTargets}.cmake
-            └─ toolchain/mfc/params/generators/cmake_gen.py (writes 18 generated .fpp includes)
+            └─ toolchain/mfc/params/generators/cmake_gen.py (writes 21 generated .fpp includes)
 ```
 
 **`mfc.sh` → `build.py`.**  `mfc.sh` is a thin shell wrapper that activates the Python
@@ -100,8 +100,8 @@ mode, debug, chemistry, MPI.  Staging and install trees are namespaced by slug u
 **CMake layer.**  `cmake/Fypp.cmake` defines `HANDLE_SOURCES`, which sets up one
 `add_custom_command` per `.fpp` file to run Fypp at build time.  `cmake/ParamsCodegen.cmake`
 registers a single ninja-tracked `add_custom_command` (DEPENDS all `params/*.py`) that
-invokes `cmake_gen.py` and writes the 18 generated includes under
-`build/include/<target>/`.  There is no configure-time generation: all 18 files are build
+invokes `cmake_gen.py` and writes the 21 generated includes under
+`build/include/<target>/`.  There is no configure-time generation: all 21 files are build
 outputs, so changing any `params/*.py` triggers only a targeted rebuild, not a full
 reconfigure.
 

@@ -650,14 +650,12 @@ helper receives only scalars or small arrays with explicit-shape dimensioning. T
 is required because the `SF` indexing lambda used in solver loops is defined locally
 inside each solver's `#:for NORM_DIR` block and cannot be referenced from a helper.
 
-**AMD case-opt compatibility.** Under `--case-optimization` with the AMD backend,
-arrays that are sized by runtime parameters at compile time must be declared with an
-explicit constant bound. Use an explicit `n` argument (e.g., `integer, intent(in) ::
-nf`) and dimension helpers as `dimension(nf)` rather than `dimension(num_fluids)`.
-See `s_compute_interface_reynolds` in `src/simulation/m_riemann_state.fpp` for the
-`#:if not MFC_CASE_OPTIMIZATION and USING_AMD` guard pattern: the guard sits on the
-dummy-argument declaration in the helper's definition, with matching guards on the
-callers' own local declarations so the actual and dummy bounds agree.
+**Per-thread array extents.** Size per-thread arrays with ``${BOUND('name')}$``
+(`src/common/include/shared_parallel_macros.fpp`), e.g. ``dimension(${BOUND('num_fluids')}$)``,
+not `dimension(num_fluids)`. Use it on both a helper's dummy arguments and its callers' locals
+so the bounds agree. `BOUND` gives the exact value under `--case-optimization`, a fixed maximum
+in GPU simulation builds (`MFC_FIXED_BOUNDS`), and the runtime extent otherwise. A new extent
+needs an entry in `BOUND_MAX` and, if it can exceed it, a check in `s_check_fixed_bounds`.
 
 **Declare scoping.** The `GPU_ROUTINE` directive must appear in the source file
 that defines the routine. Helpers added to `m_riemann_state.fpp` are automatically
@@ -943,13 +941,11 @@ answer is wrong, or one backend diverges from all the others.
   passes a `parameter` array from `m_thermochem`, such as `molecular_weights`, into a
   declare-target routine. Read such arrays directly in the kernel, or pass a plain local
   computed from them.
-- **The `USING_AMD` fypp guards are load-bearing for performance.** They swap a device-global
-  array bound for a literal in `src/common/include/shared_parallel_macros.fpp` and its use
-  sites. On AFAR 24.3, `USING_AMD = False` (no case optimization) gives correct results but runs
-  4-5x slower: private arrays sized by runtime globals move from registers to scratch (the
-  WENO kernel goes from 0 to 400 B of scratch per thread and runs 11x slower; HLLC 4.9x). On
-  AFAR 23.2.x it also produced NaNs in CBC, `wave_speeds=2`, IBM, surface tension, QBMM,
-  viscous, and MHD HLLD cases. Removing them needs a benchmark, not just the tests.
+- **Runtime-sized per-thread arrays are slow on GPUs.** A private array sized by a runtime value
+  (`num_fluids`, `sys_size`) cannot live in registers and spills to scratch: on amdflang (AFAR
+  24.3) the WENO kernel ran 11x slower and HLLC 4.9x. This is why `BOUND` gives fixed maxima in
+  GPU builds; turning them off (`-DMFC_FIXED_BOUNDS=OFF`) needs a benchmark, not just the tests.
+  On AFAR 23.2.x the runtime-sized form also gave NaNs, so passing tests alone is not enough.
 - `@:ACC_SETUP_VFs` and `@:ACC_SETUP_SFs` compile only under Cray. Around MPI, use
   `GPU_UPDATE(host=...)` before a send and `GPU_UPDATE(device=...)` after a receive.
 

@@ -11,9 +11,25 @@
 #! NUM_SPECIES (m_thermochem's species count) and CHEMISTRY, written per build by the toolchain.
 #:include 'thermochem.fpp'
 
-#! Fallback extent the USING_AMD guards substitute for sys_size arrays when case optimization is off.
+#! Compile-time extents for per-thread arrays, since runtime-sized private arrays spill to scratch
+#! on GPUs. Case optimization gives exact ones (parameters, plus CASE_OPT_SIZES from the layout
+#! table); otherwise GPU simulation builds (MFC_FIXED_BOUNDS, set per target by CMake) use maxima.
 #! Chemistry pins num_fluids to 1, leaving at most 10 flow variables beside the species.
-#:set AMD_SYS_SIZE_MAX = 10 + NUM_SPECIES if CHEMISTRY else 70
+#:set MFC_FIXED_BOUNDS = defined('MFC_FIXED_BOUNDS') and MFC_FIXED_BOUNDS
+#:set FIXED_BOUNDS = MFC_FIXED_BOUNDS and not MFC_CASE_OPTIMIZATION
+#:set NUM_FLUIDS_MAX = 3
+#:set NB_MAX = 3
+#:set SYS_SIZE_MAX = 10 + NUM_SPECIES if CHEMISTRY else 70
+#:set BOUND_MAX = {'num_fluids': NUM_FLUIDS_MAX, 'nb': NB_MAX, 'num_dims': 3, 'num_vels': 3, &
+    & 'weno_polyn': 3, 'weno_num_stencils': 4, 'nterms': 32, 'n_stress': 6, 'sys_size': SYS_SIZE_MAX}
+#:set BOUND_RUNTIME = {'n_stress': 'eqn_idx%stress%end - eqn_idx%stress%beg + 1'}
+
+#! Extent for a per-thread array sized by `name`: exact under case optimization, else the fixed
+#! maximum in GPU simulation builds, else the runtime expression.
+#:def BOUND(name)
+    $:CASE_OPT_SIZES.get(name, &
+                         & name) if MFC_CASE_OPTIMIZATION else BOUND_MAX[name] if MFC_FIXED_BOUNDS else BOUND_RUNTIME.get(name, name)
+#:enddef
 
 #:def ASSERT_LIST(data, datatype)
     #:assert data is not None

@@ -26,8 +26,8 @@ contains
         integer(kind=8), intent(in) :: n_global
 
         if (check_total_cells) call s_check_total_cells(n_global)
-        #:if USING_AMD
-            call s_check_amd
+        #:if MFC_FIXED_BOUNDS
+            call s_check_fixed_bounds
         #:endif
 
     end subroutine s_check_inputs_common
@@ -48,17 +48,16 @@ contains
 
     end subroutine s_check_total_cells
 
-    !> Check that simulation parameters stay within AMD GPU compiler limits when case optimization is disabled.
-    impure subroutine s_check_amd
+    !> Check that the case fits the fixed per-thread array bounds of a GPU simulation build.
+    impure subroutine s_check_fixed_bounds
 
         #:if not MFC_CASE_OPTIMIZATION
-            @:PROHIBIT(num_fluids > 3, "num_fluids <= 3 for AMDFLang when Case optimization is off")
-            @:PROHIBIT((bubbles_euler .or. bubbles_lagrange) .and. nb > 3, "nb <= 3 for AMDFLang when Case optimization is off")
-            ! HLLC's star states and the CBC L vectors have no other bound.
-            @:PROHIBIT(sys_size > ${AMD_SYS_SIZE_MAX}$, &
-                       & "sys_size <= ${AMD_SYS_SIZE_MAX}$ for AMDFLang when Case optimization is off")
+            @:PROHIBIT(num_fluids > ${NUM_FLUIDS_MAX}$, &
+                       & "num_fluids <= ${NUM_FLUIDS_MAX}$ in GPU builds; rebuild with --case-optimization")
+            @:PROHIBIT((bubbles_euler .or. bubbles_lagrange) .and. nb > ${NB_MAX}$, &
+                       & "nb <= ${NB_MAX}$ in GPU builds; rebuild with --case-optimization")
         #:endif
 
-    end subroutine s_check_amd
+    end subroutine s_check_fixed_bounds
 
 end module m_checker_common
@@ -305,15 +305,11 @@ contains
 
         $:GPU_ROUTINE(function_name='f_mixture_temperature', parallelism='[seq]', cray_inline=True)
 
-        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
-            real(wp), dimension(3), intent(in) :: alpha_rho_K
-        #:else
-            real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K
-        #:endif
-        real(wp), intent(in) :: pres, gamma_K, pi_inf_K
-        real(wp)             :: T
-        real(wp)             :: mCP  !< sum of alpha_rho_i*cp_i; cp_i = n_i*cv_i
-        integer              :: i
+        real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: alpha_rho_K
+        real(wp), intent(in)                                     :: pres, gamma_K, pi_inf_K
+        real(wp)                                                 :: T
+        real(wp)                                                 :: mCP  !< sum of alpha_rho_i*cp_i; cp_i = n_i*cv_i
+        integer                                                  :: i
 
         mCP = 0._wp
         $:GPU_LOOP(parallelism='[seq]')
@@ -569,15 +565,11 @@ contains
 
         $:GPU_ROUTINE(function_name='s_compute_mixture_coefficients', parallelism='[seq]', cray_inline=True)
 
-        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
-            real(wp), dimension(3), intent(in) :: alpha_rho_K, alpha_K
-        #:else
-            real(wp), dimension(num_fluids), intent(in) :: alpha_rho_K, alpha_K
-        #:endif
-        real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K
-        real(wp)              :: gamma_i, pi_inf_i, dpi_i, dgamma_i
-        real(wp)              :: rho_i, alpha_i, alpha_rho_i
-        integer               :: i  !< Loop iterator over fluids
+        real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: alpha_rho_K, alpha_K
+        real(wp), intent(out)                                    :: rho_K, gamma_K, pi_inf_K, qv_K
+        real(wp)                                                 :: gamma_i, pi_inf_i, dpi_i, dgamma_i
+        real(wp)                                                 :: rho_i, alpha_i, alpha_rho_i
+        integer                                                  :: i  !< Loop iterator over fluids
 
         ! The bubbly closure is written for one carrier liquid, which keeps its own coefficients
         ! undiluted: Gamma_l*p_l = (E - rho|u|^2/2)/(1 - alf) - Pi_inf_l, the void entering only through
@@ -614,14 +606,10 @@ contains
 
         $:GPU_ROUTINE(function_name='s_compute_mixture_coefficients_dt', parallelism='[seq]', cray_inline=True)
 
-        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
-            real(wp), dimension(3), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv
-        #:else
-            real(wp), dimension(num_fluids), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv
-        #:endif
-        real(wp), intent(out) :: drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt
-        real(wp)              :: rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i, alpha_i, alpha_rho_i
-        integer               :: i  !< Loop iterator over fluids
+        real(wp), dimension(${BOUND('num_fluids')}$), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv
+        real(wp), intent(out)                                    :: drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt
+        real(wp)                                                 :: rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i, alpha_i, alpha_rho_i
+        integer                                                  :: i  !< Loop iterator over fluids
 
         dgamma_dt = 0._wp
         dpi_inf_dt = 0._wp
@@ -654,21 +642,13 @@ contains
 
         $:GPU_ROUTINE(parallelism='[seq]')
 
-        real(wp), intent(in) :: pres, rho, gamma, pi_inf
-        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
-            real(wp), dimension(3), intent(in) :: adv
-        #:else
-            real(wp), dimension(num_fluids), intent(in) :: adv
-        #:endif
-        real(wp), intent(out) :: c
-        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
-            real(wp), dimension(3), intent(in), optional :: alpha_rho
-        #:else
-            real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho
-        #:endif
-        real(wp) :: alf  !< Subgrid void fraction; dilute by construction
-        real(wp) :: blkmod_q, alpha_q, alpha_rho_q, gamma_q, pi_inf_q
-        integer  :: q
+        real(wp), intent(in)                                               :: pres, rho, gamma, pi_inf
+        real(wp), dimension(${BOUND('num_fluids')}$), intent(in)           :: adv
+        real(wp), intent(out)                                              :: c
+        real(wp), dimension(${BOUND('num_fluids')}$), intent(in), optional :: alpha_rho
+        real(wp)                                                           :: alf  !< Subgrid void fraction; dilute by construction
+        real(wp)                                                           :: blkmod_q, alpha_q, alpha_rho_q, gamma_q, pi_inf_q
+        integer                                                            :: q
 
         if (chemistry) then  ! Reacting mixture sound speed
             c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho)
@@ -741,18 +721,10 @@ contains
 
         $:GPU_ROUTINE(parallelism='[seq]')
 
-        real(wp), intent(in) :: pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c
-        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
-            real(wp), dimension(3), intent(in) :: adv
-        #:else
-            real(wp), dimension(num_fluids), intent(in) :: adv
-        #:endif
-        real(wp), intent(out) :: c
-        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
-            real(wp), dimension(3), intent(in), optional :: alpha_rho
-        #:else
-            real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho
-        #:endif
+        real(wp), intent(in)                                               :: pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c
+        real(wp), dimension(${BOUND('num_fluids')}$), intent(in)           :: adv
+        real(wp), intent(out)                                              :: c
+        real(wp), dimension(${BOUND('num_fluids')}$), intent(in), optional :: alpha_rho
 
         if (chemistry) then  ! Reacting mixture sound speed
             if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then