diff --git a/rmgpy/tools/canteramodel.py b/rmgpy/tools/canteramodel.py index 99aed79de1..a23e14705d 100644 --- a/rmgpy/tools/canteramodel.py +++ b/rmgpy/tools/canteramodel.py @@ -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) @@ -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 @@ -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. @@ -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) @@ -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") @@ -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 @@ -208,16 +221,31 @@ 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: @@ -225,7 +253,7 @@ def __init__(self, species_list=None, reaction_list=None, canteraFile='', output 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. ======================= ==================================================== @@ -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 = {} @@ -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 @@ -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 @@ -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: @@ -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 @@ -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]) @@ -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. @@ -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 = [] diff --git a/rmgpy/tools/observablesregression.py b/rmgpy/tools/observablesregression.py index 1cce9e1750..cb966aa432 100644 --- a/rmgpy/tools/observablesregression.py +++ b/rmgpy/tools/observablesregression.py @@ -112,32 +112,74 @@ def __init__(self, title='', old_dir='', new_dir='', observables=None, expt_data if os.path.exists(os.path.join(new_dir, 'tran.dat')): new_transport_path = os.path.join(new_dir, 'tran.dat') - # load the species and reactions from each model - old_species_list, old_reaction_list = load_chemkin_file(os.path.join(old_dir, 'chem_annotated.inp'), - os.path.join(old_dir, 'species_dictionary.txt'), - old_transport_path) - - new_species_list, new_reaction_list = load_chemkin_file(os.path.join(new_dir, 'chem_annotated.inp'), - os.path.join(new_dir, 'species_dictionary.txt'), - new_transport_path) + # define chemkin file path for old and new models + old_chemkin_path = os.path.join(old_dir, 'chem_annotated.inp') + new_chemkin_path = os.path.join(new_dir, 'chem_annotated.inp') + old_species_dict_path = os.path.join(old_dir, 'species_dictionary.txt') + new_species_dict_path = os.path.join(new_dir, 'species_dictionary.txt') + + surface = False + if os.path.exists(os.path.join(old_dir, 'chem_annotated-surface.inp')) and \ + os.path.exists(os.path.join(old_dir, 'chem_annotated-gas.inp')): + surface = True + old_chemkin_path = os.path.join(old_dir, 'chem_annotated-gas.inp') + new_chemkin_path = os.path.join(new_dir, 'chem_annotated-gas.inp') + old_surface_chemkin_path = os.path.join(old_dir, 'chem_annotated-surface.inp') + new_surface_chemkin_path = os.path.join(new_dir, 'chem_annotated-surface.inp') + ck2cti = True # can't create a surface model from memory yet, must load from yaml + + old_species_list, old_reaction_list = load_chemkin_file( + old_chemkin_path, + old_species_dict_path, + old_transport_path + ) + new_species_list, new_reaction_list = load_chemkin_file( + new_chemkin_path, + new_species_dict_path, + new_transport_path + ) + + old_surface_species_list = None + new_surface_species_list = None + new_surface_reaction_list = None + old_surface_reaction_list = None + if surface: + old_surface_species_list, old_surface_reaction_list = load_chemkin_file( + old_surface_chemkin_path, + old_species_dict_path, + old_transport_path + ) + new_surface_species_list, new_surface_reaction_list = load_chemkin_file( + new_surface_chemkin_path, + new_species_dict_path, + new_transport_path + ) self.old_sim = Cantera(species_list=old_species_list, - reaction_list=old_reaction_list, - output_directory=old_dir) + reaction_list=old_reaction_list, + output_directory=old_dir, + surface=surface, + surface_species_list=old_surface_species_list, + ) self.new_sim = Cantera(species_list=new_species_list, - reaction_list=new_reaction_list, - output_directory=new_dir) + reaction_list=new_reaction_list, + output_directory=new_dir, + surface=surface, + surface_species_list=new_surface_species_list, + ) # load each chemkin file into the cantera model - if not ck2cti: + if not ck2cti: # TODO make this work for surfaces self.old_sim.load_model() self.new_sim.load_model() else: - self.old_sim.load_chemkin_model(os.path.join(old_dir, 'chem_annotated.inp'), + self.old_sim.load_chemkin_model(old_chemkin_path, transport_file=old_transport_path, + surface_file=old_surface_chemkin_path, quiet=True) - self.new_sim.load_chemkin_model(os.path.join(new_dir, 'chem_annotated.inp'), + self.new_sim.load_chemkin_model(new_chemkin_path, transport_file=new_transport_path, + surface_file=new_surface_chemkin_path, quiet=True) def __str__(self): @@ -146,7 +188,7 @@ def __str__(self): """ return 'Observables Test Case: {0}'.format(self.title) - def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_list, + def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_list, surface_mol_frac_list=[None], Tlist=None, Plist=None, Vlist=None): """ Creates a list of conditions from from the lists provided. @@ -170,6 +212,7 @@ def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_li """ # Store the conditions in the observables test case, for bookkeeping self.conditions = generate_cantera_conditions(reactor_type_list, reaction_time_list, mol_frac_list, + surface_mol_frac_list=surface_mol_frac_list, Tlist=Tlist, Plist=Plist, Vlist=Vlist) # Map the mole fractions dictionaries to species objects from the old and new models @@ -192,10 +235,14 @@ def generate_conditions(self, reactor_type_list, reaction_time_list, mol_frac_li old_mol_frac_list.append(old_condition) new_mol_frac_list.append(new_condition) + # TODO have to do the same for surface species + # Generate the conditions in each simulation self.old_sim.generate_conditions(reactor_type_list, reaction_time_list, old_mol_frac_list, + surface_mol_frac_list=surface_mol_frac_list, Tlist=Tlist, Plist=Plist, Vlist=Vlist) self.new_sim.generate_conditions(reactor_type_list, reaction_time_list, new_mol_frac_list, + surface_mol_frac_list=surface_mol_frac_list, Tlist=Tlist, Plist=Plist, Vlist=Vlist) def compare(self, tol, plot=False): @@ -205,7 +252,7 @@ def compare(self, tol, plot=False): `plot`: if set to True, it will comparison plots of the two models comparing their species. Returns a list of variables failed in a list of tuples in the format: - + (CanteraCondition, variable label, variable_old, variable_new) """ @@ -220,10 +267,14 @@ def compare(self, tol, plot=False): print('{0} Comparison'.format(self)) # Check the species profile observables if 'species' in self.observables: - old_species_dict = get_rmg_species_from_user_species(self.observables['species'], self.old_sim.species_list) - new_species_dict = get_rmg_species_from_user_species(self.observables['species'], self.new_sim.species_list) - - # Check state variable observables + if self.old_sim.surface_species_list: + old_species_dict = get_rmg_species_from_user_species(self.observables['species'], self.old_sim.species_list + self.old_sim.surface_species_list) + new_species_dict = get_rmg_species_from_user_species(self.observables['species'], self.new_sim.species_list + self.new_sim.surface_species_list) + else: + old_species_dict = get_rmg_species_from_user_species(self.observables['species'], self.old_sim.species_list) + new_species_dict = get_rmg_species_from_user_species(self.observables['species'], self.new_sim.species_list) + + # Check state variable observables implemented_variables = ['temperature', 'pressure'] if 'variable' in self.observables: for item in self.observables['variable']: diff --git a/rmgpy/tools/regression.py b/rmgpy/tools/regression.py index ad2c952592..d3bd94cbba 100644 --- a/rmgpy/tools/regression.py +++ b/rmgpy/tools/regression.py @@ -98,6 +98,15 @@ def read_input_file(path): new_imfs.append(new_dict) setups[3] = new_imfs + # TODO repeat this for surface species, make sure it's happening for them too + if setups[5][0]: + ismfs = setups[5] + new_ismfs = [] + for ismf in ismfs: + new_dict = convert(ismf, initialSpecies) + new_ismfs.append(new_dict) + setups[5] = new_ismfs + return casetitle, observables, setups, tol @@ -112,14 +121,16 @@ def species(label, structure): return spc -def reactorSetups(reactorTypes, temperatures, pressures, initialMoleFractionsList, terminationTimes): +def reactorSetups(reactorTypes, temperatures, pressures, initialMoleFractionsList, terminationTimes, initialSurfaceMoleFractionsList=[None]): global setups + # TODO add initialsufrace + terminationTimes = Quantity(terminationTimes) temperatures = Quantity(temperatures) pressures = Quantity(pressures) - setups = [reactorTypes, temperatures, pressures, initialMoleFractionsList, terminationTimes] + setups = [reactorTypes, temperatures, pressures, initialMoleFractionsList, terminationTimes, initialSurfaceMoleFractionsList] def SMILES(string): @@ -159,11 +170,12 @@ def run(benchmarkDir, testDir, title, observables, setups, tol): ck2cti=False, ) - reactor_types, temperatures, pressures, initial_mole_fractions_list, termination_times = setups + reactor_types, temperatures, pressures, initial_mole_fractions_list, termination_times, initial_surface_mole_fractions_list = setups case.generate_conditions( reactor_type_list=reactor_types, reaction_time_list=termination_times, mol_frac_list=initial_mole_fractions_list, + surface_mol_frac_list=initial_surface_mole_fractions_list, Tlist=temperatures, Plist=pressures )