Skip to content

Commit

Permalink
Merge pull request #1 from ITV-RWTH/ceptr_development
Browse files Browse the repository at this point in the history 
ceptr: Added functionality for reading and converting PLOG reactions.
  • Loading branch information
@ThomasHowarth
ThomasHowarth authored Nov 15, 2024
2 parents 1e96572 + 8ef83c4 commit e159dda
Show file tree
Hide file tree
Showing 3 changed files with 99 additions and 11 deletions.
16 changes: 14 additions & 2 deletions Support/ceptr/ceptr/jacobian.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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,
Expand All @@ -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(
Expand Down
32 changes: 24 additions & 8 deletions Support/ceptr/ceptr/production.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -553,29 +561,37 @@ 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()
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 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
Expand Down
62 changes: 61 additions & 1 deletion Support/ceptr/ceptr/utilities.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,13 @@

import copy
from collections import Counter
from math import isclose
from math import exp, isclose, log

import ceptr.constants as cc

# Global constant for plog evaluation
plog_pressure = None


def intersection(lst1, lst2):
"""Return intersection of two lists."""
Expand Down Expand Up @@ -451,3 +454,60 @@ def enhancement_d(mechanism, species_info, reaction, syms=None):
return " + ".join(alpha).replace("+ -", "- "), enhancement_smp
else:
return " + ".join(alpha).replace("+ -", "- ")


def evaluate_plog(rates):
"""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]
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]:
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

0 comments on commit e159dda

Please sign in to comment.