Now saving core species concentration profiles for each reaction system.

We again use the xlwt package to save the concentration profiles as
Excel spreadsheets, which lets us apply a bit of formatting to the data.
Each round of simulation gets its own spreadsheet; these are placed in
the new solver subdirectory. Each reaction system gets a sheet in the
spreadsheet for that iteration. Some information about the reaction
system is printed in the top few lines of each corresponding sheet.

I tried to design it so that RMG would carry on if xlwt was not available.
In the future we may want to make xlwt a required package.

Author: jwallen

rmg.py

@@ -242,11 +242,20 @@ def execute(args):
     except ImportError:
         logging.info('Optional package dependency "psutil" not found; memory profiling information will not be saved.')
 
+    # See if spreadsheet writing package is available
+    saveConcentrationProfiles = False
+    try:
+        import xlwt
+        saveConcentrationProfiles = True
+    except ImportError:
+        logging.info('Optional package dependency "xlwt" not found; reaction system concentration profiles will not be saved.')
+
     # Make output subdirectories
     makeOutputSubdirectory('plot')
     makeOutputSubdirectory('species')
     makeOutputSubdirectory('pdep')
     makeOutputSubdirectory('chemkin')
+    makeOutputSubdirectory('solver')
 
     # Read input file
     reactionModel, coreSpecies, reactionSystems, database, seedMechanisms = readInputFile(inputFile)
@@ -291,11 +300,19 @@ def execute(args):
     done = False
     while not done:
 
+        if saveConcentrationProfiles:
+            workbook = xlwt.Workbook()
+            
         done = True
         objectsToEnlarge = []
         allTerminated = True
         for index, reactionSystem in enumerate(reactionSystems):
 
+            if saveConcentrationProfiles:
+                worksheet = workbook.add_sheet('#{0:d}'.format(index+1))
+            else:
+                worksheet = None
+            
             # Conduct simulation
             logging.info('Conducting simulation of reaction system %s...' % (index+1))
             terminated, obj = reactionSystem.simulate(
@@ -308,6 +325,7 @@ def execute(args):
                 toleranceInterruptSimulation = reactionModel.fluxToleranceInterrupt,
                 termination = reactionModel.termination,
                 pdepNetworks = reactionModel.unirxnNetworks,
+                worksheet = worksheet,
             )
             allTerminated = allTerminated and terminated
             logging.info('')
@@ -323,6 +341,9 @@ def execute(args):
                 objectsToEnlarge.append(obj)
                 done = False
 
+        if saveConcentrationProfiles:
+            workbook.save(os.path.join(settings.outputDirectory, 'solver', 'simulation_{0:d}.xls'.format(len(reactionModel.core.species))))
+
         if not done:
 
             # If we reached our termination conditions, then try to prune

rmgpy/solver/base.pxd

@@ -32,6 +32,7 @@ from pydas cimport DASSL
 
 cdef class ReactionSystem(DASSL):
 
+    cdef public numpy.ndarray coreSpeciesConcentrations
     cdef public numpy.ndarray coreSpeciesRates
     cdef public numpy.ndarray coreReactionRates
     cdef public numpy.ndarray edgeSpeciesRates
@@ -44,9 +45,11 @@ cdef class ReactionSystem(DASSL):
 
     cpdef initializeModel(self, list coreSpecies, list coreReactions, list edgeSpecies, list edgeReactions, list pdepNetworks=?)
 
+    cpdef writeWorksheetHeader(self, worksheet)
+    
     cpdef simulate(self, list coreSpecies, list coreReactions, list edgeSpecies, list edgeReactions,
         double toleranceKeepInEdge, double toleranceMoveToCore, double toleranceInterruptSimulation,
-        list termination, list pdepNetworks=?)
+        list termination, list pdepNetworks=?, worksheet=?)
 
     cpdef logRates(self, double charRate, object species, double speciesRate, object network, double networkRate)
 

rmgpy/solver/base.pyx

@@ -50,6 +50,7 @@ cdef class ReactionSystem(DASSL):
     def __init__(self):
         DASSL.__init__(self)
         # The reaction and species rates at the current time (in mol/m^3*s)
+        self.coreSpeciesConcentrations = None
         self.coreSpeciesRates = None
         self.coreReactionRates = None
         self.edgeSpeciesRates = None
@@ -76,6 +77,7 @@ cdef class ReactionSystem(DASSL):
         numEdgeReactions = len(edgeReactions)
         numPdepNetworks = len(pdepNetworks)
 
+        self.coreSpeciesConcentrations = numpy.zeros((numCoreSpecies), numpy.float64)
         self.coreReactionRates = numpy.zeros((numCoreReactions), numpy.float64)
         self.edgeReactionRates = numpy.zeros((numEdgeReactions), numpy.float64)
         self.coreSpeciesRates = numpy.zeros((numCoreSpecies), numpy.float64)
@@ -85,10 +87,20 @@ cdef class ReactionSystem(DASSL):
         self.maxEdgeSpeciesRates = numpy.zeros((numEdgeSpecies), numpy.float64)
         self.maxNetworkLeakRates = numpy.zeros((numPdepNetworks), numpy.float64)
 
