Ceptr development

Mechanisms/POLIMI2020/1_000atm/Make.package

@@ -0,0 +1 @@
+CEXE_sources+=mechanism.cpp

Mechanisms/POLIMI2020/Make.package

@@ -1,9 +1 @@
-CEXE_headers+=mechanism.H
-
 CEXE_sources+=mechanism.cpp
-
-LIBRARIES +=
-
-INCLUDE_LOCATIONS += $(CHEM_DIR)
-
-VPATH_LOCATIONS +=  $(CHEM_DIR)/PMFs

Mechanisms/list_mech

@@ -35,4 +35,5 @@ nitrogens/mechanism.yaml
 ndodecane_35/mechanism.yaml
 SootReaction/mechanism.yaml
 sCO2/mechanism.yaml
-H2-CO-CO2-3spec/mechanism.yaml
+H2-CO-CO2-3spec/mechanism.yaml
+POLIMI2020/mechanism.yaml --plog=101325.0

Support/ceptr/ceptr/ceptr.py

@@ -15,11 +15,17 @@ def parse_lst_file(lst):
     """Return mechanism paths give a file containing a list of mechanism files."""
     lpath = pathlib.Path(lst)
     fnames = []
+    plog_pressures = []
     with open(lst, "r") as f:
         for line in f:
             if not line.startswith("#"):
-                fnames.append(line)
-    return [lpath.parents[0] / fn.strip() for fn in fnames]
+                parts = line.split("--plog=")
+                fnames.append(parts[0])
+                if len(parts) > 1:
+                    plog_pressures.append(float(parts[1].strip()))
+                else:
+                    plog_pressures.append(None)
+    return [lpath.parents[0] / fn.strip() for fn in fnames], plog_pressures
 
 
 def parse_qss_lst_file(lst):
@@ -44,6 +50,7 @@ def convert(
     chemistry,
     gas_name,
     interface_name,
+    plog_pressure,
 ):
     """Convert a mechanism file."""
     print(f"""Converting file {fname}""")
@@ -67,6 +74,7 @@ def convert(
         jacobian,
         qss_format_input,
         qss_symbolic_jac,
+        plog_pressure,
     )
     conv.writer()
     conv.formatter()
@@ -83,7 +91,7 @@ def convert_lst(
     interface_name,
 ):
     """Convert mechanisms from a file containing a list of directories."""
-    mechnames = parse_lst_file(lst)
+    mechnames, plog_pressures = parse_lst_file(lst)
     print(f"Using {ncpu} processes")
     with Pool(ncpu) as pool:
         pool.starmap(
@@ -96,6 +104,7 @@ def convert_lst(
                 repeat(chemistry),
                 repeat(gas_name),
                 repeat(interface_name),
+                plog_pressures,
             ),
         )
 
@@ -190,6 +199,14 @@ def main():
         "-n", "--ncpu", help="Number of processes to use", type=int, default=cpu_count()
     )
 
+    parser.add_argument(
+        "-p",
+        "--plog_pressure",
+        help="Pressure in Pascal to evaluate PLOG reactions",
+        type=float,
+        default=None,
+    )
+
     args = parser.parse_args()
 
     if args.chemistry == "heterogeneous":
@@ -206,6 +223,7 @@ def main():
             args.chemistry,
             args.gas_name,
             args.interface_name,
