Skip to content

Commit

Permalink
Merge pull request #2015 from ReactionMechanismGenerator/cantera_ther…
Browse files Browse the repository at this point in the history 
…mo_sa

Cantera thermo sa
  • Loading branch information
@mjohnson541
mjohnson541 committed Sep 2, 2020
2 parents 215f9cb + 47dfca1 commit 3c13a22
Show file tree
Hide file tree
Showing 3 changed files with 294 additions and 36 deletions.
247 changes: 228 additions & 19 deletions ipython/cantera_sensitivity_comparison.ipynb
View file

Large diffs are not rendered by default.

79 changes: 64 additions & 15 deletions rmgpy/tools/canteramodel.py
View file
Expand Up @@ -35,7 +35,7 @@
from rmgpy.chemkin import get_species_identifier
from rmgpy.quantity import Quantity
from rmgpy.tools.data import GenericData
from rmgpy.tools.plot import GenericPlot, SimulationPlot, ReactionSensitivityPlot
from rmgpy.tools.plot import GenericPlot, SimulationPlot, ReactionSensitivityPlot, ThermoSensitivityPlot


class CanteraCondition(object):
Expand Down Expand Up @@ -198,14 +198,16 @@ class Cantera(object):
"""

def __init__(self, species_list=None, reaction_list=None, canteraFile='', output_directory='', conditions=None,
sensitive_species=None):
sensitive_species=None, thermo_SA=False):
"""
`species_list`: list of RMG species objects
`reaction_list`: list of RMG reaction objects
`reaction_map`: dict mapping the RMG reaction index within the `reaction_list` to cantera model reaction(s) indices
`canteraFile` path of the chem.cti file associated with this job
`conditions`: a list of `CanteraCondition` objects
`sensitive_species`: a list of RMG species objects for conductng sensitivity analysis on
`thermo_SA`: a boolean indicating whether or not to run thermo SA. By default, if sensitive_species is given,
only kinetic_SA will be calculated and it must be additionally specified to perform thermo SA.
"""
self.species_list = species_list
self.reaction_list = reaction_list
Expand All @@ -214,6 +216,7 @@ def __init__(self, species_list=None, reaction_list=None, canteraFile='', output
self.output_directory = output_directory if output_directory else os.getcwd()
self.conditions = conditions if conditions else []
self.sensitive_species = sensitive_species if sensitive_species else []
self.thermo_SA = thermo_SA

# Make output directory if it does not yet exist:
if not os.path.exists(self.output_directory):
Expand Down Expand Up @@ -327,7 +330,7 @@ def modify_species_thermo(self, rmg_species_index, rmg_species, use_chemkin_iden
ct_species = self.model.species(rmg_species_index)
ct_species.thermo = modified_ct_species.thermo

def plot(self, data, top_species=10, top_sensitive_reactions=10):
def plot(self, data, top_species=10, top_sensitive_reactions=10, top_sensitive_species=10):
"""
Plots data from the simulations from this cantera job.
Takes data in the format of a list of tuples containing (time, [list of temperature, pressure, and species data])
Expand All @@ -342,8 +345,9 @@ def plot(self, data, top_species=10, top_sensitive_reactions=10):
"""
num_ct_reactions = len(self.model.reactions())
num_ct_species = len(self.model.species())
for i, condition_data in enumerate(data):
time, data_list, reaction_sensitivity_data = condition_data
time, data_list, reaction_sensitivity_data, thermodynamic_sensitivity_data = condition_data
# In RMG, any species with an index of -1 is an inert and should not be plotted
inert_list = [species for species in self.species_list if species.index == -1]

Expand All @@ -363,7 +367,13 @@ def plot(self, data, top_species=10, top_sensitive_reactions=10):
ReactionSensitivityPlot(x_var=time,
y_var=reaction_sensitivity_data[j * num_ct_reactions:(j + 1) * num_ct_reactions],
num_reactions=top_sensitive_reactions).barplot(
os.path.join(self.output_directory, '{0}_{1}_sensitivity.png'.format(i + 1, species.to_chemkin())))
os.path.join(self.output_directory, '{0}_{1}_reaction_sensitivity.png'.format(i + 1, species.to_chemkin())))
if self.thermo_SA:
ThermoSensitivityPlot(x_var=time,
y_var=thermodynamic_sensitivity_data[j * num_ct_species:(j + 1) * num_ct_species]*4184000,
xlabel='dln(c)/d(H_i) [(kcal/mol)^-1]',
num_species=top_sensitive_species).barplot(
os.path.join(self.output_directory, '{0}_{1}_thermo_sensitivity.png'.format(i + 1, species.to_chemkin())))

def simulate(self):
"""
Expand Down Expand Up @@ -407,13 +417,18 @@ def simulate(self):
cantera_simulation = ct.ReactorNet([cantera_reactor])

num_ct_reactions = len(self.model.reactions())
num_ct_species = len(self.model.species())
if self.sensitive_species:
if ct.__version__ == '2.2.1':
print('Warning: Cantera version 2.2.1 may not support sensitivity analysis unless SUNDIALS was used during compilation.')
print('Warning: Upgrade to newer of Cantera in anaconda using the command "conda update -c rmg cantera"')
# Add all the reactions as part of the analysis
for i in range(num_ct_reactions):
cantera_reactor.add_sensitivity_reaction(i)
# if thermo SA is requested, add all species enthalpies
if self.thermo_SA:
for i in range(num_ct_species):
cantera_reactor.add_sensitivity_species_enthalpy(i)
# Set the tolerances for the sensitivity coefficients
cantera_simulation.rtol_sensitivity = 1e-4
cantera_simulation.atol_sensitivity = 1e-6
Expand All @@ -423,7 +438,8 @@ def simulate(self):
temperature = []
pressure = []
species_data = []
sensitivity_data = []
kinetic_sensitivity_data = []
thermo_sensitivity_data = []

# Begin integration
time = 0.0
Expand Down Expand Up @@ -458,21 +474,40 @@ def simulate(self):
for i in range(len(mass_frac_sensitivity_array)):
mass_frac_sensitivity_array[i] *= species_data[-1][i]

# extract kinetics SA
kinetics_mass_frac_sa = mass_frac_sensitivity_array[:, 0:num_ct_reactions]
sensitivity_array = np.zeros(len(self.sensitive_species) * len(self.model.reactions()))
for index, species in enumerate(self.sensitive_species):
for j in range(num_ct_reactions):
sensitivity_array[num_ct_reactions * index + j] = cantera_simulation.sensitivity(
species.to_chemkin(), j)

for i in range(len(mass_frac_sensitivity_array)):
for i in range(len(kinetics_mass_frac_sa)):
if i not in inert_index_list:
# massFracSensitivity for inerts are returned as nan in Cantera, so we must not include them here
sensitivity_array[num_ct_reactions * index + j] -= mass_frac_sensitivity_array[i][j]
sensitivity_data.append(sensitivity_array)

# Convert species_data and sensitivity_data to a numpy array
# massFracSensitivity for inerts are returned as 0.0 in Cantera, so we do not include them here
sensitivity_array[num_ct_reactions * index + j] -= kinetics_mass_frac_sa[i][j]
kinetic_sensitivity_data.append(sensitivity_array)

# extract thermo SA if requested
if self.thermo_SA:
# extract thermo SA
thermo_mass_frac_sa = mass_frac_sensitivity_array[:, num_ct_reactions:]
sensitivity_array = np.zeros(len(self.sensitive_species) * num_ct_species)
for index, species in enumerate(self.sensitive_species):
for j in range(num_ct_species):
sensitivity_array[num_ct_species * index + j] = cantera_simulation.sensitivity(
species.to_chemkin(), j + num_ct_reactions)

for i in range(len(mass_frac_sensitivity_array)):
if i not in inert_index_list:
# massFracSensitivity for inerts are returned as 0.0 in Cantera, so we must not include them here
sensitivity_array[num_ct_species * index + j] -= thermo_mass_frac_sa[i][j]
thermo_sensitivity_data.append(sensitivity_array)

# Convert species_data and sensitivity data to numpy arrays
species_data = np.array(species_data)
sensitivity_data = np.array(sensitivity_data)
kinetic_sensitivity_data = np.array(kinetic_sensitivity_data)
thermo_sensitivity_data = np.array(thermo_sensitivity_data)

# Resave data into generic data objects
time = GenericData(label='Time',
Expand All @@ -497,19 +532,33 @@ def simulate(self):
)
condition_data.append(species_generic_data)

# save kinetic data as generic data object
reaction_sensitivity_data = []
for index, species in enumerate(self.sensitive_species):
for j in range(num_ct_reactions):
reaction_sensitivity_generic_data = GenericData(
label='dln[{0}]/dln[k{1}]: {2}'.format(species.to_chemkin(), j + 1, self.model.reactions()[j]),
species=species,
reaction=self.model.reactions()[j],
data=sensitivity_data[:, num_ct_reactions * index + j],
data=kinetic_sensitivity_data[:, num_ct_reactions * index + j],
index=j + 1,
)
reaction_sensitivity_data.append(reaction_sensitivity_generic_data)

all_data.append((time, condition_data, reaction_sensitivity_data))
# save thermo data as generic data object
thermodynamic_sensitivity_data = []
if self.thermo_SA:
for index, species in enumerate(self.sensitive_species):
for j in range(num_ct_species):
thermo_sensitivity_generic_data = GenericData(
label='dln[{0}]/dH[{1}]'.format(species, self.model.species()[j].name),
species=species,
data=thermo_sensitivity_data[:, num_ct_species * index + j],
index=j + 1,
)
thermodynamic_sensitivity_data.append(thermo_sensitivity_generic_data)

all_data.append((time, condition_data, reaction_sensitivity_data, thermodynamic_sensitivity_data))

return all_data

Expand Down
4 changes: 2 additions & 2 deletions rmgpy/tools/observablesregression.py
View file
Expand Up @@ -235,8 +235,8 @@ def compare(self, tol, plot=False):
fail_header = '\nThe following observables did not match:\n'
fail_header_printed = False
for i in range(len(old_condition_data)):
time_old, data_list_old, reaction_sensitivity_data_old = old_condition_data[i]
time_new, data_list_new, reaction_sensitivity_data_new = new_condition_data[i]
time_old, data_list_old, reaction_sensitivity_data_old, thermodynamic_sensitivity_data_old = old_condition_data[i]
time_new, data_list_new, reaction_sensitivity_data_new, thermodynamic_sensitivity_data_new = new_condition_data[i]

# Compare species observables
if 'species' in self.observables:
Expand Down

0 comments on commit 3c13a22

Please sign in to comment.