Skip to content

Add surface mechanisms to regression tool #2563

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
sevyharris wants to merge 1 commit into from
Closed

Add surface mechanisms to regression tool #2563

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
96 changes: 80 additions & 16 deletions rmgpy/tools/canteramodel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ class CanteraCondition(object):

"""

def __init__(self, reactor_type, reaction_time, mol_frac, T0=None, P0=None, V0=None):
def __init__(self, reactor_type, reaction_time, mol_frac, surface_mol_frac=None, T0=None, P0=None, V0=None):
self.reactor_type = reactor_type
self.reaction_time = Quantity(reaction_time)

Expand All @@ -75,6 +75,14 @@ def __init__(self, reactor_type, reaction_time, mol_frac, T0=None, P0=None, V0=N

self.mol_frac = mol_frac

if surface_mol_frac:
# normalize surface mol fractions
if sum(surface_mol_frac.values()) != 1.00:
total = sum(surface_mol_frac.values())
for species, value in surface_mol_frac.items():
surface_mol_frac[species] = value / total
self.surface_mol_frac = surface_mol_frac

# Check to see that one of the three attributes T0, P0, and V0 is less unspecified
props = [T0, P0, V0]
total = 0
Expand Down Expand Up @@ -121,7 +129,7 @@ def __str__(self):
return string


def generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_list, Tlist=None, Plist=None, Vlist=None):
def generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_list, surface_mol_frac_list=[None], Tlist=None, Plist=None, Vlist=None):
"""
Creates a list of cantera conditions from from the arguments provided.

Expand All @@ -137,6 +145,8 @@ def generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_
`reaction_time_list` A tuple object giving the ([list of reaction times], units)
`mol_frac_list` A list of molfrac dictionaries with species object keys
and mole fraction values
`surface_mol_frac_list` A list of molfrac dictionaries with surface species object keys
and mole fraction values
To specify the system for an ideal gas, you must define 2 of the following 3 parameters:
`T0List` A tuple giving the ([list of initial temperatures], units)
'P0List' A tuple giving the ([list of initial pressures], units)
Expand Down Expand Up @@ -167,25 +177,28 @@ def generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_
for reactor_type in reactor_type_list:
for reaction_time in reaction_time_list:
for mol_frac in mol_frac_list:
for P in Plist:
for V in Vlist:
conditions.append(CanteraCondition(reactor_type, reaction_time, mol_frac, P0=P, V0=V))
for surface_mol_frac in surface_mol_frac_list:
for P in Plist:
for V in Vlist:
conditions.append(CanteraCondition(reactor_type, reaction_time, mol_frac, surface_mol_frac=surface_mol_frac, P0=P, V0=V))

elif Plist is None:
for reactor_type in reactor_type_list:
for reaction_time in reaction_time_list:
for mol_frac in mol_frac_list:
for T in Tlist:
for V in Vlist:
conditions.append(CanteraCondition(reactor_type, reaction_time, mol_frac, T0=T, V0=V))
for surface_mol_frac in surface_mol_frac_list:
for T in Tlist:
for V in Vlist:
conditions.append(CanteraCondition(reactor_type, reaction_time, mol_frac, surface_mol_frac=surface_mol_frac, T0=T, V0=V))

elif Vlist is None:
for reactor_type in reactor_type_list:
for reaction_time in reaction_time_list:
for mol_frac in mol_frac_list:
for T in Tlist:
for P in Plist:
conditions.append(CanteraCondition(reactor_type, reaction_time, mol_frac, T0=T, P0=P))
for surface_mol_frac in surface_mol_frac_list:
for T in Tlist:
for P in Plist:
conditions.append(CanteraCondition(reactor_type, reaction_time, mol_frac, surface_mol_frac=surface_mol_frac, T0=T, P0=P))

