Merge remote-tracking branch 'upstream/master' into MPIRefactor

Author: wilfonba

.cursor/rules/mfc-agent-rules.mdc

@@ -0,0 +1,82 @@
+---
+description: Full MFC project rules – consolidated for Agent Mode
+alwaysApply: true
+---
+
+# 0  Purpose & Scope
+Consolidated guidance for the MFC exascale, many-physics solver.  
+Written primarily for Fortran/Fypp; the OpenACC and style sections matter only when
+`.fpp` / `.f90` files are in view.
+
+---
+
+# 1  Global Project Context (always)
+- **Project**: *MFC* is modern Fortran 2008+ generated with **Fypp**.  
+  - Sources `src/`, tests `tests/`, examples `examples/`.  
+  - Most sources are `.fpp`; CMake transpiles them to `.f90`.  
+- **Fypp macros** live in `src/<subprogram>/include/` you should scan these first.  
+  `<subprogram>` ∈ {`simulation`,`common`,`pre_process`,`post_process`}.  
+- Only `simulation` (+ its `common` calls) is GPU-accelerated via **OpenACC**.  
+- Assume free-form Fortran 2008+, `implicit none`, explicit `intent`, and modern
+  intrinsics.  
+- Prefer `module … contains … subroutine foo()`; avoid `COMMON` blocks and
+  file-level `include` files.  
+- **Read the full codebase and docs *before* changing code.**  
+  Docs: <https://mflowcode.github.io/documentation/md_readme.html> and the respository root `README.md`.  
+
+### Incremental-change workflow
+1. Draft a step-by-step plan.  
+2. After each step, build:  
+   ```bash
+   ./mfc.sh build -t pre_process simulation -j $(nproc)
+    ```
+3. If it compiles, run focused tests:
+   ```bash
+   ./mfc.sh test -j $(nproc) -f EA8FA07E -t 9E2CA336
+   ```
+4. Roll back & fix if a step fails.
+
+* Do not run ./mfc.sh test -j $(nproc) without any other arguments (it takes too long to run all tests).
+
+---
+
+# 2  Style & Naming Conventions (for \*.fpp / \*.f90)
+
+* **Indent 2 spaces**; continuation lines align under `&`.
+* Lower-case keywords and intrinsics (`do`, `end subroutine`, …).
+* **Modules**: `m_<feature>` (e.g. `m_transport`).
+* **Public procedures**:
+  * Subroutine → `s_<verb>_<noun>` (e.g. `s_compute_flux`)
+  * Function   → `f_<verb>_<noun>`
+* Private helpers stay in the module; avoid nested procedures.
+* **Size limits**: subroutine ≤ 500 lines, helper ≤ 150, function ≤ 100,
+  module/file ≤ 1000.
+* ≤ 6 arguments per routine; otherwise pass a derived-type “params” struct.
+* No `goto` (except unavoidable legacy); no global state (`COMMON`, `save`).
+* Every variable: `intent(in|out|inout)` + appropriate `dimension` / `allocatable`
+  / `pointer`.
+* Use `s_mpi_abort(<msg>)` for errors, not `stop`.
+* Mark OpenACC-callable helpers that are called from OpenACC parallel loops immediately after declaration:
+  ```fortran
+  subroutine s_flux_update(...)
+    !$acc routine seq
+    ...
+  end subroutine
+  ```
+
+---
+
+# 3  OpenACC Programming Guidelines (for kernels)
+
+Wrap tight loops with
+
+```fortran
+!$acc parallel loop gang vector default(present) reduction(...)
+```
+* Add `collapse(n)` to merge nested loops when safe.
+* Declare loop-local variables with `private(...)`.
+* Allocate large arrays with `managed` or move them into a persistent
+  `!$acc enter data` region at start-up.
+* **Do not** place `stop` / `error stop` inside device code.
+* Must compile with Cray `ftn` and NVIDIA `nvfortran` for GPU offloading; also build CPU-only with
+  GNU `gfortran` and Intel `ifx`/`ifort`.

.github/workflows/phoenix/bench.sh

@@ -8,8 +8,13 @@ if [ "$job_device" == "gpu" ]; then
     device_opts="--gpu -g $gpu_ids"
 fi
 
+mkdir -p /storage/scratch1/6/sbryngelson3/mytmp_build
+export TMPDIR=/storage/scratch1/6/sbryngelson3/mytmp_build
+
 if ["$job_device" == "gpu"]; then
     ./mfc.sh bench --mem 12 -j $(nproc) -o "$job_slug.yaml" -- -c phoenix-bench $device_opts -n $n_ranks
 else
     ./mfc.sh bench --mem 1 -j $(nproc) -o "$job_slug.yaml" -- -c phoenix-bench $device_opts -n $n_ranks
 fi
+
+unset TMPDIR

