In [ ]:
import os
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
from rmgpy.species import Species
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 [ ]:
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)
#species.generateResonanceIsomers()
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)
#species.generateResonanceIsomers()
label = species.label
speciesDict[label] = species
return speciesDict
In [ ]:
directory = '/Users/mjliu/Dropbox (MIT)/Research/Models/Naphthalene/RMG_sim/acetylene/run10'
filename = 'chem_mod'
chemkin_file = os.path.join(directory, filename + '.inp')
cantera_file = os.path.join(directory, filename + '.cti')
In [ ]:
species_dict = loadSpeciesDictionary(os.path.join(directory, 'species_dictionary.txt'))
In [ ]:
# Convert chemkin to cantera
from cantera import ck2cti
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_df = pd.DataFrame()
species_df = pd.DataFrame()
T = 1200
P = 0.5 * 101325
mol_frac = {'C2H2(2)':0.8, 'A2I1':0.2}
conditions, species = simulate(model, T, P, mol_frac, 0.0001)
In [ ]:
weights = []
for label, _ in species.iteritems():
try:
weights.append(1000*species_dict[label].molecule[0].getMolecularWeight())
except KeyError:
weights.append(0)
print label
In [ ]:
species.sort_values(ascending=False)[0:10]
In [ ]:
plt.style.use('seaborn-poster')
#plt.rcParams['axes.labelsize'] = 24
#plt.rcParams['xtick.labelsize'] = 20
#plt.rcParams['ytick.labelsize'] = 20
fig = plt.figure()
plt.stem(weights, species, markerfmt=' ')
plt.xlabel('Mass')
plt.ylabel('Abundance')
In [ ]: