In [ ]:

    

import os

from IPython.display import display
from rmgpy.species import Species
from collections import OrderedDict
import cantera as ct
import numpy as np
import pandas as pd

%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt


    

In [ ]:

    

def simulate(model, T, P, mol_frac, time):

    model.TPX = T, P, mol_frac

    reactor = ct.IdealGasReactor(model)
    
    simulation = ct.ReactorNet([reactor])

    simulation.advance(time)

    conditions = pd.Series()
    conditions['time'] = simulation.time
    conditions['temperature0'] = T
    conditions['temperature'] = model.T
    conditions['pressure0'] = P
    conditions['pressure'] = model.P
        
    species = pd.Series()
    for key, value in model.mole_fraction_dict().iteritems():
        species[key] = value
        
    return conditions, species


    

In [ ]:

    

# Convert chemkin to cantera
from cantera import ck2cti
directory = '/home/mjliu/Documents/models/faa_v2/pah/run4'
filename = 'chem_annotated'
chemkin_file = os.path.join(directory, filename + '.inp')
cantera_file = os.path.join(directory, filename + '.cti')
if os.path.exists(cantera_file):
    raise Exception('File already exists')
ck2cti.Parser().convertMech(chemkin_file, outName=cantera_file)


    

In [ ]:

    

model = ct.Solution(cantera_file)


    

In [ ]:

    

conditions = [800, 1000, 1200, 1400, 1600, 1800, 2000]


    

In [ ]:

    

conditions_df = pd.DataFrame()
species_df = pd.DataFrame()

for T in conditions:
    P = 2 * 101325
    mol_frac = {
        'A1(1)': 0.08,
        'A1C(3)': .01,
        'A2(2)': .01,
        'C2H2(12)': 0.2,
        'N2': 0.7,
    }

    conditions, species = simulate(model, T, P, mol_frac, 0.002)
    
    conditions_df = conditions_df.append(conditions, ignore_index=True)
    species_df = species_df.append(species, ignore_index=True)
    
    print 'T = {0:4} Completed!'.format(T)

data = pd.concat([conditions_df, species_df], axis=1)
print 'Simulation Completed!'


    

In [ ]:

    

plt.style.use('seaborn-poster')

fig = plt.figure()

colormap = mpl.cm.tab20
i = 0
for label in species_df.columns:
    x = data['temperature0']
    y = data[label]
    if any(y > 1e-4) and label not in ['N2', 'C2H2(12)']:
        plt.plot(x, y, c=colormap(i), label=label)
        i += 1

plt.yscale('log')
plt.xlabel('Temperature (K)')
plt.ylabel('Mole Fraction')
plt.legend(bbox_to_anchor=(1.02, 1), loc=2, borderaxespad=0.)


    

In [ ]:

    

# slightly modified version of rmgpy.chemkin.loadSpeciesDictionary
def loadSpeciesDictionary(path):
    speciesDict = OrderedDict()
    with open(path, 'r') as f:
        adjlist = ''
        for line in f:
            if line.strip() == '' and adjlist.strip() != '':
                # Finish this adjacency list
                species = Species().fromAdjacencyList(adjlist)
                label = species.label
                speciesDict[label] = species
                adjlist = ''
            else:
                if "InChI" in line:
                    line = line.split()[0] + '\n'
                if '//' in line:
                    index = line.index('//')
                    line = line[0:index]
                adjlist += line
        else: #reach end of file
            if adjlist.strip() != '':
                species = Species().fromAdjacencyList(adjlist)
                label = species.label
                speciesDict[label] = species

    return speciesDict


    

In [ ]:

    

# set path here
path = os.path.join(directory, 'species_dictionary.txt')

# load
species_dict = loadSpeciesDictionary(path)


    

In [ ]:

    

for label in species_df.columns:
    y = data[label]
    if any(y > 1e-4):
        print label
        display(species_dict[label])


    

In [ ]:

    

for label in species_df.columns:
    y = data[label]
    if all(y < 1e-8) and any(y > 1e-10):
        display(species_dict[label])


    

In [ ]: