add option to dynamically read pyro mech from input file

githooks/pre-commit

@@ -4,7 +4,7 @@ set -o pipefail
 
 echo "Running flake8..."
 if [[ $(command -v "flake8") ]]; then
-    flake8 --ignore=E126,E127,E128,E123,E226,E241,E242,E265,N802,W503,E402,N803,N806,N814,N817,W504 --show-source --statistics setup.py driver.py y3prediction 
+    flake8 --ignore=E126,E127,E128,E123,E226,E241,E242,E265,N802,W503,E402,N803,N806,N814,N817,W504 --show-source --statistics setup.py driver.py y3prediction/*.py
     res=$?
     if [[ $res -ne 0 ]]; then
         echo "Error: flake8 check failed. Fix the errors (or run git with --no-verify to bypass the check)."

y3prediction/prediction.py

@@ -122,7 +122,6 @@
     WallVars,
     WallModel,
 )
-from y3prediction.uiuc_sharp import Thermochemistry
 from y3prediction.shock1d import PlanarDiscontinuityMulti
 
 from dataclasses import dataclass
@@ -1144,6 +1143,14 @@ def _compiled_stepper_wrapper(state, t, dt, rhs):
             av_mu=av2_mu0, av_beta=av2_beta0, av_kappa=av2_kappa0,
             av_prandtl=av2_prandtl0)
 
+    pyro_mech = configurate("pyro_mech", input_data, "uiuc_sharp")
+    pyro_mech_name = f"y3prediction.pyro_mechs.{pyro_mech}"
+
+    import importlib
+    pyromechlib = importlib.import_module(pyro_mech_name)
+
+    #from y3prediction.uiuc_sharp import Thermochemistry
+
     chem_source_tol = 1.e-10
     # make the eos
     if eos_type == 0:
@@ -1153,13 +1160,13 @@ def _compiled_stepper_wrapper(state, t, dt, rhs):
     else:
         from mirgecom.thermochemistry import get_pyrometheus_wrapper_class
         pyro_mech = get_pyrometheus_wrapper_class(
-            pyro_class=Thermochemistry, temperature_niter=pyro_temp_iter,
+            pyro_class=pyromechlib.Thermochemistry, temperature_niter=pyro_temp_iter,
             zero_level=chem_source_tol)(actx.np)
         eos = PyrometheusMixture(pyro_mech, temperature_guess=init_temperature)
         # seperate gas model for initialization,
         # just to make sure we get converged temperature
         pyro_mech_init = get_pyrometheus_wrapper_class(
-            pyro_class=Thermochemistry, temperature_niter=5,
+            pyro_class=pyromechlib.Thermochemistry, temperature_niter=5,
             zero_level=chem_source_tol)(actx.np)
         eos_init = PyrometheusMixture(pyro_mech_init,
                                       temperature_guess=init_temperature)

y3prediction/pyro_mechs/check_mechanism.py

@@ -0,0 +1,79 @@
+import cantera as ct
+import numpy as np
+from matplotlib import pyplot as plt
+
+gas = ct.Solution("my_mechanism.yaml")
+gas.TP = 300, 101325
+gas.X = "C2H4:1"
+
+numpts = 43
+visc = np.zeros(numpts,)
+enthalpy = np.zeros(numpts,)
+cp = np.zeros(numpts,)
+cv = np.zeros(numpts,)
+temperature_list = np.linspace(300, 2400, numpts)
+ii = 0
+for temp in temperature_list:
+    gas.TP = temp, 101325
+    visc[ii] = gas.viscosity
+    enthalpy[ii] = gas.enthalpy_mass
+    cp[ii] = gas.cp_mass
+    cv[ii] = gas.cv_mass
+    ii += 1
+
+gamma = cp/cv
+
+plt.close("all")
+fig = plt.figure(1, figsize=[6.4, 4.8])
+ax1 = fig.add_subplot(111)
+ax1.set_position([0.16, 0.12, 0.75, 0.83])
+ax1.plot(temperature_list, visc, "--", label="Current")
+ax1.set_ylabel(r"$\mathbf{Viscosity \, (Pa-s)}$")
+ax1.set_xlabel(r"$\mathbf{Temperature \, (K)}$")
+plt.savefig("viscosity.png", dpi=200)
+ax1.legend()
+plt.show()
+
+plt.close("all")
+fig = plt.figure(1, figsize=[6.4, 4.8])
+ax1 = fig.add_subplot(111)
+ax1.set_position([0.16, 0.12, 0.75, 0.83])
+ax1.plot(temperature_list, enthalpy, "--", label="Current")
+ax1.set_ylabel(r"$\mathbf{Enthalpy \, (J/kg)}$")
+ax1.set_xlabel(r"$\mathbf{Temperature \, (K)}$")
+plt.savefig("enthalpy.png", dpi=200)
+ax1.legend()
+plt.show()
+
+plt.close("all")
+fig = plt.figure(1, figsize=[6.4, 4.8])
+ax1 = fig.add_subplot(111)
+ax1.set_position([0.16, 0.12, 0.75, 0.83])
+ax1.plot(temperature_list, cp, "--", label="Current")
+ax1.set_ylabel(r"$\mathbf{Heat \; Capacity (Cp)\, (J/kg-K)}$")
+ax1.set_xlabel(r"$\mathbf{Temperature \, (K)}$")
+plt.savefig("heatCapacityCp.png", dpi=200)
+ax1.legend()
+plt.show()
+
+plt.close("all")
+fig = plt.figure(1, figsize=[6.4, 4.8])
+ax1 = fig.add_subplot(111)
+ax1.set_position([0.16, 0.12, 0.75, 0.83])
+ax1.plot(temperature_list, cv, "--", label="Current")
+ax1.set_ylabel(r"$\mathbf{Heat \; Capacity (Cv)\, (J/kg-K)}$")
+ax1.set_xlabel(r"$\mathbf{Temperature \, (K)}$")
+plt.savefig("heatCapacityCv.png", dpi=200)
+ax1.legend()
+plt.show()
+plt.close("all")
+
+fig = plt.figure(1, figsize=[6.4, 4.8])
+ax1 = fig.add_subplot(111)
+ax1.set_position([0.16, 0.12, 0.75, 0.83])
+ax1.plot(temperature_list, gamma, "--", label="Current")
+ax1.set_ylabel(r"$\mathbf{gamma}$")
+ax1.set_xlabel(r"$\mathbf{Temperature \, (K)}$")
+plt.savefig("gamma.png", dpi=200)
+ax1.legend()
+plt.show()

y3prediction/pyro_mechs/generate_pyro_mech.py

@@ -0,0 +1,45 @@
+import cantera as ct
+import pyrometheus as pyro
+
+
+def generate_mechfile(mech_file):
+    """This produces the mechanism codes."""
+
+    sol = ct.Solution(f"{mech_file}", "gas")
+    from pathlib import Path
+
+    mech_output_file = Path(f"{mech_file}").stem
+    with open(f"{mech_output_file}.py", "w") as file:
+        code = pyro.codegen.python.gen_thermochem_code(sol)
+        print(code, file=file)
+
+
+import logging
+import sys
+import argparse
+
+if __name__ == "__main__":
+
+    logging.basicConfig(
+        format="%(asctime)s - %(levelname)s - %(name)s - %(message)s",
+        level=logging.INFO)
+
+    parser = argparse.ArgumentParser(
+        description="Generate Pyrometheus mechanisms")
+    parser.add_argument("-f", "--mech_file", type=ascii, dest="mech_file",
+                        nargs="?", action="store", help="mechanism file")
+
+    args = parser.parse_args()
+
+    # get mechanism name from the arguments
+    from mirgecom.simutil import ApplicationOptionsError
+    mech_file = ""
+    if args.mech_file:
+        print(f"Processing Cantera mechanism file {args.mech_file}")
+        mech_file = args.mech_file.replace("'", "")
+    else:
+        raise ApplicationOptionsError("Missing Cantera mechanism file from input")
+
+    print(f"Running {sys.argv[0]}\n")
+
+    generate_mechfile(mech_file)