Address PR comments

src/common/m_derived_types.fpp

@@ -314,6 +314,9 @@ module m_derived_types
         logical :: diffusion
         logical :: reactions
 
+        !> Method of determining gamma.
+        !> gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
+        !> gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
         integer :: gamma_method
     end type chemistry_parameters
 

src/common/m_variables_conversion.fpp

@@ -942,7 +942,7 @@ contains
                             qK_prim_vf(i)%sf(j, k, l) = max(0d0, qK_cons_vf(i)%sf(j, k, l)/rho_K)
                         end do
 
-                        qK_prim_vf(tempxb)%sf(j, k, l) = qK_cons_vf(tempxb)%sf(j, k, l)
+                        qK_prim_vf(T_idx)%sf(j, k, l) = qK_cons_vf(T_idx)%sf(j, k, l)
                     else
                         !$acc loop seq
                         do i = 1, contxe
@@ -980,7 +980,7 @@ contains
 
                     qK_prim_vf(E_idx)%sf(j, k, l) = pres
                     if (chemistry) then
-                        qK_prim_vf(tempxb)%sf(j, k, l) = T
+                        qK_prim_vf(T_idx)%sf(j, k, l) = T
                     end if
 
                     if (bubbles) then
@@ -1133,13 +1133,13 @@ contains
                         end do
 
                         call get_mixture_molecular_weight(Ys, mix_mol_weight)
-                        T = q_prim_vf(tempxb)%sf(j, k, l)
+                        T = q_prim_vf(T_idx)%sf(j, k, l)
                         call get_mixture_energy_mass(T, Ys, e_mix)
 
                         q_cons_vf(E_idx)%sf(j, k, l) = &
                             dyn_pres + rho*e_mix
 
-                        q_cons_vf(tempxb)%sf(j, k, l) = q_prim_vf(tempxb)%sf(j, k, l)
+                        q_cons_vf(T_idx)%sf(j, k, l) = q_prim_vf(T_idx)%sf(j, k, l)
                     #:else
                         ! Computing the energy from the pressure
                         if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then
@@ -1398,7 +1398,7 @@ contains
     end subroutine s_finalize_variables_conversion_module
 
 #ifndef MFC_PRE_PROCESS
-    subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_avggg, c)
+    subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_sum, c_c, c)
 #ifdef CRAY_ACC_WAR
         !DIR$ INLINEALWAYS s_compute_speed_of_sound
 #else
@@ -1409,16 +1409,14 @@ contains
         real(kind(0d0)), intent(in) :: H
         real(kind(0d0)), dimension(num_fluids), intent(in) :: adv
         real(kind(0d0)), intent(in) :: vel_sum
-        real(kind(0d0)), intent(in) :: c_avggg
+        real(kind(0d0)), intent(in) :: c_c
         real(kind(0d0)), intent(out) :: c
 
         real(kind(0d0)) :: blkmod1, blkmod2
-        real(kind(0d0)) :: Tolerance, c_c
+        real(kind(0d0)) :: Tolerance
         integer :: q
 
         if (chemistry) then
-            c_c = c_avggg
-
             if (avg_state == 1 .and. abs(c_c) > Tolerance) then
                 c = sqrt(c_c - (gamma - 1.0d0)*(vel_sum - H))
             else

src/post_process/m_global_parameters.fpp

@@ -129,14 +129,16 @@ module m_global_parameters
     type(int_bounds_info) :: stress_idx            !< Indices of elastic stresses
     integer :: c_idx                               !< Index of color function
     type(int_bounds_info) :: species_idx           !< Indexes of first & last concentration eqns.
-    type(int_bounds_info) :: temperature_idx       !< Indexes of first & last temperature eqns.
+    integer :: T_idx                               !< Index of temperature eqn.
     !> @}
 
     ! Cell Indices for the (local) interior points (O-m, O-n, 0-p).
+    ! Stands for "InDices With BUFFer".
     type(int_bounds_info) :: idwint(1:3)
 
     ! Cell Indices for the entire (local) domain. In simulation, this includes
     ! the buffer region. idwbuff and idwint are the same otherwise.
