GPU workarounds for OLCF Frontier

Author: henryleberre

CMakeLists.txt

@@ -133,6 +133,8 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
             -fimplicit-none
             #-ffpe-trap=invalid,zero,denormal,overflow
             -fsignaling-nans
+            -finit-real=snan
+            -finit-integer=-99999999
 	    )
     endif()
 

src/common/m_mpi_common.fpp

@@ -236,8 +236,8 @@ contains
         real(kind(0d0)), intent(out) :: ccfl_max_glb
         real(kind(0d0)), intent(out) :: Rc_min_glb
 
-#ifdef MFC_MPI
 #ifdef MFC_SIMULATION
+#ifdef MFC_MPI
 
         ! Reducing local extrema of ICFL, VCFL, CCFL and Rc numbers to their
         ! global extrema and bookkeeping the results on the rank 0 processor
@@ -254,6 +254,15 @@ contains
                             MPI_COMM_WORLD, ierr)
         end if
 
+#else
+
+        icfl_max_glb = icfl_max_loc
+
+        if (any(Re_size > 0)) then
+            vcfl_max_glb = vcfl_max_loc
+            Rc_min_glb = Rc_min_loc
+        end if
+
 #endif
 #endif
 

src/common/m_phase_change.fpp

@@ -298,7 +298,11 @@ contains
         !!  @param rhoe mixture energy
         !!  @param TS equilibrium temperature at the interface
     subroutine s_infinite_pt_relaxation_k(j, k, l, MFL, pS, p_infpT, rM, q_cons_vf, rhoe, TS)
+#ifdef CRAY_ACC_WAR
+        !DIR$ INLINEALWAYS s_compute_speed_of_sound
+#else
         !$acc routine seq
+#endif
 
         ! initializing variables
         integer, intent(in) :: j, k, l, MFL
@@ -402,7 +406,11 @@ 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)
 
+#ifdef CRAY_ACC_WAR
+        !DIR$ INLINEALWAYS s_infinite_ptg_relaxation_k
+#else
         !$acc routine seq
+#endif
 
         integer, intent(in) :: j, k, l
         real(kind(0.0d0)), intent(inout) :: pS
@@ -522,7 +530,11 @@ contains
         !!  @param k generic loop iterator for y direction
         !!  @param l generic loop iterator for z direction
     subroutine s_correct_partial_densities(MCT, q_cons_vf, rM, j, k, l)
+#ifdef CRAY_ACC_WAR
+        !DIR$ INLINEALWAYS s_correct_partial_densities
+#else
         !$acc routine seq
+#endif
 
         !> @name variables for the correction of the reacting partial densities
         !> @{
@@ -580,7 +592,12 @@ contains
         !!  @param q_cons_vf Cell-average conservative variables
         !!  @param TJac Transpose of the Jacobian Matrix
     subroutine s_compute_jacobian_matrix(InvJac, j, Jac, k, l, mCPD, mCVGP, mCVGP2, pS, q_cons_vf, TJac)
+
+#ifdef CRAY_ACC_WAR
+        !DIR$ INLINEALWAYS s_compute_jacobian_matrix
+#else
         !$acc routine seq
+#endif
 
         real(kind(0.0d0)), dimension(2, 2), intent(out) :: InvJac
         integer, intent(in) :: j
@@ -682,7 +699,12 @@ contains
         !!  @param rhoe mixture energy
         !!  @param R2D (2D) residue array
     subroutine s_compute_pTg_residue(j, k, l, mCPD, mCVGP, mQD, q_cons_vf, pS, rhoe, R2D)
+
+#ifdef CRAY_ACC_WAR
+        !DIR$ INLINEALWAYS s_compute_pTg_residue
+#else
         !$acc routine seq
+#endif
 
         integer, intent(in) :: j, k, l
         real(kind(0.0d0)), intent(in) :: mCPD, mCVGP, mQD
@@ -728,7 +750,11 @@ contains
         !!  @param TSat Saturation Temperature
         !!  @param TSIn equilibrium Temperature
     subroutine s_TSat(pSat, TSat, TSIn)
+#ifdef CRAY_ACC_WAR
+        !DIR$ INLINEALWAYS s_compute_speed_of_sound
+#else
         !$acc routine seq
+#endif
 
         real(kind(0.0d0)), intent(in) :: pSat
         real(kind(0.0d0)), intent(out) :: TSat

src/common/m_variables_conversion.fpp

@@ -126,7 +126,11 @@ contains
         !! @param mom Momentum
     subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoYks, pres, T, stress, mom, G)
 
