Work around CCE 19.0.0 compiler bugs for Cray+OpenACC builds

Three distinct CCE 19.0.0 compiler bugs required fixes:

Bug 1: InstCombine ICE in matmul() in m_phase_change.fpp
- Replace matmul() with explicit 2x2 arithmetic

Bug 2: IPA bring_routine_resident SIGSEGV in m_phase_change.fpp
- Add -Oipa0 per-file in CMakeLists.txt (Cray+OpenACC only)
- Use cray_noinline=True on 4 GPU_ROUTINE calls in m_phase_change.fpp
  and 4 in m_variables_conversion.fpp

Bug 3: IPA castIsValid ICE in m_bubbles_EL.fpp
- Change proc_bubble_counts from VLA to allocatable
- Add -Oipa0 per-file in CMakeLists.txt (Cray+OpenACC only)

Bug 4: m_chemistry.fpp VLA ICE in case-optimized pre_process builds
- Guard 4 dimension(num_species) local arrays with USING_CCE

Bug 5: Pyrometheus GPU_ROUTINE macro missing !acc routine seq on Cray+ACC
- Post-process generated m_thermochem.f90 in toolchain/mfc/run/input.py
  to replace the broken Cray INLINEALWAYS-only macro with plain
  #define GPU_ROUTINE(name) !acc routine seq

Also fix uninitialized FT in s_TSat (use huge(1.0_wp) not huge(FT)).

See PR #1286.

Author: Spencer Bryngelson

CMakeLists.txt

@@ -397,6 +397,7 @@ HANDLE_SOURCES(simulation ON)
 HANDLE_SOURCES(post_process ON)
 HANDLE_SOURCES(syscheck OFF)
 
+
 # MFC_SETUP_TARGET: Given a target (herein <target>), this macro creates a new
 # executable <target> with the appropriate sources, compiler definitions, and
 # linked libraries (assuming HANDLE_SOURCES was called on <target>).
@@ -633,6 +634,19 @@ if (MFC_SIMULATION)
     MFC_SETUP_TARGET(TARGET  simulation
                      SOURCES "${simulation_SRCs}"
                      MPI FFTW OpenACC OpenMP)
+    # CCE 19.0.0 IPA workaround: two files trigger IPA crashes:
+    #   m_bubbles_EL: castIsValid assertion (InstCombine/foldIntegerTypedPHI)
+    #   m_phase_change: bring_routine_resident SIGSEGV
+    # Disabling IPA per-file avoids the crashes while preserving IPA for
+    # the rest of simulation (needed for thermochem INLINEALWAYS inlining).
+    # See PR #1286.
+    if (CMAKE_Fortran_COMPILER_ID STREQUAL "Cray" AND MFC_OpenACC)
+        set_source_files_properties(
+            "${CMAKE_BINARY_DIR}/fypp/simulation/m_bubbles_EL.fpp.f90"
+            "${CMAKE_BINARY_DIR}/fypp/simulation/m_phase_change.fpp.f90"
+            PROPERTIES COMPILE_FLAGS "-Oipa0"
+        )
+    endif()
 endif()
 
 if (MFC_POST_PROCESS)

src/common/include/parallel_macros.fpp

@@ -48,18 +48,44 @@
 
 #:enddef
 
-#:def GPU_ROUTINE(function_name=None, parallelism=None, nohost=False, cray_inline=False, extraAccArgs=None, extraOmpArgs=None)
+#:def GPU_ROUTINE(function_name=None, parallelism=None, nohost=False, cray_inline=False, cray_noinline=False, extraAccArgs=None, extraOmpArgs=None)
     #:assert isinstance(cray_inline, bool)
+    #:assert isinstance(cray_noinline, bool)
+    #:assert not (cray_inline and cray_noinline), "cray_inline and cray_noinline are mutually exclusive"
     #:set acc_directive = ACC_ROUTINE(function_name=function_name, parallelism=parallelism, nohost=nohost, extraAccArgs=extraAccArgs)
     #:set omp_directive = OMP_ROUTINE(function_name=function_name, nohost=nohost, extraOmpArgs=extraOmpArgs)
 