+    ! Stands for "InDices With BUFFer".
     type(int_bounds_info) :: idwbuff(1:3)
 
     !> @name Boundary conditions in the x-, y- and z-coordinate directions
@@ -284,7 +286,6 @@ module m_global_parameters
     integer :: bubxb, bubxe
     integer :: strxb, strxe
     integer :: chemxb, chemxe
-    integer :: tempxb, tempxe
     !> @}
 
 contains
@@ -624,14 +625,12 @@ contains
             species_idx%end = sys_size + num_species
             sys_size = species_idx%end
 
-            temperature_idx%beg = sys_size + 1
-            temperature_idx%end = sys_size + 1
-            sys_size = temperature_idx%end
+            T_idx = sys_size + 1
+            sys_size = T_idx
         else
             species_idx%beg = 1
             species_idx%end = 1
-            temperature_idx%beg = 1
-            temperature_idx%end = 1
+            T_idx = 1
         end if
 
         momxb = mom_idx%beg
@@ -648,8 +647,6 @@ contains
         intxe = internalEnergies_idx%end
         chemxb = species_idx%beg
         chemxe = species_idx%end
-        tempxb = temperature_idx%beg
-        tempxe = temperature_idx%end
 
 #ifdef MFC_MPI
         allocate (MPI_IO_DATA%view(1:sys_size))

src/post_process/m_start_up.f90

@@ -312,9 +312,9 @@ subroutine s_save_data(t_step, varname, pres, c, H)
             end do
 
             if (chem_wrt_T) then
-                q_sf = q_prim_vf(tempxb)%sf(-offset_x%beg:m + offset_x%end, &
-                                            -offset_y%beg:n + offset_y%end, &
-                                            -offset_z%beg:p + offset_z%end)
+                q_sf = q_prim_vf(T_idx)%sf(-offset_x%beg:m + offset_x%end, &
+                                           -offset_y%beg:n + offset_y%end, &
+                                           -offset_z%beg:p + offset_z%end)
 
                 write (varname, '(A)') 'T'
                 call s_write_variable_to_formatted_database_file(varname, t_step)

src/pre_process/m_assign_variables.fpp

@@ -191,7 +191,7 @@ contains
             end block
 
             call get_mixture_molecular_weight(Ys, mean_molecular_weight)
-            q_prim_vf(tempxb)%sf(j, k, l) = &
+            q_prim_vf(T_idx)%sf(j, k, l) = &
                 q_prim_vf(E_idx)%sf(j, k, l)*mean_molecular_weight &
                 /(gas_constant*q_prim_vf(1)%sf(j, k, l))
         #:endif
@@ -570,7 +570,7 @@ contains
             end block
 
             call get_mixture_molecular_weight(Ys, mean_molecular_weight)
-            q_prim_vf(tempxb)%sf(j, k, l) = &
+            q_prim_vf(T_idx)%sf(j, k, l) = &
                 q_prim_vf(E_idx)%sf(j, k, l)*mean_molecular_weight/(gas_constant*q_prim_vf(1)%sf(j, k, l))
         #:endif
 

src/pre_process/m_data_output.fpp

@@ -261,7 +261,7 @@ contains
                                  .or. &
                                  ((i >= chemxb) .and. (i <= chemxe)) &
                                  .or. &
-                                 ((i == tempxb)) &
+                                 ((i == T_idx)) &
                                  ) then
                             write (2, FMT) x_cb(j), q_cons_vf(i)%sf(j, 0, 0)
                         else if (i == mom_idx%beg) then !u
@@ -825,7 +825,7 @@ contains
                 write (1, '(I3,A20,A20)') chemxb + i - 1, "Y_{"//trim(species_names(i))//"} \rho", "Y_{"//trim(species_names(i))//"}"
             end do
 
-            write (1, '(I3,A20,A20)') tempxb, "T", "T"
+            write (1, '(I3,A20,A20)') T_idx, "T", "T"
         end if
 
         write (1, '(A)') ""
@@ -837,7 +837,7 @@ contains
         if (strxb /= 0) write (1, '("[",I2,",",I2,"]",A)') strxb, strxe, " Stress"
         if (intxb /= 0) write (1, '("[",I2,",",I2,"]",A)') intxb, intxe, " Internal Energies"
         if (chemxb /= 0) write (1, '("[",I2,",",I2,"]",A)') chemxb, chemxe, " Chemistry"
-        if (tempxb /= 0) write (1, '("[",I2,",",I2,"]",A)') tempxb, tempxe, " Temperature"
+        if (T_idx /= 0) write (1, '("[",I2,",",I2,"]",A)') T_idx, T_idx, " Temperature"
 
         close (1)
 

src/pre_process/m_global_parameters.fpp

@@ -106,13 +106,15 @@ module m_global_parameters
     type(int_bounds_info) :: stress_idx            !< Indexes of elastic shear stress eqns.
     integer :: c_idx                               !< Index of the color function
     type(int_bounds_info) :: species_idx           !< Indexes of first & last concentration eqns.
-    type(int_bounds_info) :: temperature_idx       !< Indexes of first & last temperature eqns.
+    integer :: T_idx                               !< Index of temperature eqn.
 
     ! Cell Indices for the (local) interior points (O-m, O-n, 0-p).
+    ! Stands for "InDices With BUFFer".
     type(int_bounds_info) :: idwint(1:3)
 
     ! Cell Indices for the entire (local) domain. In simulation and post_process,
     ! this includes the buffer region. idwbuff and idwint are the same otherwise.
+    ! Stands for "InDices With BUFFer".
     type(int_bounds_info) :: idwbuff(1:3)
 
     type(int_bounds_info) :: bc_x, bc_y, bc_z !<
@@ -237,7 +239,6 @@ module m_global_parameters
     integer :: bubxb, bubxe
     integer :: strxb, strxe
     integer :: chemxb, chemxe
-    integer :: tempxb, tempxe
     !> @}
 
     integer, allocatable, dimension(:, :, :) :: logic_grid
@@ -707,9 +708,8 @@ contains
             species_idx%end = sys_size + num_species
             sys_size = species_idx%end
 
-            temperature_idx%beg = sys_size + 1
-            temperature_idx%end = sys_size + 1
-            sys_size = temperature_idx%end
+            T_idx = sys_size + 1
+            sys_size = T_idx
         end if
 
         momxb = mom_idx%beg
@@ -726,8 +726,6 @@ contains
         intxe = internalEnergies_idx%end
         chemxb = species_idx%beg
         chemxe = species_idx%end
-        tempxb = temperature_idx%beg
-        tempxe = temperature_idx%end
 
         ! Configuring Coordinate Direction Indexes =========================
         idwint(1)%beg = 0; idwint(2)%beg = 0; idwint(3)%beg = 0

src/simulation/include/inline_riemann.fpp

@@ -58,7 +58,7 @@
         gamma_avg = Cp_avg/Cv_avg
 
         Phi_avg(:) = (gamma_avg - 1.d0)*(vel_avg_rms/2.0d0 - h_avg_2(:)) + gamma_avg*gas_constant/mol_weights(:)*T_avg
-        c_avggg = sum(Yi_avg(:)*Phi_avg(:))
+        c_sum_Yi_Phi = sum(Yi_avg(:)*Phi_avg(:))
     #:endif
 
 #:enddef roe_avg

src/simulation/m_chemistry.fpp

@@ -88,7 +88,7 @@ contains
                                 rhs_vf(eqn)%sf(x, y, z) = flux_x + flux_y + flux_z
                             end do
 
-                            rhs_vf(tempxb)%sf(x, y, z) = 0d0
+                            rhs_vf(T_idx)%sf(x, y, z) = 0d0
                         end do
                     end do
                 end do
@@ -130,7 +130,7 @@ contains
                         end do
 
                         rho = q_cons_qp(contxe)%sf(x, y, z)
-                        T = q_prim_qp(tempxb)%sf(x, y, z)
+                        T = q_prim_qp(T_idx)%sf(x, y, z)
 
                         call get_net_production_rates(rho, T, Ys, omega)
 

src/simulation/m_global_parameters.fpp

