Correct initial conditions from Cantera solver

Get the y0 vector from the cantera solver, so comparisons of later y
vectors are on an equal footing.

examples/diesel/input.xml

@@ -68,14 +68,14 @@
 	<!-- Termination -->
 	<termination>
 		<target type="conversion" speciesID="oxygen">0.5</target>
-		<target type="time" units="s">1e12</target>
+		<target type="time" units="s">3600</target>
 	</termination>
 	<!-- Dynamic simulator -->
 	<simulator atol="1e-16" rtol="1e-8" />
 	<fluxTolerance>
-		<keepInEdge>1e-5</keepInEdge>
+		<keepInEdge>0.0001</keepInEdge>
 		<moveToCore>0.01</moveToCore>
-		<interruptSimulation>0.1</interruptSimulation>
+		<interruptSimulation>0.5</interruptSimulation>
 	</fluxTolerance>
 	<!-- Options -->
 	<drawMolecules>on</drawMolecules>

source/rmg/model.py

@@ -902,17 +902,16 @@ def simulate(self, model):
 		for target in model.termination:
 			if target.__class__ == TerminationTime:
 				endtime = target.time
-				
+
 		# Set up initial conditions
-		P = float(self.pressureModel.getPressure(0))
-		T = float(self.temperatureModel.getTemperature(0))
-		V = 1.0 # [=] m**3
-		Ni = numpy.zeros(len(model.core.species), float)
-		for i, spec in enumerate(model.core.species):
-			if spec in self.initialConcentration:
-				Ni[i] = self.initialConcentration[spec] * V
-		Ni0 = Ni
+		P = gas.pressure()
+		V = sim.reactors()[0].volume()
+		T = gas.temperature()
+		# recall that Cantera returns molarDensity() in units of kmol/m3
+		# and this program thinks in mol/m3
+		Ni = gas.molarDensity()*1000.0 * gas.moleFractions() * V 
 		y = [P, V, T]; y.extend(Ni)
+		Ni0 = Ni
 		y0 = y
 
 #		# Output information about simulation at current time
@@ -925,12 +924,7 @@ def simulate(self, model):
 #		tlist.append(0.0); ylist.append(y0)
 #		dydtlist.append(self.getResidual(0.0, y0, model, stoichiometry))
 
-		# Set up solver
-		#solver = scipy.integrate.ode(self.getResidual,None)
-		#solver.set_integrator('vode', method='bdf', with_jacobian=True, atol=model.absoluteTolerance, rtol=model.relativeTolerance)
-		#solver.set_f_params(model, stoichiometry)
-		#solver.set_initial_value(y0,0.0)
-		#
+
 		solver = FakeSolver()
 
 		done = False
@@ -992,6 +986,7 @@ def simulate(self, model):
 				#		print model.core.species[i], maxSpeciesFluxes[i]
 				#	else:
 				#		print model.edge.species[i-len(model.core.species)], maxSpeciesFluxes[i]
+				print gas
 				return tlist, ylist, dydtlist, False, maxSpecies
 
 			# Test for simulation completion
@@ -1003,6 +998,8 @@ def simulate(self, model):
 				elif target.__class__ == TerminationTime:
 					if time > target.time: done = True
 
+		print gas 
+
 		# Test for model validity once simulation complete
 		maxSpecies = None
 		maxRelativeFlux = 0.0