-    #:if cray_inline == True
+    #:if cray_noinline == True
+        #:if not isinstance(function_name, str)
+            #:stop "When using cray_noinline, function name must be given and given as a string"
+        #:endif
+        #:set cray_noinline_directive = ('!DIR$ NOINLINE ' + function_name).strip('\n')
+#ifdef _CRAYFTN
+#if MFC_OpenACC
+        $:acc_directive
+#elif MFC_OpenMP
+        $:omp_directive
+#else
+        $:cray_noinline_directive
+#endif
+#elif MFC_OpenACC
+        $:acc_directive
+#elif MFC_OpenMP
+        $:omp_directive
+#endif
+    #:elif cray_inline == True
         #:if not isinstance(function_name, str)
             #:stop "When inlining for Cray Compiler, function name must be given and given as a string"
         #:endif
         #:set cray_directive = ('!DIR$ INLINEALWAYS ' + function_name).strip('\n')
 #ifdef _CRAYFTN
+#if MFC_OpenACC
+        $:acc_directive
+#elif MFC_OpenMP
+        $:omp_directive
+#else
         $:cray_directive
+#endif
 #elif MFC_OpenACC
         $:acc_directive
 #elif MFC_OpenMP

src/common/m_chemistry.fpp

@@ -63,7 +63,11 @@ contains
 
         integer :: x, y, z, eqn
         real(wp) :: energy, T_in
-        real(wp), dimension(num_species) :: Ys
+        #:if USING_CCE
+            real(wp), dimension(10) :: Ys
+        #:else
+            real(wp), dimension(num_species) :: Ys
+        #:endif
 
         do z = bounds(3)%beg, bounds(3)%end
             do y = bounds(2)%beg, bounds(2)%end
@@ -101,7 +105,11 @@ contains
         type(int_bounds_info), dimension(1:3), intent(in) :: bounds
 
         integer :: x, y, z, i
-        real(wp), dimension(num_species) :: Ys
+        #:if USING_CCE
+            real(wp), dimension(10) :: Ys
+        #:else
+            real(wp), dimension(num_species) :: Ys
+        #:endif
         real(wp) :: mix_mol_weight
 
         do z = bounds(3)%beg, bounds(3)%end
@@ -131,7 +139,7 @@ contains
         integer :: eqn
         real(wp) :: T
         real(wp) :: rho, omega_m
-        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
+        #:if (not MFC_CASE_OPTIMIZATION and USING_AMD) or USING_CCE
             real(wp), dimension(10) :: Ys
             real(wp), dimension(10) :: omega
         #:else
@@ -180,7 +188,7 @@ contains
         type(int_bounds_info), intent(in) :: irx, iry, irz
 
         integer, intent(in) :: idir
-        #:if not MFC_CASE_OPTIMIZATION and USING_AMD
+        #:if (not MFC_CASE_OPTIMIZATION and USING_AMD) or USING_CCE
             real(wp), dimension(10) :: Xs_L, Xs_R, Xs_cell, Ys_L, Ys_R, Ys_cell
             real(wp), dimension(10) :: mass_diffusivities_mixavg1, mass_diffusivities_mixavg2
             real(wp), dimension(10) :: mass_diffusivities_mixavg_Cell, dXk_dxi, h_l, h_r, h_k

src/common/m_phase_change.fpp

@@ -104,6 +104,16 @@ contains
         !< Generic loop iterators
         integer :: i, j, k, l
 
+#ifdef _CRAYFTN
+#ifdef MFC_OpenACC
+        ! CCE 19 IPA workaround: prevent bring_routine_resident SIGSEGV
+        !DIR$ NOINLINE s_infinite_pt_relaxation_k
+        !DIR$ NOINLINE s_infinite_ptg_relaxation_k
+        !DIR$ NOINLINE s_correct_partial_densities
+        !DIR$ NOINLINE s_TSat
+#endif
+#endif
+
         ! starting equilibrium solver
         $:GPU_PARALLEL_LOOP(collapse=3, private='[i,j,k,l,p_infOV, p_infpT, p_infSL, sk, hk, gk, ek, rhok,pS, pSOV, pSSL, TS, TSOV, TSatOV, TSatSL, TSSL, rhoe, dynE, rhos, rho, rM, m1, m2, MCT, TvF]')
         do j = 0, m