@@ -239,17 +239,19 @@ module m_global_parameters
     type(int_bounds_info) :: stress_idx                !< Indexes of first and last shear stress eqns.
     integer :: c_idx         ! Index of the color function
     type(int_bounds_info) :: species_idx           !< Indexes of first & last concentration eqns.
-    type(int_bounds_info) :: temperature_idx       !< Indexes of first & last temperature eqns.
+    integer :: T_idx       !< Index of the temperature equation
     !> @}
 
     !$acc declare create(bub_idx)
 
     ! Cell Indices for the (local) interior points (O-m, O-n, 0-p).
+    ! Stands for "InDices With INTerior".
     type(int_bounds_info) :: idwint(1:3)
     !$acc declare create(idwint)
 
     ! Cell Indices for the entire (local) domain. In simulation and post_process,
     ! this includes the buffer region. idwbuff and idwint are the same otherwise.
+    ! Stands for "InDices With BUFFer".
     type(int_bounds_info) :: idwbuff(1:3)
     !$acc declare create(idwbuff)
 
@@ -301,7 +303,7 @@ module m_global_parameters
 
     integer :: startx, starty, startz
 
-    !$acc declare create(sys_size, buff_size, startx, starty, startz, E_idx, gamma_idx, pi_inf_idx, alf_idx, n_idx, stress_idx, species_idx)
+    !$acc declare create(sys_size, buff_size, startx, starty, startz, E_idx, T_idx, gamma_idx, pi_inf_idx, alf_idx, n_idx, stress_idx, species_idx)
 
     ! END: Simulation Algorithm Parameters =====================================
 
@@ -467,9 +469,8 @@ module m_global_parameters
     integer :: bubxb, bubxe
     integer :: strxb, strxe
     integer :: chemxb, chemxe
-    integer :: tempxb, tempxe
 
-!$acc declare create(momxb, momxe, advxb, advxe, contxb, contxe, intxb, intxe, bubxb, bubxe, strxb, strxe,  chemxb, chemxe, tempxb, tempxe)
+!$acc declare create(momxb, momxe, advxb, advxe, contxb, contxe, intxb, intxe, bubxb, bubxe, strxb, strxe,  chemxb, chemxe)
 
 #ifdef CRAY_ACC_WAR
     @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:), gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps)
@@ -1004,9 +1005,8 @@ contains
             species_idx%end = sys_size + num_species
             sys_size = species_idx%end
 
-            temperature_idx%beg = sys_size + 1
-            temperature_idx%end = sys_size + 1
-            sys_size = temperature_idx%end
+            T_idx = sys_size + 1
+            sys_size = T_idx
         end if
 
         if (qbmm .and. .not. polytropic) then
@@ -1118,10 +1118,8 @@ contains
         intxe = internalEnergies_idx%end
         chemxb = species_idx%beg
         chemxe = species_idx%end
-        tempxb = temperature_idx%beg
-        tempxe = temperature_idx%end
 
-        !$acc update device(momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe, intxb, intxe, sys_size, buff_size, E_idx, alf_idx, n_idx, adv_n, adap_dt, pi_fac, strxb, strxe, chemxb, chemxe, tempxb, tempxe)
+        !$acc update device(momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe, intxb, intxe, sys_size, buff_size, E_idx, T_idx, alf_idx, n_idx, adv_n, adap_dt, pi_fac, strxb, strxe, chemxb, chemxe)
         !$acc update device(species_idx)
         !$acc update device(cfl_target, m, n, p)
 

src/simulation/m_ibm.fpp

@@ -106,6 +106,8 @@ contains
         @:ALLOCATE_GLOBAL(ghost_points(num_gps))
         @:ALLOCATE_GLOBAL(inner_points(num_inner_gps))
 
+        print *, "num_gps was ", num_gps
+
         !$acc enter data copyin(ghost_points, inner_points)
 
         call s_find_ghost_points(ghost_points, inner_points)
@@ -451,6 +453,8 @@ contains
             :: subsection_3D
         integer :: i, j, k, l, q !< Iterator variables
 
+        num_gps = 0
+
         do i = 0, m
             do j = 0, n
                 if (p == 0) then