Dummy Chemistry Mechanism (#680) (#689)

henryleberre · web-flow · commit c36ce91bf691 · 2024-11-06T09:09:32.000-05:00
diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp
@@ -9,7 +9,8 @@
 module m_derived_types
 
     use m_constants  !< Constants
-    use m_thermochem !< Thermodynamic properties
+
+    use m_thermochem, only: num_species
 
     implicit none
 
diff --git a/src/common/m_thermochem.fpp b/src/common/m_thermochem.fpp
diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp
@@ -22,7 +22,9 @@ module m_variables_conversion
 
     use m_helper
 
-    use m_thermochem
+    use m_thermochem, only: &
+        num_species, get_temperature, get_pressure, &
+        get_mixture_molecular_weight, get_mixture_energy_mass
 
     ! ==========================================================================
 
@@ -1126,7 +1128,7 @@ contains
                                    q_prim_vf(i)%sf(j, k, l)/2d0
                     end do
 
-                    #:if chemistry
+                    if (chemistry) then
                         do i = chemxb, chemxe
                             Ys(i - chemxb + 1) = q_prim_vf(i)%sf(j, k, l)
                             q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l)
@@ -1140,7 +1142,7 @@ contains
                             dyn_pres + rho*e_mix
 
                         q_cons_vf(T_idx)%sf(j, k, l) = q_prim_vf(T_idx)%sf(j, k, l)
-                    #:else
+                    else
                         ! Computing the energy from the pressure
                         if ((model_eqns /= 4) .and. (bubbles .neqv. .true.)) then
                             ! E = Gamma*P + \rho u u /2 + \pi_inf + (\alpha\rho qv)
@@ -1156,7 +1158,7 @@ contains
                             !Tait EOS, no conserved energy variable
                             q_cons_vf(E_idx)%sf(j, k, l) = 0.
                         end if
-                    #:endif
+                    end if
 
                     ! Computing the internal energies from the pressure and continuities
                     if (model_eqns == 3) then
diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp
@@ -18,7 +18,7 @@ module m_global_parameters
 
     use m_helper_basic          !< Functions to compare floating point numbers
 
-    use m_thermochem            !< Thermodynamic and chemical properties module
+    use m_thermochem, only: num_species, species_names
 
     ! ==========================================================================
 
diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90
@@ -35,8 +35,7 @@ module m_start_up
 
     use m_checker
 
-    use m_thermochem            !< Procedures used to compute thermodynamic
-                                !! quantities
+    use m_thermochem, only: num_species, species_names
 
     use m_finite_differences
 
diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp
@@ -15,7 +15,7 @@ module m_assign_variables
 
     use m_helper_basic          !< Functions to compare floating point numbers
 
-    use m_thermochem            !< Thermodynamic and chemical properties
+    use m_thermochem, only: num_species, gas_constant, get_mixture_molecular_weight
 
     ! one form to another
     ! ==========================================================================
@@ -166,7 +166,7 @@ contains
             + (1d0 - eta)*patch_icpp(smooth_patch_id)%pi_inf
 
         ! Species Concentrations
-        #:if chemistry
+        if (chemistry) then
             block
                 real(kind(0d0)) :: sum, term
 
@@ -194,7 +194,7 @@ contains
             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
+        end if
 
         ! Updating the patch identities bookkeeping variable
         if (1d0 - eta < 1d-16) patch_id_fp(j, k, l) = patch_id
@@ -543,7 +543,7 @@ contains
         end do
 
         ! Species Concentrations
-        #:if chemistry
+        if (chemistry) then
             block
                 real(kind(0d0)) :: sum, term
 
@@ -572,7 +572,7 @@ contains
             call get_mixture_molecular_weight(Ys, mean_molecular_weight)
             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
+        end if
 
         ! Set streamwise velocity to hyperbolic tangent function of y
         if (mixlayer_vel_profile) then
diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp
@@ -29,6 +29,9 @@ module m_data_output
     use m_helper
 
     use m_delay_file_access
+
+    use m_thermochem, only: species_names
+
     ! ==========================================================================
 
     implicit none
diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp
@@ -18,7 +18,7 @@ module m_global_parameters
 
     use m_helper_basic          ! Functions to compare floating point numbers
 
-    use m_thermochem            ! Thermodynamic and chemical properties
+    use m_thermochem, only: num_species
 
     ! ==========================================================================
 