@@ -296,7 +306,7 @@ contains
         !!  @param TS equilibrium temperature at the interface
     subroutine s_infinite_pt_relaxation_k(j, k, l, MFL, pS, p_infpT, q_cons_vf, rhoe, TS)
         $:GPU_ROUTINE(function_name='s_infinite_pt_relaxation_k', &
-            & parallelism='[seq]', cray_inline=True)
+            & parallelism='[seq]', cray_noinline=True)
 
         ! initializing variables
         integer, intent(in) :: j, k, l, MFL
@@ -411,7 +421,7 @@ contains
         !!  @param TS equilibrium temperature at the interface
     subroutine s_infinite_ptg_relaxation_k(j, k, l, pS, p_infpT, rhoe, q_cons_vf, TS)
         $:GPU_ROUTINE(function_name='s_infinite_ptg_relaxation_k', &
-            & parallelism='[seq]', cray_inline=True)
+            & parallelism='[seq]', cray_noinline=True)
 
         integer, intent(in) :: j, k, l
         real(wp), intent(inout) :: pS
@@ -579,7 +589,8 @@ contains
             InvJac = InvJac/(Jac(1, 1)*Jac(2, 2) - Jac(1, 2)*Jac(2, 1))
 
             ! calculating correction array for Newton's method
-            DeltamP = -1.0_wp*(matmul(InvJac, R2D))
+            DeltamP(1) = -1.0_wp*(InvJac(1, 1)*R2D(1) + InvJac(1, 2)*R2D(2))
+            DeltamP(2) = -1.0_wp*(InvJac(2, 1)*R2D(1) + InvJac(2, 2)*R2D(2))
 
             ! updating two reacting 'masses'. Recall that inert 'masses' do not change during the phase change
             ! liquid
@@ -638,7 +649,7 @@ contains
         !!  @param l generic loop iterator for z direction
     subroutine s_correct_partial_densities(MCT, q_cons_vf, rM, j, k, l)
         $:GPU_ROUTINE(function_name='s_correct_partial_densities', &
-            & parallelism='[seq]', cray_inline=True)
+            & parallelism='[seq]', cray_noinline=True)
 
         !> @name variables for the correction of the reacting partial densities
         !> @{
@@ -689,7 +700,7 @@ contains
         !!  @param TSIn equilibrium Temperature
     elemental subroutine s_TSat(pSat, TSat, TSIn)
         $:GPU_ROUTINE(function_name='s_TSat',parallelism='[seq]', &
-            & cray_inline=True)
+            & cray_noinline=True)
 
         real(wp), intent(in) :: pSat
         real(wp), intent(out) :: TSat
@@ -716,6 +727,12 @@ contains
 
             ! underrelaxation factor
             Om = 1.0e-3_wp
+
+            ! FT must be initialized before the do while condition is evaluated.
+            ! Fortran .or. is not short-circuit: abs(FT) is always evaluated even
+            ! when ns == 0, so FT must have a defined value here.
+            FT = huge(1.0_wp)
+
             do while ((abs(FT) > ptgalpha_eps) .or. (ns == 0))
                 ! increasing counter
                 ns = ns + 1

src/common/m_variables_conversion.fpp

@@ -116,7 +116,7 @@ contains
         !! @param pres_mag Magnetic pressure (optional)
     subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoYks, pres, T, stress, mom, G, pres_mag)
         $:GPU_ROUTINE(function_name='s_compute_pressure',parallelism='[seq]', &
-            & cray_inline=True)
+            & cray_noinline=True)
 
         real(stp), intent(in) :: energy, alf
         real(wp), intent(in) :: dyn_p
