Skip to content

Commit

Permalink
Merge pull request #1669 from ReactionMechanismGenerator/simulate_MB_…
Browse files Browse the repository at this point in the history
…sampled_reactor

Created a new reactor type MBSampledReactor for modeling laser lab Mass-Spec experiments
  • Loading branch information
mliu49 authored Aug 5, 2019
2 parents d66fbfe + 206f149 commit f0131fe
Show file tree
Hide file tree
Showing 9 changed files with 1,400 additions and 6 deletions.
12 changes: 10 additions & 2 deletions rmgpy/chemkin.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -1411,7 +1411,11 @@ def getSpeciesIdentifier(species):
# Try the chemical formula
name = '{0}({1:d})'.format(species.molecule[0].getFormula(), species.index)
if len(name) <= 10:
return name
if 'obs' in label:
# For MBSampledReactor, keep observed species tag
return name + '_obs'
else:
return name

# As a last resort, just use the index
if species.index >= 0:
Expand All @@ -1421,7 +1425,11 @@ def getSpeciesIdentifier(species):
else:
name = 'S({0:d})'.format(species.index)
if len(name) <= 10:
return name
if 'obs' in label:
# For MBSampledReactor, keep observed species tag
return name + '_obs'
else:
return name

# If we're here then we just can't come up with a valid Chemkin name
# for this species, so raise an exception
Expand Down
67 changes: 66 additions & 1 deletion rmgpy/rmg/input.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,10 +38,11 @@
from rmgpy import settings

from rmgpy.molecule import Molecule
from rmgpy.quantity import Quantity, Energy, SurfaceConcentration
from rmgpy.quantity import Quantity, Energy, RateCoefficient, SurfaceConcentration
from rmgpy.solver.base import TerminationTime, TerminationConversion, TerminationRateRatio
from rmgpy.solver.simple import SimpleReactor
from rmgpy.solver.liquid import LiquidReactor
from rmgpy.solver.mbSampled import MBSampledReactor
from rmgpy.solver.surface import SurfaceReactor
from rmgpy.rmg.settings import ModelSettings, SimulatorSettings
from model import CoreEdgeReactionModel
Expand Down Expand Up @@ -447,6 +448,69 @@ def surfaceReactor(temperature,
rmg.reactionSystems.append(system)
system.log_initial_conditions(number=len(rmg.reactionSystems))


# Reaction systems
def mbsampledReactor(temperature,
pressure,
initialMoleFractions,
mbsamplingRate,
terminationConversion=None,
terminationTime=None,
sensitivity=None,
sensitivityThreshold=1e-3,
constantSpecies=None,
):
logging.debug('Found MBSampledReactor reaction system')

for value in initialMoleFractions.values():
if value < 0:
raise InputError('Initial mole fractions cannot be negative.')

for spec in initialMoleFractions:
initialMoleFractions[spec] = float(initialMoleFractions[spec])

totalInitialMoles = sum(initialMoleFractions.values())
if totalInitialMoles != 1:
logging.warning('Initial mole fractions do not sum to one; normalizing.')
logging.info('')
logging.info('Original composition:')
for spec, molfrac in initialMoleFractions.iteritems():
logging.info("{0} = {1}".format(spec, molfrac))
for spec in initialMoleFractions:
initialMoleFractions[spec] /= totalInitialMoles
logging.info('')
logging.info('Normalized mole fractions:')
for spec, molfrac in initialMoleFractions.iteritems():
logging.info("{0} = {1}".format(spec, molfrac))

T = Quantity(temperature)
P = Quantity(pressure)

k_sampling = RateCoefficient(mbsamplingRate, 's^-1')

constantSpeciesList = []

for spec in constantSpecies:
constantSpeciesList.append(speciesDict[spec])

termination = []
if terminationConversion is not None:
for spec, conv in terminationConversion.iteritems():
termination.append(TerminationConversion(speciesDict[spec], conv))
if terminationTime is not None:
termination.append(TerminationTime(Quantity(terminationTime)))
if len(termination) == 0:
raise InputError(
'No termination conditions specified for reaction system #{0}.'.format(len(rmg.reactionSystems) + 2))

sensitiveSpecies = []
if sensitivity:
if isinstance(sensitivity, str): sensitivity = [sensitivity]
for spec in sensitivity:
sensitiveSpecies.append(speciesDict[spec])
system = MBSampledReactor(T, P, initialMoleFractions, k_sampling, constantSpeciesList, termination, sensitiveSpecies, sensitivityThreshold)
rmg.reactionSystems.append(system)

def simulator(atol, rtol, sens_atol=1e-6, sens_rtol=1e-4):
rmg.simulatorSettingsList.append(SimulatorSettings(atol, rtol, sens_atol, sens_rtol))

Expand Down Expand Up @@ -756,6 +820,7 @@ def readInputFile(path, rmg0):
'simpleReactor': simpleReactor,
'liquidReactor': liquidReactor,
'surfaceReactor': surfaceReactor,
'mbsampledReactor': mbsampledReactor,
'simulator': simulator,
'solvation': solvation,
'model': model,
Expand Down
Loading

0 comments on commit f0131fe

Please sign in to comment.