+#ifdef CRAY_ACC_WAR
+        !DIR$ INLINEALWAYS s_compute_pressure
+#else
         !$acc routine seq
+#endif
 
         real(kind(0d0)), intent(in) :: energy, alf
         real(kind(0d0)), intent(in) :: dyn_p
@@ -135,7 +139,6 @@ contains
         real(kind(0d0)), intent(in), optional :: stress, mom, G
 
         ! Chemistry
-        integer :: i
         real(kind(0d0)), dimension(1:num_species), intent(in) :: rhoYks
         real(kind(0d0)) :: E_e
         real(kind(0d0)) :: e_Per_Kg, Pdyn_Per_Kg
@@ -182,11 +185,8 @@ contains
             end if
 
         #:else
-            !$acc loop seq
-            do i = 1, num_species
-                Y_rs(i) = rhoYks(i)/rho
-            end do
 
+            Y_rs(:) = rhoYks(:)/rho
             e_Per_Kg = energy/rho
             Pdyn_Per_Kg = dyn_p/rho
 
@@ -886,7 +886,9 @@ contains
             end if
         #:endif
 
-        !$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, Re_K, nRtmp, rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K, R3tmp, rhoyks)
+        !$acc parallel loop collapse(3) gang vector default(present) &
+        !$acc private(alpha_K, alpha_rho_K, Re_K, nRtmp, rho_K, gamma_K, &
+        !$acc pi_inf_K, qv_K, dyn_pres_K, R3tmp, rhoYks)
         do l = ibounds(3)%beg, ibounds(3)%end
             do k = ibounds(2)%beg, ibounds(2)%end
                 do j = ibounds(1)%beg, ibounds(1)%end
@@ -927,7 +929,6 @@ contains
                         rho_K = 0d0
                         !$acc loop seq
                         do i = chemxb, chemxe
-                            !print*, j,k,l, qK_cons_vf(i)%sf(j, k, l)
                             rho_K = rho_K + max(0d0, qK_cons_vf(i)%sf(j, k, l))
                         end do
 
@@ -936,16 +937,9 @@ contains
                             qK_prim_vf(i)%sf(j, k, l) = rho_K
                         end do
 
-                        Yksum = 0d0
                         !$acc loop seq
                         do i = chemxb, chemxe
                             qK_prim_vf(i)%sf(j, k, l) = max(0d0, qK_cons_vf(i)%sf(j, k, l)/rho_K)
-                            Yksum = Yksum + qK_prim_vf(i)%sf(j, k, l)
-                        end do
-
-                        !$acc loop seq
-                        do i = chemxb, chemxe
-                            qK_prim_vf(i)%sf(j, k, l) = qK_prim_vf(i)%sf(j, k, l)/Yksum
                         end do
 
                         qK_prim_vf(tempxb)%sf(j, k, l) = qK_cons_vf(tempxb)%sf(j, k, l)

src/simulation/m_chemistry.fpp

@@ -63,16 +63,16 @@ contains
             if (num_dims == ${num_dims}$) then
                 !$acc parallel loop collapse(4) gang vector default(present) &
                 !$acc private(flux_x, flux_y, flux_z)