else:
raise Exception("Cantera conditions must leave one of T0, P0, and V0 state variables unspecified")
Expand All @@ -198,7 +211,7 @@ class Cantera(object):
"""

def __init__(self, species_list=None, reaction_list=None, canteraFile='', output_directory='', conditions=None,
sensitive_species=None, thermo_SA=False):
sensitive_species=None, thermo_SA=False, surface=False, surface_species_list=None, surface_reaction_list=None):
"""
`species_list`: list of RMG species objects
`reaction_list`: list of RMG reaction objects
Expand All @@ -208,24 +221,39 @@ def __init__(self, species_list=None, reaction_list=None, canteraFile='', output
`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.
`surface`: a boolean indicating whether or not to run a surface simulation
"""
self.species_list = species_list
self.reaction_list = reaction_list
self.reaction_map = {}
self.model = ct.Solution(canteraFile) if canteraFile else None
self.surface = surface
self.surface_species_list = surface_species_list
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

# need a way to make a surface simulation... for now, we'll assume there's a Cantera surface yaml file
if surface:
# assert canteraFile, "Cantera file must be specified to run a surface simulation" # or we add the capability to load it later
assert surface_species_list, "Surface species list must be specified to run a surface simulation"
# TODO some kind of assertion for starting surface mol fractions
# self.surface_species_list = surface_species_list

if canteraFile:
self.surface = ct.Interface(canteraFile, 'surface1')
self.model = self.surface.adjacent['gas']
# also need starting surface mol fractions

# Make output directory if it does not yet exist:
if not os.path.exists(self.output_directory):
try:
os.makedirs(self.output_directory)
except:
raise Exception('Cantera output directory could not be created.')

def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_list, Tlist=None, Plist=None, Vlist=None):
def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_list, surface_mol_frac_list=[None], Tlist=None, Plist=None, Vlist=None):
"""
This saves all the reaction conditions into the Cantera class.
======================= ====================================================
Expand All @@ -240,19 +268,21 @@ def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_li
`reaction_time_list` A tuple object giving the ([list of reaction times], units)
`mol_frac_list` A list of molfrac dictionaries with species object keys
and mole fraction values
`surface_mol_frac_list` A list of molfrac dictionaries with surface species object keys and mole fraction values
To specify the system for an ideal gas, you must define 2 of the following 3 parameters:
`T0List` A tuple giving the ([list of initial temperatures], units)
'P0List' A tuple giving the ([list of initial pressures], units)
'V0List' A tuple giving the ([list of initial specific volumes], units)
"""

self.conditions = generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_list, Tlist, Plist, Vlist)
self.conditions = generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_list, surface_mol_frac_list, Tlist, Plist, Vlist)

def load_model(self):
"""
Load a cantera Solution model from the job's own species_list and reaction_list attributes
"""

# TODO add surface version - this is limited by the rmgpy.kinetics.arrhenius.Arrhenius.to_cantera_kinetics()
# function not being able to handle surface reactions yet
ct_species = [spec.to_cantera(use_chemkin_identifier=True) for spec in self.species_list]

self.reaction_map = {}
Expand Down Expand Up @@ -292,6 +322,7 @@ def load_chemkin_model(self, chemkin_file, transport_file=None, **kwargs):
Convert a chemkin mechanism chem.inp file to a cantera mechanism file chem.yaml
and save it in the output_directory
Then load it into self.model
Note that if this is a surface mechanism, the calling function much include surface_file as a keyword argument
"""
from cantera import ck2yaml

Expand All @@ -301,9 +332,14 @@ def load_chemkin_model(self, chemkin_file, transport_file=None, **kwargs):
if os.path.exists(out_name):
os.remove(out_name)
parser = ck2yaml.Parser()

parser.convert_mech(chemkin_file, transport_file=transport_file, out_name=out_name, **kwargs)
self.model = ct.Solution(out_name)

if self.surface:
self.surface = ct.Interface(out_name, 'surface1')
self.model = self.surface.adjacent['gas']

def modify_reaction_kinetics(self, rmg_reaction_index, rmg_reaction):
"""
Modify the corresponding cantera reaction's kinetics to match
Expand Down Expand Up @@ -385,6 +421,9 @@ def simulate(self):
species_names_list = [get_species_identifier(species) for species in self.species_list]
inert_index_list = [self.species_list.index(species) for species in self.species_list if species.index == -1]

if self.surface:
surface_species_names_list = [get_species_identifier(species) for species in self.surface_species_list]

all_data = []
for condition in self.conditions:

Expand All @@ -394,6 +433,12 @@ def simulate(self):
newkey = get_species_identifier(key)
new_mol_frac[newkey] = value

if self.surface:
new_surface_mol_frac = {}
for key, value in condition.surface_mol_frac.items(): # this needs to be aspecies object
newkey = get_species_identifier(key)
new_surface_mol_frac[newkey] = value

# Set Cantera simulation conditions
if condition.V0 is None:
self.model.TPX = condition.T0.value_si, condition.P0.value_si, new_mol_frac
Expand All @@ -413,6 +458,11 @@ def simulate(self):
else:
raise Exception('Other types of reactor conditions are currently not supported')

# add the surface/gas adjacent thing...
if self.surface:
rsurf = ct.ReactorSurface(self.surface, cantera_reactor)
self.surface.TPX = condition.T0.value_si, condition.P0.value_si, new_surface_mol_frac

# Run this individual condition as a simulation
cantera_simulation = ct.ReactorNet([cantera_reactor])

Expand Down Expand Up @@ -451,7 +501,12 @@ def simulate(self):
times.append(cantera_simulation.time)
temperature.append(cantera_reactor.T)
pressure.append(cantera_reactor.thermo.P)
species_data.append(cantera_reactor.thermo[species_names_list].X)

if self.surface:
species_data.append(np.concatenate((cantera_reactor.thermo[species_names_list].X, rsurf.kinetics.coverages)))
N_gas = len(cantera_reactor.thermo[species_names_list].X)
else:
species_data.append(cantera_reactor.thermo[species_names_list].X)

if self.sensitive_species:
# Cantera returns mass-based sensitivities rather than molar concentration or mole fraction based sensitivities.
Expand Down Expand Up @@ -531,6 +586,15 @@ def simulate(self):
index=species.index
)
condition_data.append(species_generic_data)
if self.surface:
for index, species in enumerate(self.surface_species_list):
# Create generic data object that saves the surface species object into the species object.
species_generic_data = GenericData(label=surface_species_names_list[index],
species=species,
data=species_data[:, index + N_gas],
index=species.index + N_gas
)
condition_data.append(species_generic_data)

# save kinetic data as generic data object
reaction_sensitivity_data = []
Expand Down
Loading