+            args.plog_pressure,
         )
     elif args.lst:
         convert_lst(

Support/ceptr/ceptr/converter.py

@@ -1,12 +1,14 @@
 """Generate C++ files for a mechanism."""
 
+import os
 import pathlib
 import shutil
 import subprocess as spr
 
 import numpy as np
 
 import ceptr.ck as cck
+import ceptr.constants as cc
 import ceptr.formatter as cf
 import ceptr.gjs as cgjs
 import ceptr.jacobian as cj
@@ -32,6 +34,7 @@ def __init__(
         jacobian=True,
         qss_format_input=None,
         qss_symbolic_jacobian=False,
+        plog_pressure=None,
     ):
         self.mechIsAHetMech = chemistry == "heterogeneous"
 
@@ -49,9 +52,6 @@ def __init__(
             else pathlib.Path(self.mechanism.source)
         )
 
-        self.rootname = "mechanism"
-        self.hdrname = self.mechpath.parents[0] / f"{self.rootname}.H"
-        self.cppname = self.mechpath.parents[0] / f"{self.rootname}.cpp"
         self.species_info = csi.SpeciesInfo()
 
         self.set_species()
@@ -70,6 +70,47 @@ def __init__(
         # 6/interface/sticking
         # 6/7        /8
         self.reaction_info = cri.sort_reactions(self.mechanism, self.interface)
+
+        # Set up folder structure for PLOG reactions
+        if self.reaction_info.has_plog_reactions:
+            print(
+                "\nWARNING: Your mechanism contains at least one PLOG reaction.\n"
+                "WARNING: The compiled mechanism will only be valid for the given "
+                "constant pressure. It is not applicable for compressible solvers.\n\n"
+            )
+            # The pressure necessary for the plog evaluation is either provided via
+            # the command line argument or during runtime
+            if plog_pressure is None:  # runtime option
+                plog_pressure = float(
+                    input(
+                        "Please specify the pressure, at which you want to evaluate"
+                        f" the rates in Pascal (1 atm = {cc.Patm_pa} Pa, 1 bar = 1e5"
+                        " Pa):\n"
+                    )
+                )
+            print(
+                f"plog_pressure set to {plog_pressure} Pa /"
+                f" {plog_pressure / 1e5} bar /."
+            )
+            mechanism.TP = mechanism.T, plog_pressure
+            # To catch confusion about Pascal and bar/atm :
+            if mechanism.P < 1e3:
+                raise ValueError("Provided plog_pressure too low.")
+
+            # Now, create the pressure-specific folder
+            plog_folder = f"{plog_pressure/cc.Patm_pa:0.3f}atm".replace(".", "_")
+            if not os.path.isdir(self.mechpath.parents[0] / plog_folder):
+                os.makedirs(self.mechpath.parents[0] / plog_folder)
+                # Copy the Make.package file into the new folder
+                source = self.mechpath.parents[0] / "Make.package"
+                destination = self.mechpath.parents[0] / plog_folder / "Make.package"
+                shutil.copy(source, destination)
+            self.rootname = f"{plog_folder}/mechanism"
+        else:
+            self.rootname = "mechanism"
+        self.hdrname = self.mechpath.parents[0] / f"{self.rootname}.H"
+        self.cppname = self.mechpath.parents[0] / f"{self.rootname}.cpp"
+
         # QSS  -- sort reactions/networks/check validity of QSSs
         if self.species_info.n_qssa_species > 0:
             print("QSSA information")

Support/ceptr/ceptr/jacobian.py

@@ -595,7 +595,9 @@ def ajac_reaction_d(
             cw.comment("a non-third-body and non-pressure-fall-off reaction (PLOG)"),
         )
         ctuc = cu.prefactor_units(cc.ureg("kmol/m**3"), 1 - dim)
-        plog_pef, plog_beta, plog_ae = cu.evaluate_plog(reaction.rate.rates)
+        plog_pef, plog_beta, plog_ae = cu.evaluate_plog(
+            reaction.rate.rates, mechanism.P
+        )
         pef = (plog_pef * ctuc).to_base_units()
         beta = plog_beta
         ae = (plog_ae * cc.ureg.joule / cc.ureg.kmol).to(aeuc)

Support/ceptr/ceptr/production.py

@@ -213,7 +213,9 @@ def production_rate(
             elif plog:
                 # Case 4 PLOG
                 ctuc = cu.prefactor_units(cc.ureg("kmol/m**3"), 1 - dim)
-                plog_pef, plog_beta, plog_ae = cu.evaluate_plog(reaction.rate.rates)
+                plog_pef, plog_beta, plog_ae = cu.evaluate_plog(
+                    reaction.rate.rates, mechanism.P
+                )
                 pef = (plog_pef * ctuc).to_base_units()
                 beta = plog_beta
                 ae = (plog_ae * cc.ureg.joule / cc.ureg.kmol).to(aeuc)
@@ -586,7 +588,9 @@ def production_rate(
             elif not third_body and not falloff and plog:
                 # Case 4 PLOG
                 ctuc = cu.prefactor_units(cc.ureg("kmol/m**3"), 1 - dim)
-                plog_pef, plog_beta, plog_ae = cu.evaluate_plog(reaction.rate.rates)
+                plog_pef, plog_beta, plog_ae = cu.evaluate_plog(
+                    reaction.rate.rates, mechanism.P
+                )
                 pef = (plog_pef * ctuc).to_base_units()
                 beta = plog_beta
                 ae = (plog_ae * cc.ureg.joule / cc.ureg.kmol).to(aeuc)

Support/ceptr/ceptr/reaction_info.py

@@ -27,6 +27,8 @@ def __init__(self, mechanism, interface):
         self.n_qssa_reactions = 0
         self.qfqr_co_idx_map = []
 
+        self.has_plog_reactions = False
+
         if interface is not None:
             self.rs_unsorted += interface.reactions()
 
@@ -86,6 +88,9 @@ def sort_reactions(mechanism, interface):
                 continue
 
             if r.third_body is None:
+                # Check for PLOG reactions on the fly
+                if r.reaction_type == "pressure-dependent-Arrhenius":
+                    reaction_info.has_plog_reactions = True
                 reaction_info.idxmap[k] = i
                 reaction_info.rs.append(r)
                 i += 1

Support/ceptr/ceptr/utilities.py

@@ -6,9 +6,6 @@
 
 import ceptr.constants as cc
 
-# Global constant for plog evaluation
-plog_pressure = None
-
 
 def intersection(lst1, lst2):
     """Return intersection of two lists."""
@@ -456,25 +453,8 @@ def enhancement_d(mechanism, species_info, reaction, syms=None):
         return " + ".join(alpha).replace("+ -", "- ")
 
 
-def evaluate_plog(rates):
+def evaluate_plog(rates, plog_pressure):
     """Evaluate rate constants for a PLOG reaction."""
-    # Ask for pressure on command line if not already done:
-    global plog_pressure
-    if plog_pressure is None:
-        plog_pressure = float(
-            input(
-                "WARNING: Your mechanism contains at least one PLOG reaction.\n"
-                "WARNING: The compiled mechanism will only be valid for the given "
-                "constant pressure. It is not applicable for compressible solvers.\n\n"
-                "Please specify the pressure, at which you want to evaluate the rates "
-                f"in Pascal (1 atm = {cc.Patm_pa} Pa, 1 bar = 1e5 Pa):\n"
-            )
-        )
-        print(f"plog_pressure set to {plog_pressure} Pa / {plog_pressure / 1e5} bar /.")
-        if plog_pressure < 1e3:
-            raise ValueError(
-                "Provided plog_pressure too low."
-            )  # To catch confusion about Pascal and bar/atm
     if plog_pressure <= rates[0][0]:
         # Case 1: plog_pressure <= lower bound -> take lower bound:
         plog_reaction = rates[0][1]