diff --git a/src/simulation/include/inline_riemann.fpp b/src/simulation/include/inline_riemann.fpp
@@ -31,7 +31,7 @@
     vel_avg_rms = (sqrt(rho_L)*vel_L(1) + sqrt(rho_R)*vel_R(1))**2d0/ &
                   (sqrt(rho_L) + sqrt(rho_R))**2d0
 
-    #:if chemistry
+    if (chemistry) then
         eps = 0.001d0
         call get_species_enthalpies_rt(T_L, h_iL)
         call get_species_enthalpies_rt(T_R, h_iR)
@@ -59,7 +59,7 @@
 
         Phi_avg(:) = (gamma_avg - 1.d0)*(vel_avg_rms/2.0d0 - h_avg_2(:)) + gamma_avg*gas_constant/mol_weights(:)*T_avg
         c_sum_Yi_Phi = sum(Yi_avg(:)*Phi_avg(:))
-    #:endif
+    end if
 
 #:enddef roe_avg
 
diff --git a/src/simulation/m_chemistry.fpp b/src/simulation/m_chemistry.fpp
@@ -11,8 +11,11 @@ module m_chemistry
     use ieee_arithmetic
 
     use m_mpi_proxy
-    use m_thermochem
+    use m_thermochem, only: &
+        num_species, mol_weights, get_net_production_rates
+
     use m_global_parameters
+
     use m_finite_differences
 
     implicit none
@@ -116,8 +119,7 @@ contains
         real(kind(0d0)), dimension(num_species) :: omega
         real(kind(0d0)) :: cp_mix
 
-        #:if chemistry
-
+        if (chemistry) then
             !$acc parallel loop collapse(3) gang vector default(present) &
             !$acc private(Ys, omega)
             do z = idwint(3)%beg, idwint(3)%end
@@ -147,11 +149,11 @@ contains
                 end do
             end do
 
-        #:else
+        else
 
             @:ASSERT(.false., "Chemistry is not enabled")
 
-        #:endif
+        end if
 
     end subroutine s_compute_chemistry_reaction_flux
 
diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp
@@ -106,8 +106,6 @@ 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)
diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp
@@ -961,11 +961,11 @@ contains
         end do
         ! END: Dimensional Splitting Loop =================================
 
-        #:if chemistry
+        if (chemistry) then
             call nvtxStartRange("RHS_Chem_Advection")
             call s_compute_chemistry_advection_flux(flux_n, rhs_vf)
             call nvtxEndRange
-        #:endif
+        end if
 
         if (ib) then
             !$acc parallel loop collapse(3) gang vector default(present)
@@ -1000,13 +1000,13 @@ contains
             rhs_vf)
         call nvtxEndRange
 
-        #:if chemistry
+        if (chemistry) then
             if (chem_params%reactions) then
                 call nvtxStartRange("RHS_Chem_Reactions")
                 call s_compute_chemistry_reaction_flux(rhs_vf, q_cons_qp%vf, q_prim_qp%vf)
                 call nvtxEndRange
             end if
-        #:endif
+        end if
 
         ! END: Additional pphysics and source terms ============================
 
diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp
@@ -38,6 +38,12 @@ module m_riemann_solvers
     use m_surface_tension      !< To get the capilary fluxes
 
     use m_chemistry
+
+    use m_thermochem, only: &
+        gas_constant, get_mixture_molecular_weight, &
+        get_mixture_specific_heat_cv_mass, get_mixture_energy_mass, &
+        get_species_specific_heats_r, get_species_enthalpies_rt, &
+        get_mixture_specific_heat_cp_mass
     ! ==========================================================================
 
     implicit none
@@ -483,7 +489,7 @@ contains
                                 end do
                             end if
 
-                            #:if chemistry
+                            if (chemistry) then
                                 !$acc loop seq
                                 do i = chemxb, chemxe
                                     Ys_L(i - chemxb + 1) = qL_prim_rs${XYZ}$_vf(j, k, l, i)
@@ -532,12 +538,12 @@ contains
                                 E_R = rho_R*E_R + 5d-1*rho_R*vel_R_rms
                                 H_L = (E_L + pres_L)/rho_L
                                 H_R = (E_R + pres_R)/rho_R
-                            #:else
+                            else
                                 E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms + qv_L
                                 E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms + qv_R
                                 H_L = (E_L + pres_L)/rho_L
                                 H_R = (E_R + pres_R)/rho_R
