Add surface mechanisms to regression tool

Author: sevyharris

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,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.
         ======================= ====================================================
@@ -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 = []