CMakeLists.txt

@@ -135,13 +135,17 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
     if (CMAKE_BUILD_TYPE STREQUAL "Debug")
         add_compile_options(
             -Wall
+            -Wextra	
             -fcheck=all,no-array-temps
             -fbacktrace
             -fimplicit-none
-            #-ffpe-trap=invalid,zero,denormal,overflow
             -fsignaling-nans
             -finit-real=snan
             -finit-integer=-99999999
+            -Wintrinsic-shadow	
+            -Wunderflow	
+            -Wrealloc-lhs	
+            -Wsurprising	
 	    )
     endif()
 

misc/length-subroutines.sh

@@ -0,0 +1,29 @@
+#!/bin/bash
+
+# Use gawk if available, otherwise fall back to awk
+if command -v gawk > /dev/null; then
+    AWK_CMD="gawk"
+    IGNORECASE_BLOCK='BEGIN { IGNORECASE = 1 }'
+else
+    AWK_CMD="awk"
+    IGNORECASE_BLOCK=''
+    echo "Warning: gawk not found. Case-insensitive matching may not work as expected." >&2
+fi
+
+find . -type f \( -name "*.f90" -o -name "*.fpp" \) | while read file; do
+    "$AWK_CMD" "
+    $IGNORECASE_BLOCK
+    /^[ \t]*((pure|elemental|impure)[ \t]+)*subroutine[ \t]+[a-zA-Z_][a-zA-Z0-9_]*[ \t]*\\(/ {
+        in_sub = 1
+        match(\$0, /subroutine[ \t]+([a-zA-Z_][a-zA-Z0-9_]*)/, arr)
+        sub_name = arr[1]
+        start_line = NR
+        next
+    }
+    /^[ \t]*end[ \t]+subroutine[ \t]*([a-zA-Z_][a-zA-Z0-9_]*)?[ \t]*\$/ && in_sub {
+        end_line = NR
+        print (end_line - start_line + 1) \"\t\" FILENAME \": \" sub_name
+        in_sub = 0
+    }
+    " "$file"
+done | sort -nr | awk -F'\t' '{print $2 " : " $1 " lines"}' | head -20

src/common/m_boundary_common.fpp

@@ -156,7 +156,7 @@ contains
                     case (BC_CHAR_SUP_OUTFLOW:BC_GHOST_EXTRAP)
                         call s_ghost_cell_extrapolation(q_prim_vf, 2, -1, k, l)
                     case (BC_AXIS)
-                        call s_axis(q_prim_vf, 2, -1, k, l, pb, mv)
+                        call s_axis(q_prim_vf, pb, mv, k, l)
                     case (BC_REFLECTIVE)
                         call s_symmetry(q_prim_vf, 2, -1, k, l, pb, mv)
                     case (BC_PERIODIC)
@@ -735,11 +735,10 @@ contains
 
     end subroutine s_periodic
 
-    pure subroutine s_axis(q_prim_vf, bc_dir, bc_loc, k, l, pb, mv)
-        !$acc routine seq
+    pure subroutine s_axis(q_prim_vf, pb, mv, k, l)
+!$acc routine seq
         type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
-        real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv
-        integer, intent(in) :: bc_dir, bc_loc
+        real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv
         integer, intent(in) :: k, l
 
         integer :: j, q, i

src/common/m_finite_differences.fpp

@@ -15,8 +15,7 @@ contains
         integer :: x, y, z !< Generic loop iterators
 
         real(wp) :: divergence
-
-        !$acc parallel loop collapse(3) private(divergence)
+        !$acc parallel loop collapse(3) gang vector default(present) private(divergence)
         do x = ix_s%beg, ix_s%end
             do y = iy_s%beg, iy_s%end
                 do z = iz_s%beg, iz_s%end

src/common/m_mpi_common.fpp

@@ -350,21 +350,17 @@ contains
         !!  @param Rc_min_glb Global minimum Rc stability criterion
     impure subroutine s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, &
                                                               vcfl_max_loc, &
-                                                              ccfl_max_loc, &
                                                               Rc_min_loc, &
                                                               icfl_max_glb, &
                                                               vcfl_max_glb, &
-                                                              ccfl_max_glb, &
                                                               Rc_min_glb)
 
         real(wp), intent(in) :: icfl_max_loc
         real(wp), intent(in) :: vcfl_max_loc
-        real(wp), intent(in) :: ccfl_max_loc
         real(wp), intent(in) :: Rc_min_loc
 
         real(wp), intent(out) :: icfl_max_glb
         real(wp), intent(out) :: vcfl_max_glb
-        real(wp), intent(out) :: ccfl_max_glb
         real(wp), intent(out) :: Rc_min_glb
 
 #ifdef MFC_SIMULATION

src/common/m_phase_change.fpp