-                            #:endif
+                            end if
 
                             if (hypoelasticity) then
                                 !$acc loop seq
@@ -1002,8 +1008,10 @@ contains
                 if (model_eqns == 3) then
                     !ME3
 
-                    !$acc parallel loop collapse(3) gang vector default(present) private(vel_L, vel_R, Re_L, Re_R, &
-                    !$acc rho_avg, h_avg, gamma_avg, s_L, s_R, s_S, vel_avg_rms, alpha_L, alpha_R)
+                    !$acc parallel loop collapse(3) gang vector default(present) &
+                    !$acc private(vel_L, vel_R, Re_L, Re_R, rho_avg, h_avg, gamma_avg, &
+                    !$acc s_L, s_R, s_S, vel_avg_rms, alpha_L, alpha_R, Ys_L, Ys_R, Xs_L, Xs_R, &
+                    !$acc Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2)
 
                     do l = is3%beg, is3%end
                         do k = is2%beg, is2%end
@@ -2182,7 +2190,7 @@ contains
                                     end do
                                 end if
 
-                                #:if chemistry
+                                if (chemistry) then
                                     c_sum_Yi_Phi = 0.0d0
                                     !$acc loop seq
                                     do i = chemxb, chemxe
@@ -2232,14 +2240,14 @@ contains
                                     E_R = rho_R*E_R + 5d-1*rho_R*vel_R_rms
                                     H_L = (E_L + pres_L)/rho_L
                                     H_R = (E_R + pres_R)/rho_R
-                                #:else
+                                else
                                     E_L = gamma_L*pres_L + pi_inf_L + 5d-1*rho_L*vel_L_rms + qv_L
 
                                     E_R = gamma_R*pres_R + pi_inf_R + 5d-1*rho_R*vel_R_rms + qv_R
 
                                     H_L = (E_L + pres_L)/rho_L
                                     H_R = (E_R + pres_R)/rho_R
-                                #:endif
+                                end if
 
                                 @:compute_average_state()
 
diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp
@@ -37,7 +37,7 @@ module m_time_steppers
 
     use m_nvtx
 
-    use m_thermochem
+    use m_thermochem, only: num_species
 
     use m_body_forces
     ! ==========================================================================
diff --git a/toolchain/mfc/run/input.py b/toolchain/mfc/run/input.py
@@ -34,13 +34,19 @@ def __save_fpp(self, target, contents: str) -> None:
         common.file_write(fpp_path, contents, True)
 
     def get_cantera_solution(self) -> ct.Solution:
-        cantera_file = self.params["cantera_file"]
-
-        candidates = [
-            cantera_file,
-            os.path.join(self.dirpath, cantera_file),
-            os.path.join(common.MFC_MECHANISMS_DIR, cantera_file),
-        ]
+        if self.params.get("chemistry", 'F') == 'T':
+            cantera_file = self.params["cantera_file"]
+
+            candidates = [
+                cantera_file,
+                os.path.join(self.dirpath, cantera_file),
+                os.path.join(common.MFC_MECHANISMS_DIR, cantera_file),
+            ]
+        else:
+            # If Chemistry is turned off, we return a default (dummy) solution
+            # that will not be used in the simulation, so that MFC can still
+            # be compiled.
+            candidates = ["h2o2.yaml"]
 
         for candidate in candidates:
             try:
@@ -60,19 +66,17 @@ def generate_fpp(self, target) -> None:
         # Case FPP file
         self.__save_fpp(target, self.get_fpp(target))
 
-        # Chemistry Rates FPP file
+        # (Thermo)Chemistry source file
         modules_dir = os.path.join(target.get_staging_dirpath(self), "modules", target.name)
         common.create_directory(modules_dir)
-
-        if self.params.get("chemistry", 'F') == 'T':
-            common.file_write(
-                os.path.join(modules_dir, "m_pyrometheus.f90"),
-                pyro.codegen.fortran90.gen_thermochem_code(
-                    self.get_cantera_solution(),
-                    module_name="m_pyrometheus"
-                ),
-                True
-            )
+        common.file_write(
+            os.path.join(modules_dir, "m_thermochem.f90"),
+            pyro.codegen.fortran90.gen_thermochem_code(
+                self.get_cantera_solution(),
+                module_name="m_thermochem"
+            ),
+            True
+        )
 
         cons.unindent()
 