@@ -326,7 +326,7 @@ contains
                                                           alpha_K, alpha_rho_K, Re_K, &
                                                           G_K, G)
         $:GPU_ROUTINE(function_name='s_convert_species_to_mixture_variables_acc', &
-            & parallelism='[seq]', cray_inline=True)
+            & parallelism='[seq]', cray_noinline=True)
 
         real(wp), intent(out) :: rho_K, gamma_K, pi_inf_K, qv_K
         #:if not MFC_CASE_OPTIMIZATION and USING_AMD
@@ -1335,7 +1335,7 @@ contains
     !>  This subroutine computes partial densities and volume fractions
     subroutine s_compute_species_fraction(q_vf, k, l, r, alpha_rho_K, alpha_K)
         $:GPU_ROUTINE(function_name='s_compute_species_fraction', &
-            & parallelism='[seq]', cray_inline=True)
+            & parallelism='[seq]', cray_noinline=True)
         type(scalar_field), dimension(sys_size), intent(in) :: q_vf
         integer, intent(in) :: k, l, r
         #:if not MFC_CASE_OPTIMIZATION and USING_AMD
@@ -1480,7 +1480,7 @@ contains
     !> @brief Computes the fast magnetosonic wave speed from the sound speed, density, and magnetic field components.
     subroutine s_compute_fast_magnetosonic_speed(rho, c, B, norm, c_fast, h)
         $:GPU_ROUTINE(function_name='s_compute_fast_magnetosonic_speed', &
-            & parallelism='[seq]', cray_inline=True)
+            & parallelism='[seq]', cray_noinline=True)
 
         real(wp), intent(in) :: B(3), rho, c
         real(wp), intent(in) :: h ! only used for relativity

src/simulation/m_bubbles_EL.fpp

@@ -1532,8 +1532,7 @@ contains
         integer(KIND=MPI_OFFSET_KIND) :: disp
         integer :: view
         integer, dimension(2) :: gsizes, lsizes, start_idx_part
-        integer, dimension(num_procs) :: part_order, part_ord_mpi
-        integer, dimension(num_procs) :: proc_bubble_counts
+        integer, allocatable :: proc_bubble_counts(:)
         real(wp), dimension(1:1, 1:lag_io_vars) :: dummy
         dummy = 0._wp
 
@@ -1548,6 +1547,8 @@ contains
 
         if (.not. parallel_io) return
 
+        allocate (proc_bubble_counts(num_procs))
+
         lsizes(1) = bub_id
         lsizes(2) = lag_io_vars
 
@@ -1659,6 +1660,8 @@ contains
             call MPI_FILE_CLOSE(ifile, ierr)
         end if
 
+        deallocate (proc_bubble_counts)
+
 #endif
 
     end subroutine s_write_restart_lag_bubbles

toolchain/mfc/run/input.py

@@ -90,13 +90,25 @@ def generate_fpp(self, target) -> None:
             directive_str = None
 
         # Write the generated Fortran code to the m_thermochem.f90 file with the chosen precision
+        thermochem_code = pyro.FortranCodeGenerator().generate(
+            "m_thermochem",
+            self.get_cantera_solution(),
+            pyro.CodeGenerationOptions(scalar_type = real_type, directive_offload = directive_str)
+        )
+
+        # CCE 19.0.0 workaround: pyrometheus generates !DIR$ INLINEALWAYS for Cray+ACC
+        # but omits !$acc routine seq, so thermochem routines are not registered as
+        # OpenACC device routines. Replace with plain !$acc routine seq (no INLINEALWAYS).
+        if directive_str == 'acc':
+            thermochem_code = thermochem_code.replace(
+                "#ifdef _CRAYFTN\n#define GPU_ROUTINE(name) !DIR$ INLINEALWAYS name\n"
+                "#else\n#define GPU_ROUTINE(name) !$acc routine seq\n#endif",
+                "#define GPU_ROUTINE(name) !$acc routine seq"
+            )
+
         common.file_write(
             os.path.join(modules_dir, "m_thermochem.f90"),
-            pyro.FortranCodeGenerator().generate(
-                "m_thermochem",
-                self.get_cantera_solution(),
-                pyro.CodeGenerationOptions(scalar_type = real_type, directive_offload = directive_str)
-            ),
+            thermochem_code,
             True
         )