@@ -143,7 +143,7 @@ contains
 
                     ! Calling pT-equilibrium for either finishing phase-change module, or as an IC for the pTg-equilibrium
                     ! for this case, MFL cannot be either 0 or 1, so I chose it to be 2
-                    call s_infinite_pt_relaxation_k(j, k, l, 2, pS, p_infpT, rM, q_cons_vf, rhoe, TS)
+                    call s_infinite_pt_relaxation_k(j, k, l, 2, pS, p_infpT, q_cons_vf, rhoe, TS)
 
                     ! check if pTg-equilibrium is required
                     ! NOTE that NOTHING else needs to be updated OTHER than the individual partial densities
@@ -164,7 +164,7 @@ contains
                         q_cons_vf(vp + contxb - 1)%sf(j, k, l) = (1.0_wp - mixM)*rM
 
                         ! calling pT-equilibrium for overheated vapor, which is MFL = 0
-                        call s_infinite_pt_relaxation_k(j, k, l, 0, pSOV, p_infOV, rM, q_cons_vf, rhoe, TSOV)
+                        call s_infinite_pt_relaxation_k(j, k, l, 0, pSOV, p_infOV, q_cons_vf, rhoe, TSOV)
 
                         ! calculating Saturation temperature
                         call s_TSat(pSOV, TSatOV, TSOV)
@@ -177,7 +177,7 @@ contains
                         q_cons_vf(vp + contxb - 1)%sf(j, k, l) = mixM*rM
 
                         ! calling pT-equilibrium for subcooled liquid, which is MFL = 1
-                        call s_infinite_pt_relaxation_k(j, k, l, 1, pSSL, p_infSL, rM, q_cons_vf, rhoe, TSSL)
+                        call s_infinite_pt_relaxation_k(j, k, l, 1, pSSL, p_infSL, q_cons_vf, rhoe, TSSL)
 
                         ! calculating Saturation temperature
                         call s_TSat(pSSL, TSatSL, TSSL)
@@ -281,7 +281,7 @@ contains
         !!  @param q_cons_vf Cell-average conservative variables
         !!  @param rhoe mixture energy
         !!  @param TS equilibrium temperature at the interface
-    pure subroutine s_infinite_pt_relaxation_k(j, k, l, MFL, pS, p_infpT, rM, q_cons_vf, rhoe, TS)
+    pure subroutine s_infinite_pt_relaxation_k(j, k, l, MFL, pS, p_infpT, q_cons_vf, rhoe, TS)
 
 #ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_infinite_pt_relaxation_k
@@ -293,7 +293,6 @@ contains
         integer, intent(in) :: j, k, l, MFL
         real(wp), intent(out) :: pS
         real(wp), dimension(num_fluids), intent(out) :: p_infpT
-        real(wp), intent(in) :: rM
         type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
         real(wp), intent(in) :: rhoe
         real(wp), intent(out) :: TS

src/common/m_variables_conversion.fpp

@@ -91,11 +91,11 @@ contains
 
         if (model_eqns == 1) then        ! Gamma/pi_inf model
             call s_convert_mixture_to_mixture_variables(q_vf, i, j, k, &
-                                                        rho, gamma, pi_inf, qv, Re_K, G_K, G)
+                                                        rho, gamma, pi_inf, qv)
 
         else if (bubbles_euler) then
             call s_convert_species_to_mixture_variables_bubbles(q_vf, i, j, k, &
-                                                                rho, gamma, pi_inf, qv, Re_K, G_K, G)
+                                                                rho, gamma, pi_inf, qv, Re_K)
         else
             ! Volume fraction model
             call s_convert_species_to_mixture_variables(q_vf, i, j, k, &
@@ -206,7 +206,7 @@ contains
         !! @param pi_inf liquid stiffness
         !! @param qv fluid reference energy
     subroutine s_convert_mixture_to_mixture_variables(q_vf, i, j, k, &
-                                                      rho, gamma, pi_inf, qv, Re_K, G_K, G)
+                                                      rho, gamma, pi_inf, qv)
 
         type(scalar_field), dimension(sys_size), intent(in) :: q_vf
         integer, intent(in) :: i, j, k
@@ -216,11 +216,6 @@ contains
         real(wp), intent(out), target :: pi_inf
         real(wp), intent(out), target :: qv
 
-        real(wp), optional, dimension(2), intent(out) :: Re_K
-
-        real(wp), optional, intent(out) :: G_K
-        real(wp), optional, dimension(num_fluids), intent(in) :: G
-
         ! Transferring the density, the specific heat ratio function and the
         ! liquid stiffness function, respectively
         rho = q_vf(1)%sf(i, j, k)