+    
+    cpdef writeWorksheetHeader(self, worksheet):
+        """
+        Write some descriptive information about the reaction system to the
+        first two rows of the given `worksheet`.
+        """
+        import xlwt
+        style0 = xlwt.easyxf('font: bold on')
+        worksheet.write(0, 0, 'Reaction System', style0)
+
     @cython.boundscheck(False)
     cpdef simulate(self, list coreSpecies, list coreReactions, list edgeSpecies, list edgeReactions,
         double toleranceKeepInEdge, double toleranceMoveToCore, double toleranceInterruptSimulation,
-        list termination, list pdepNetworks=None):
+        list termination, list pdepNetworks=None, worksheet=None):
         """
         Simulate the reaction system with the provided reaction model,
         consisting of lists of core species, core reactions, edge species, and
@@ -109,6 +121,7 @@ cdef class ReactionSystem(DASSL):
         cdef numpy.ndarray[numpy.float64_t, ndim=1] maxCoreSpeciesRates, maxEdgeSpeciesRates, maxNetworkLeakRates
         cdef bint terminated
         cdef object maxSpecies, maxNetwork
+        cdef int iteration, i
 
         pdepNetworks = pdepNetworks or []
 
@@ -130,6 +143,7 @@ cdef class ReactionSystem(DASSL):
         maxNetworkIndex = -1
         maxNetwork = None
         maxNetworkRate = 0.0
+        iteration = 0
 
         maxCoreSpeciesRates = self.maxCoreSpeciesRates
         maxEdgeSpeciesRates = self.maxEdgeSpeciesRates
@@ -138,11 +152,28 @@ cdef class ReactionSystem(DASSL):
         # Copy the initial conditions to use in evaluating conversions
         y0 = self.y.copy()
 
+        if worksheet:
+            import xlwt
+            self.writeWorksheetHeader(worksheet)
+            style0 = xlwt.easyxf('font: bold on')
+            style1 = xlwt.easyxf(num_format_str='0.000E+00')
+            worksheet.write(3, 0, 'Time (s)', style0)
+            worksheet.write(3, 1, 'Concentrations (mol/m^3)', style0)
+            for i in range(numCoreSpecies):
+                worksheet.write(4, i+1, str(coreSpecies[i]), style0)
+            
         stepTime = 1e-12
         while not terminated:
             # Integrate forward in time by one time step
             self.step(stepTime)
 
+            iteration += 1
+
+            if worksheet:
+                worksheet.write(iteration+4, 0, self.t, style1)
+                for i in range(numCoreSpecies):
+                    worksheet.write(iteration+4, i+1, self.coreSpeciesConcentrations[i], style1)
+
             # Get the characteristic flux
             charRate = sqrt(numpy.sum(self.coreSpeciesRates * self.coreSpeciesRates))
 

rmgpy/solver/simple.pyx

@@ -144,11 +144,22 @@ cdef class SimpleReactor(ReactionSystem):
         y0 = numpy.zeros((numCoreSpecies), numpy.float64)
         for spec, moleFrac in self.initialMoleFractions.iteritems():
             y0[speciesIndex[spec]] = moleFrac * (self.P / constants.R / self.T)
+            self.coreSpeciesConcentrations[speciesIndex[spec]] = y0[speciesIndex[spec]]
 
         # Initialize the model
         dydt0 = - self.residual(t0, y0, numpy.zeros((numCoreSpecies), numpy.float64))[0]
         DASSL.initialize(self, t0, y0, dydt0)
 
+    cpdef writeWorksheetHeader(self, worksheet):
+        """
+        Write some descriptive information about the reaction system to the
+        first two rows of the given `worksheet`.
+        """
+        import xlwt
+        style0 = xlwt.easyxf('font: bold on')
+        worksheet.write(0, 0, 'Simple Reactor', style0)
+        worksheet.write(1, 0, 'T = {0:g} K, P = {1:g} bar'.format(self.T, self.P/1e5))
+
     @cython.boundscheck(False)
     def residual(self, double t, numpy.ndarray[numpy.float64_t, ndim=1] y, numpy.ndarray[numpy.float64_t, ndim=1] dydt):
 
@@ -178,12 +189,16 @@ cdef class SimpleReactor(ReactionSystem):
         numEdgeReactions = len(self.edgeReactionRates)
         numPdepNetworks = len(self.networkLeakRates)
 
+        coreSpeciesConcentrations = numpy.zeros_like(self.coreSpeciesConcentrations)
         coreSpeciesRates = numpy.zeros_like(self.coreSpeciesRates)
         coreReactionRates = numpy.zeros_like(self.coreReactionRates)
         edgeSpeciesRates = numpy.zeros_like(self.edgeSpeciesRates)
         edgeReactionRates = numpy.zeros_like(self.edgeReactionRates)
         networkLeakRates = numpy.zeros_like(self.networkLeakRates)
 
+        for j in range(y.shape[0]):
+            coreSpeciesConcentrations[j] = y[j]
+        
         for j in range(ir.shape[0]):
             k = kf[j]
             if ir[j,0] >= numCoreSpecies or ir[j,1] >= numCoreSpecies or ir[j,2] >= numCoreSpecies:
@@ -264,6 +279,7 @@ cdef class SimpleReactor(ReactionSystem):
                 reactionRate = k * y[inet[j,0]] * y[inet[j,1]] * y[inet[j,2]]
             networkLeakRates[j] = reactionRate
 
+        self.coreSpeciesConcentrations = coreSpeciesConcentrations
         self.coreSpeciesRates = coreSpeciesRates
         self.coreReactionRates = coreReactionRates
         self.edgeSpeciesRates = edgeSpeciesRates