From b8a120ae8bd18af62b37ffaab9cf5fefe9ea5f5d Mon Sep 17 00:00:00 2001 From: Henry LE BERRE Date: Tue, 5 Nov 2024 12:38:14 -0500 Subject: [PATCH] Dummy Chemistry Mechanism (#680) --- src/common/m_derived_types.fpp | 3 +- src/common/m_thermochem.fpp | 21 ------------ src/common/m_variables_conversion.fpp | 10 +++--- src/post_process/m_global_parameters.fpp | 2 +- src/post_process/m_start_up.f90 | 3 +- src/pre_process/m_assign_variables.fpp | 10 +++--- src/pre_process/m_data_output.fpp | 3 ++ src/pre_process/m_global_parameters.fpp | 2 +- src/simulation/include/inline_riemann.fpp | 4 +-- src/simulation/m_chemistry.fpp | 12 ++++--- src/simulation/m_ibm.fpp | 2 -- src/simulation/m_rhs.fpp | 8 ++--- src/simulation/m_riemann_solvers.fpp | 24 +++++++++----- src/simulation/m_time_steppers.fpp | 2 +- toolchain/mfc/run/input.py | 40 +++++++++++++---------- 15 files changed, 71 insertions(+), 75 deletions(-) delete mode 100644 src/common/m_thermochem.fpp diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index d7339ceaa..0ce09bf57 100644 --- 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 deleted file mode 100644 index adee7a15d..000000000 --- a/src/common/m_thermochem.fpp +++ /dev/null @@ -1,21 +0,0 @@ -#:include 'case.fpp' - -module m_thermochem - - #:if chemistry - use m_pyrometheus, only: & - num_species, species_names, gas_constant, mol_weights, & - get_temperature, get_net_production_rates, get_pressure, & - get_mixture_molecular_weight, get_mixture_energy_mass, & - get_mixture_specific_heat_cv_mass, get_mixture_specific_heat_cp_mass, & - get_species_enthalpies_rt, get_species_specific_heats_r - #:endif - - implicit none - - #:if not chemistry - integer, parameter :: num_species = 0 - character(len=:), allocatable, dimension(:) :: species_names - #:endif - -end module m_thermochem diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 7dfc4c07b..c82a15c97 100644 --- 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 index e74f179ce..260b27d3b 100644 --- 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 index 946d6cfc5..90a2da971 100644 --- 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 index f9280d86d..8459f36a5 100644 --- 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 index 0ddc3deaa..bc01e9b84 100644 --- 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 index 4ce08a0c4..c35698ebb 100644 --- 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 index c0805d082..1af67bfb4 100644 --- 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 index 063f38a2e..4465a8a4d 100644 --- 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 index 095cdc0a1..cf6449535 100644 --- 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 index ae541ec64..a415baf21 100644 --- 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 index 6e5123746..9930e3ee6 100644 --- 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 index 97cc9fdac..4f78bb28b 100644 --- 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 index 3437bfdf2..91e65864b 100644 --- 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()