@@ -253,7 +248,7 @@ contains
         !! @param pi_inf liquid stiffness
         !! @param qv fluid reference energy
     subroutine s_convert_species_to_mixture_variables_bubbles(q_vf, j, k, l, &
-                                                              rho, gamma, pi_inf, qv, Re_K, G_K, G)
+                                                              rho, gamma, pi_inf, qv, Re_K)
 
         type(scalar_field), dimension(sys_size), intent(in) :: q_vf
 
@@ -265,8 +260,6 @@ contains
         real(wp), intent(out), target :: qv
 
         real(wp), optional, dimension(2), intent(out) :: Re_K
-        real(wp), optional, intent(out) :: G_K
-        real(wp), optional, dimension(num_fluids), intent(in) :: G
 
         integer :: i, q
         real(wp), dimension(num_fluids) :: alpha_rho_K, alpha_K
@@ -463,7 +456,7 @@ contains
 
     pure subroutine s_convert_species_to_mixture_variables_acc(rho_K, &
                                                                gamma_K, pi_inf_K, qv_K, &
-                                                               alpha_K, alpha_rho_K, Re_K, k, l, r, &
+                                                               alpha_K, alpha_rho_K, Re_K, &
                                                                G_K, G)
 #ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_convert_species_to_mixture_variables_acc
@@ -480,8 +473,6 @@ contains
         real(wp), optional, intent(out) :: G_K
         real(wp), optional, dimension(num_fluids), intent(in) :: G
 
-        integer, intent(in) :: k, l, r
-
         integer :: i, j !< Generic loop iterators
         real(wp) :: alpha_K_sum
 
@@ -547,7 +538,7 @@ contains
 
     pure subroutine s_convert_species_to_mixture_variables_bubbles_acc(rho_K, &
                                                                        gamma_K, pi_inf_K, qv_K, &
-                                                                       alpha_K, alpha_rho_K, Re_K, k, l, r)
+                                                                       alpha_K, alpha_rho_K, Re_K)
 #ifdef _CRAYFTN
         !DIR$ INLINEALWAYS s_convert_species_to_mixture_variables_bubbles_acc
 #else
@@ -560,7 +551,6 @@ contains
             !! Partial densities and volume fractions
 
         real(wp), dimension(2), intent(out) :: Re_K
-        integer, intent(in) :: k, l, r
 
         integer :: i, j !< Generic loop iterators
 
@@ -819,16 +809,12 @@ contains
     subroutine s_convert_conservative_to_primitive_variables(qK_cons_vf, &
                                                              q_T_sf, &
                                                              qK_prim_vf, &
-                                                             ibounds, &
-                                                             gm_alphaK_vf)
+                                                             ibounds)
 
         type(scalar_field), dimension(sys_size), intent(in) :: qK_cons_vf
         type(scalar_field), intent(inout) :: q_T_sf
         type(scalar_field), dimension(sys_size), intent(inout) :: qK_prim_vf
         type(int_bounds_info), dimension(1:3), intent(in) :: ibounds
-        type(scalar_field), &
-            allocatable, optional, dimension(:), &
-            intent(in) :: gm_alphaK_vf
 
         real(wp), dimension(num_fluids) :: alpha_K, alpha_rho_K
         real(wp), dimension(2) :: Re_K
@@ -902,13 +888,13 @@ contains
                         ! If in simulation, use acc mixture subroutines
                         if (elasticity) then
                             call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, &
-                                                                            alpha_rho_K, Re_K, j, k, l, G_K, Gs)
+                                                                            alpha_rho_K, Re_K, G_K, Gs)
                         else if (bubbles_euler) then
                             call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, pi_inf_K, qv_K, &
-                                                                                    alpha_K, alpha_rho_K, Re_K, j, k, l)
+                                                                                    alpha_K, alpha_rho_K, Re_K)
                         else
                             call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, &
-                                                                            alpha_K, alpha_rho_K, Re_K, j, k, l)
+                                                                            alpha_K, alpha_rho_K, Re_K)
                         end if
 #else
                         ! If pre-processing, use non acc mixture subroutines
@@ -1513,13 +1499,13 @@ contains
                     if (elasticity) then
                         call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, &
                                                                         alpha_K, alpha_rho_K, Re_K, &
-                                                                        j, k, l, G_K, Gs)
+                                                                        G_K, Gs)
                     else if (bubbles_euler) then
                         call s_convert_species_to_mixture_variables_bubbles_acc(rho_K, gamma_K, &
-                                                                                pi_inf_K, qv_K, alpha_K, alpha_rho_K, Re_K, j, k, l)
+                                                                                pi_inf_K, qv_K, alpha_K, alpha_rho_K, Re_K)
                     else
                         call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, &
-                                                                        alpha_K, alpha_rho_K, Re_K, j, k, l)
+                                                                        alpha_K, alpha_rho_K, Re_K)
                     end if
 
                     ! Computing the energy from the pressure