From 114f54876813e0709a8cd3e150d2a7e6d9362665 Mon Sep 17 00:00:00 2001 From: "terence.lehmann" Date: Thu, 14 Nov 2024 18:16:02 +0100 Subject: [PATCH 1/3] ceptr: Added functionality for reading and converting PLOG reactions. --- Support/ceptr/ceptr/jacobian.py | 16 +++++++++-- Support/ceptr/ceptr/production.py | 32 +++++++++++++++------ Support/ceptr/ceptr/utilities.py | 47 +++++++++++++++++++++++++++++++ 3 files changed, 85 insertions(+), 10 deletions(-) diff --git a/Support/ceptr/ceptr/jacobian.py b/Support/ceptr/ceptr/jacobian.py index 592ec83da..8cab6dff2 100644 --- a/Support/ceptr/ceptr/jacobian.py +++ b/Support/ceptr/ceptr/jacobian.py @@ -562,12 +562,13 @@ def ajac_reaction_d( else: dim = cu.phase_space_units(reaction.reactants) third_body = reaction.third_body is not None + plog = reaction.rate.type == "pressure-dependent-Arrhenius" falloff = reaction.rate.type == "falloff" is_troe = reaction.rate.sub_type == "Troe" is_sri = reaction.rate.sub_type == "Sri" is_lindemann = reaction.rate.sub_type == "Lindemann" aeuc = cu.activation_energy_units() - if not third_body and not falloff: + if not third_body and not falloff and not plog: # Case 3 !PD, !TB cw.writer( fstream, @@ -577,7 +578,7 @@ def ajac_reaction_d( pef = (reaction.rate.pre_exponential_factor * ctuc).to_base_units() beta = reaction.rate.temperature_exponent ae = (reaction.rate.activation_energy * cc.ureg.joule / cc.ureg.kmol).to(aeuc) - elif not falloff: + elif not falloff and not plog: # Case 2 !PD, TB cw.writer( fstream, @@ -587,6 +588,17 @@ def ajac_reaction_d( pef = (reaction.rate.pre_exponential_factor * ctuc).to_base_units() beta = reaction.rate.temperature_exponent ae = (reaction.rate.activation_energy * cc.ureg.joule / cc.ureg.kmol).to(aeuc) + elif not third_body and not falloff and plog: + # Case 4 PLOG + cw.writer( + fstream, + 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) + pef = (plog_pef * ctuc).to_base_units() + beta = plog_beta + ae = (plog_ae * cc.ureg.joule / cc.ureg.kmol).to(aeuc) else: # Case 1 PD, TB cw.writer( diff --git a/Support/ceptr/ceptr/production.py b/Support/ceptr/ceptr/production.py index d59f241a9..2854cc637 100644 --- a/Support/ceptr/ceptr/production.py +++ b/Support/ceptr/ceptr/production.py @@ -188,27 +188,35 @@ def production_rate( else: dim = cu.phase_space_units(reaction.reactants) third_body = reaction.third_body is not None + plog = reaction.rate.type == "pressure-dependent-Arrhenius" falloff = reaction.rate.type == "falloff" is_troe = reaction.rate.sub_type == "Troe" is_sri = reaction.rate.sub_type == "Sri" is_lindemann = reaction.rate.sub_type == "Lindemann" aeuc = cu.activation_energy_units() - if not third_body and not falloff: - # Case 3 !PD, !TB + if not third_body and not falloff and not plog: + # Case 3 !PD, !TB, !PLOG ctuc = cu.prefactor_units(cc.ureg("kmol/m**3"), 1 - dim) pef = (reaction.rate.pre_exponential_factor * ctuc).to_base_units() beta = reaction.rate.temperature_exponent ae = ( reaction.rate.activation_energy * cc.ureg.joule / cc.ureg.kmol ).to(aeuc) - elif not falloff: - # Case 2 !PD, TB + elif not falloff and not plog: + # Case 2 !PD, TB, !PLOG ctuc = cu.prefactor_units(cc.ureg("kmol/m**3"), -dim) pef = (reaction.rate.pre_exponential_factor * ctuc).to_base_units() beta = reaction.rate.temperature_exponent ae = ( reaction.rate.activation_energy * cc.ureg.joule / cc.ureg.kmol ).to(aeuc) + 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) + pef = (plog_pef * ctuc).to_base_units() + beta = plog_beta + ae = (plog_ae * cc.ureg.joule / cc.ureg.kmol).to(aeuc) else: # Case 1 PD, TB ctuc = cu.prefactor_units(cc.ureg("kmol/m**3"), 1 - dim) @@ -553,12 +561,13 @@ def production_rate( dim = cu.phase_space_units(reaction.reactants) third_body = reaction.third_body is not None + plog = reaction.rate.type == "pressure-dependent-Arrhenius" falloff = reaction.rate.type == "falloff" is_troe = reaction.rate.sub_type == "Troe" is_sri = reaction.rate.sub_type == "Sri" is_lindemann = reaction.rate.sub_type == "Lindemann" aeuc = cu.activation_energy_units() - if not third_body and not falloff: + if not third_body and not falloff and not plog: # Case 3 !PD, !TB ctuc = cu.prefactor_units(cc.ureg("kmol/m**3"), 1 - dim) pef = (reaction.rate.pre_exponential_factor * ctuc).to_base_units() @@ -566,16 +575,23 @@ def production_rate( ae = ( reaction.rate.activation_energy * cc.ureg.joule / cc.ureg.kmol ).to(aeuc) - elif not falloff: - # Case 2 !PD, TB + elif not falloff and not plog: + # Case 2 !PD, TB, !PLOG ctuc = cu.prefactor_units(cc.ureg("kmol/m**3"), -dim) pef = (reaction.rate.pre_exponential_factor * ctuc).to_base_units() beta = reaction.rate.temperature_exponent ae = ( reaction.rate.activation_energy * cc.ureg.joule / cc.ureg.kmol ).to(aeuc) + 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) + pef = (plog_pef * ctuc).to_base_units() + beta = plog_beta + ae = (plog_ae * cc.ureg.joule / cc.ureg.kmol).to(aeuc) else: - # Case 1 PD, TB + # Case 1 PD, TB, !PLOG ctuc = cu.prefactor_units(cc.ureg("kmol/m**3"), 1 - dim) pef = ( reaction.rate.high_rate.pre_exponential_factor * ctuc diff --git a/Support/ceptr/ceptr/utilities.py b/Support/ceptr/ceptr/utilities.py index 555dfab75..577007fd6 100644 --- a/Support/ceptr/ceptr/utilities.py +++ b/Support/ceptr/ceptr/utilities.py @@ -3,9 +3,13 @@ import copy from collections import Counter from math import isclose +from math import log +from math import exp import ceptr.constants as cc +# Global constant for plog evaluation +plog_pressure = None def intersection(lst1, lst2): """Return intersection of two lists.""" @@ -451,3 +455,46 @@ def enhancement_d(mechanism, species_info, reaction, syms=None): return " + ".join(alpha).replace("+ -", "- "), enhancement_smp else: return " + ".join(alpha).replace("+ -", "- ") + + +def evaluate_plog(rates): + """Evaluation the 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: + sys.exit("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] + pef = plog_reaction.pre_exponential_factor + beta = plog_reaction.temperature_exponent + ae = plog_reaction.activation_energy + return pef, beta, ae + elif plog_pressure >= rates[-1][0]: + # Case 2: plog_pressure >= upper bound -> take upper bound: + plog_reaction = rates[-1][1] + pef = plog_reaction.pre_exponential_factor + beta = plog_reaction.temperature_exponent + ae = plog_reaction.activation_energy + return pef, beta, ae + else: + # Case 3: lower bound < plog_pressure < upper bound -> logarithmic interpolation: + for plog_p_i in range(len(rates) - 1): + if rates[plog_p_i][0] <= plog_pressure < rates[plog_p_i + 1][0]: + # print(f"Interpolation between {rates[plog_p_i][0]} Pa and {rates[plog_p_i + 1][0]} Pa.") + rate_low = rates[plog_p_i][1] + rate_high = rates[plog_p_i + 1][1] + interp_fac = (log(plog_pressure) - log(rates[plog_p_i][0])) / ( + log(rates[plog_p_i + 1][0]) - log(rates[plog_p_i][0])) + pef = exp(log(rate_low.pre_exponential_factor) * (1 - interp_fac) + log(rate_high.pre_exponential_factor) * interp_fac) + beta = rate_low.temperature_exponent * (1 - interp_fac) + rate_high.temperature_exponent * interp_fac + ae = rate_low.activation_energy * (1 - interp_fac) + rate_high.activation_energy * interp_fac + return pef, beta, ae + From eac86a022bd35768f237e116dee19fcae4118c3c Mon Sep 17 00:00:00 2001 From: "terence.lehmann" Date: Thu, 14 Nov 2024 23:13:36 +0100 Subject: [PATCH 2/3] Do formatting with black --- Support/ceptr/ceptr/utilities.py | 38 +++++++++++++++++++++++--------- 1 file changed, 27 insertions(+), 11 deletions(-) diff --git a/Support/ceptr/ceptr/utilities.py b/Support/ceptr/ceptr/utilities.py index 577007fd6..0b2b71a9f 100644 --- a/Support/ceptr/ceptr/utilities.py +++ b/Support/ceptr/ceptr/utilities.py @@ -11,6 +11,7 @@ # Global constant for plog evaluation plog_pressure = None + def intersection(lst1, lst2): """Return intersection of two lists.""" return list(set(lst1).intersection(lst2)) @@ -462,14 +463,20 @@ def evaluate_plog(rates): # 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")) + 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: - sys.exit("Provided plog_pressure too low.") # To catch confusion about Pascal and bar/atm + sys.exit( + "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] @@ -492,9 +499,18 @@ def evaluate_plog(rates): rate_low = rates[plog_p_i][1] rate_high = rates[plog_p_i + 1][1] interp_fac = (log(plog_pressure) - log(rates[plog_p_i][0])) / ( - log(rates[plog_p_i + 1][0]) - log(rates[plog_p_i][0])) - pef = exp(log(rate_low.pre_exponential_factor) * (1 - interp_fac) + log(rate_high.pre_exponential_factor) * interp_fac) - beta = rate_low.temperature_exponent * (1 - interp_fac) + rate_high.temperature_exponent * interp_fac - ae = rate_low.activation_energy * (1 - interp_fac) + rate_high.activation_energy * interp_fac + log(rates[plog_p_i + 1][0]) - log(rates[plog_p_i][0]) + ) + pef = exp( + log(rate_low.pre_exponential_factor) * (1 - interp_fac) + + log(rate_high.pre_exponential_factor) * interp_fac + ) + beta = ( + rate_low.temperature_exponent * (1 - interp_fac) + + rate_high.temperature_exponent * interp_fac + ) + ae = ( + rate_low.activation_energy * (1 - interp_fac) + + rate_high.activation_energy * interp_fac + ) return pef, beta, ae - From b38e6cb4bca6c66b7c8e98d4a22f823de5a6d5ca Mon Sep 17 00:00:00 2001 From: "terence.lehmann" Date: Thu, 14 Nov 2024 23:30:03 +0100 Subject: [PATCH 3/3] Run isort and flake8 --- Support/ceptr/ceptr/utilities.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/Support/ceptr/ceptr/utilities.py b/Support/ceptr/ceptr/utilities.py index 0b2b71a9f..1036a1f59 100644 --- a/Support/ceptr/ceptr/utilities.py +++ b/Support/ceptr/ceptr/utilities.py @@ -2,9 +2,7 @@ import copy from collections import Counter -from math import isclose -from math import log -from math import exp +from math import exp, isclose, log import ceptr.constants as cc @@ -459,7 +457,7 @@ def enhancement_d(mechanism, species_info, reaction, syms=None): def evaluate_plog(rates): - """Evaluation the rate constants for a PLOG reaction""" + """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: @@ -474,7 +472,7 @@ def evaluate_plog(rates): ) print(f"plog_pressure set to {plog_pressure} Pa / {plog_pressure / 1e5} bar /.") if plog_pressure < 1e3: - sys.exit( + raise ValueError( "Provided plog_pressure too low." ) # To catch confusion about Pascal and bar/atm if plog_pressure <= rates[0][0]: @@ -495,7 +493,6 @@ def evaluate_plog(rates): # Case 3: lower bound < plog_pressure < upper bound -> logarithmic interpolation: for plog_p_i in range(len(rates) - 1): if rates[plog_p_i][0] <= plog_pressure < rates[plog_p_i + 1][0]: - # print(f"Interpolation between {rates[plog_p_i][0]} Pa and {rates[plog_p_i + 1][0]} Pa.") rate_low = rates[plog_p_i][1] rate_high = rates[plog_p_i + 1][1] interp_fac = (log(plog_pressure) - log(rates[plog_p_i][0])) / (