-                do x = 0, m
-                    do y = 0, n
-                        do z = 0, p
+                do z = idwint(3)%beg, idwint(3)%end
+                    do y = idwint(2)%beg, idwint(2)%end
+                        do x = idwint(1)%beg, idwint(1)%end
                             do eqn = chemxb, chemxe
                                 ! \nabla \cdot (F)
                                 flux_x = (flux_n(1)%vf(eqn)%sf(x - 1, y, z) - &
                                           flux_n(1)%vf(eqn)%sf(x, y, z))/dx(x)
 
                                 #:if num_dims >= 2
-                                    flux_x = (flux_n(2)%vf(eqn)%sf(x, y - 1, z) - &
+                                    flux_y = (flux_n(2)%vf(eqn)%sf(x, y - 1, z) - &
                                               flux_n(2)%vf(eqn)%sf(x, y, z))/dy(y)
                                 #:else
                                     flux_y = 0d0
@@ -87,6 +87,8 @@ contains
 
                                 rhs_vf(eqn)%sf(x, y, z) = flux_x + flux_y + flux_z
                             end do
+
+                            rhs_vf(tempxb)%sf(x, y, z) = 0d0
                         end do
                     end do
                 end do
@@ -112,38 +114,27 @@ contains
         real(kind(0d0)) :: rho, omega_m
         real(kind(0d0)), dimension(num_species) :: Ys
         real(kind(0d0)), dimension(num_species) :: omega
-        real(kind(0d0)), dimension(num_species) :: enthalpies
         real(kind(0d0)) :: cp_mix
 
         #:if chemistry
 
-            !$acc parallel loop collapse(3) private(Ys)
-            do x = 0, m
-                do y = 0, n
-                    do z = 0, p
-
-                        rho = 0d0
-                        do eqn = chemxb, chemxe
-                            rho = rho + q_cons_qp(eqn)%sf(x, y, z)
-                        end do
+            !$acc parallel loop collapse(3) gang vector default(present) &
+            !$acc private(Ys, omega)
+            do z = idwint(3)%beg, idwint(3)%end
+                do y = idwint(2)%beg, idwint(2)%end
+                    do x = idwint(1)%beg, idwint(1)%end
 
+                        !$acc loop seq
                         do eqn = chemxb, chemxe
                             Ys(eqn - chemxb + 1) = q_prim_qp(eqn)%sf(x, y, z)
                         end do
 
-                        dyn_pres = 0d0
-                        !$acc loop seq
-                        do i = momxb, momxe
-                            dyn_pres = dyn_pres + q_cons_qp(i)%sf(x, y, z)* &
-                                       q_prim_qp(i)%sf(x, y, z)/2d0
-                        end do
-
-                        call get_temperature(q_cons_qp(E_idx)%sf(x, y, z)/rho - dyn_pres/rho, &
-                            & 1200d0, Ys, .true., T)
+                        rho = q_cons_qp(contxe)%sf(x, y, z)
+                        T = q_prim_qp(tempxb)%sf(x, y, z)
 
                         call get_net_production_rates(rho, T, Ys, omega)
 
-                        !$acc routine seq
+                        !$acc loop seq
                         do eqn = chemxb, chemxe
 
                             omega_m = mol_weights(eqn - chemxb + 1)*omega(eqn - chemxb + 1)

src/simulation/m_data_output.fpp

@@ -348,16 +348,13 @@ contains
                 call s_mpi_abort('ICFL is greater than 1.0. Exiting ...')
             end if
 
-            do i = chemxb, chemxe
-                !@:ASSERT(all(q_prim_vf(i)%sf(:,:,:) >= -1d0), "bad conc")
-                !@:ASSERT(all(q_prim_vf(i)%sf(:,:,:) <=  2d0), "bad conc")
-            end do
-
-            if (vcfl_max_glb /= vcfl_max_glb) then
-                call s_mpi_abort('VCFL is NaN. Exiting ...')
-            elseif (vcfl_max_glb > 1d0) then
-                print *, 'vcfl', vcfl_max_glb
-                call s_mpi_abort('VCFL is greater than 1.0. Exiting ...')
+            if (any(Re_size > 0)) then
+                if (vcfl_max_glb /= vcfl_max_glb) then
+                    call s_mpi_abort('VCFL is NaN. Exiting ...')
+                elseif (vcfl_max_glb > 1d0) then
+                    print *, 'vcfl', vcfl_max_glb
+                    call s_mpi_abort('VCFL is greater than 1.0. Exiting ...')
+                end if
             end if
         end if
 
@@ -510,7 +507,7 @@ contains
         if (.not. file_exist) call s_create_directory(trim(t_step_dir))
 
         if (prim_vars_wrt .or. (n == 0 .and. p == 0)) then
-            call s_convert_conservative_to_primitive_variables(q_cons_vf, q_prim_vf, idwbuff)
+            call s_convert_conservative_to_primitive_variables(q_cons_vf, q_prim_vf, idwint)
             do i = 1, sys_size
                 !$acc update host(q_prim_vf(i)%sf(:,:,:))
             end do

src/simulation/m_global_parameters.fpp

@@ -164,6 +164,8 @@ module m_global_parameters
     logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling
     logical :: cu_tensor
 
+    !$acc declare create(chemistry)
+
     logical :: bodyForces
     logical :: bf_x, bf_y, bf_z !< body force toggle in three directions
     !< amplitude, frequency, and phase shift sinusoid in each direction

src/simulation/m_rhs.fpp

@@ -958,15 +958,15 @@ contains
             call nvtxEndRange
             ! END: Additional physics and source terms =========================
 
-            #:if chemistry
-                call nvtxStartRange("RHS_Chem_Advection")
-                call s_compute_chemistry_advection_flux(flux_n, rhs_vf)
-                call nvtxEndRange
-            #:endif
-
         end do
         ! END: Dimensional Splitting Loop =================================
 
+        #:if chemistry
+            call nvtxStartRange("RHS_Chem_Advection")
+            call s_compute_chemistry_advection_flux(flux_n, rhs_vf)
+            call nvtxEndRange
+        #:endif
+
         if (ib) then
             !$acc parallel loop collapse(3) gang vector default(present)
             do l = 0, p

src/simulation/m_riemann_solvers.fpp

@@ -363,8 +363,11 @@ contains
         #:for NORM_DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')]
 
             if (norm_dir == ${NORM_DIR}$) then
-                !$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, vel_avg, tau_e_L, tau_e_R, G_L, G_R, Re_L, Re_R, &
-                !$acc rho_avg, h_avg, gamma_avg, s_L, s_R, s_S, Y_L, Ys_L, Y_R, Ys_R)
+                !$acc parallel loop collapse(3) gang vector default(present)    &
+                !$acc private(alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L,  &
+                !$acc alpha_R, vel_avg, tau_e_L, tau_e_R, G_L, G_R, Re_L, Re_R, &
+                !$acc rho_avg, h_avg, gamma_avg, s_L, s_R, s_S, Y_L, Ys_L, Y_R, &
+                !$acc Ys_R)
                 do l = is3%beg, is3%end
                     do k = is2%beg, is2%end
                         do j = is1%beg, is1%end

src/simulation/m_sim_helpers.f90

@@ -31,7 +31,12 @@ module m_sim_helpers
         !! @param k y index
         !! @param l z index
     subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l)
+#ifdef CRAY_ACC_WAR
+        !DIR$ INLINEALWAYS s_compute_enthalpy
+#else
         !$acc routine seq
+#endif
+
         type(scalar_field), dimension(sys_size) :: q_prim_vf
         real(kind(0d0)), dimension(num_fluids) :: alpha_rho
         real(kind(0d0)), dimension(num_fluids) :: alpha