From f20d39cf7594cfbe8a5cf8de0548232f0be1f5b8 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 30 Jan 2013 22:05:20 +0000 Subject: [PATCH] [Cython] Implemented class Mixture (including unit tests) --- interfaces/cython/cantera/_cantera.pxd | 37 +- interfaces/cython/cantera/mixture.pyx | 323 +++++++++++++++++- .../cython/cantera/test/test_mixture.py | 176 +++++++++- 3 files changed, 509 insertions(+), 27 deletions(-) diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 5400e15b2..229b8db12 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -191,12 +191,39 @@ cdef extern from "cantera/equil/MultiPhase.h" namespace "Cantera": CxxMultiPhase() void addPhase(CxxThermoPhase*, double) except + void init() except + - double nSpecies() - void setTemperature(double) + + size_t nSpecies() + size_t nElements() + size_t nPhases() + size_t elementIndex(string) except + + size_t speciesIndex(size_t, size_t) except + + string speciesName(size_t) except + + double nAtoms(size_t, size_t) except + + + double phaseMoles(size_t) except + + void setPhaseMoles(size_t, double) except + + void setMoles(double*) except + + void setMolesByName(string) except + + + double speciesMoles(size_t) except + + double elementMoles(size_t) except + + + void setTemperature(double) except + double temperature() - void setPressure(double) + void setPressure(double) except + double pressure() + double minTemp() except + + double maxTemp() except + + double charge() except + + double phaseCharge(size_t) except + + void getChemPotentials(double*) except + + double enthalpy() except + + double entropy() except + + double gibbs() except + + double cp() except + + double volume() except + + cdef extern from "cantera/equil/equil.h" namespace "Cantera": int equilibrate(CxxThermoPhase&, char*, int, double, int, int, int) except + @@ -497,10 +524,6 @@ cdef class _SolutionBase: cdef np.ndarray _selectedSpecies cdef object parent -cdef class Mixture: - cdef CxxMultiPhase* mix - cdef list _phases - cdef class Kinetics(_SolutionBase): pass diff --git a/interfaces/cython/cantera/mixture.pyx b/interfaces/cython/cantera/mixture.pyx index 8dd97362d..79d6d20aa 100644 --- a/interfaces/cython/cantera/mixture.pyx +++ b/interfaces/cython/cantera/mixture.pyx @@ -1,4 +1,36 @@ cdef class Mixture: + """ + + Class Mixture represents mixtures of one or more phases of matter. To + construct a mixture, supply a list of phases to the constructor, each + paired with the number of moles for that phase:: + + >>> gas = cantera.Solution('gas.cti') + >>> gas.speciesNames + ['H2', 'H', 'O2', 'O', 'OH'] + >>> graphite = cantera.Solution('graphite.cti') + >>> graphite.speciesNames + ['C(g)'] + >>> mix = Mixture([(gas, 1.0), (graphite, 0.1)]) + >>> mix.speciesNames + ['H2', 'H', 'O2', 'O', 'OH', 'C(g)'] + + Note that the objects representing each phase compute only the intensive + state of the phase -- they do not store any information on the amount of + this phase. Mixture objects, on the other hand, represent the full + extensive state. + + Mixture objects are 'lightweight' in the sense that they do not store + parameters needed to compute thermodynamic or kinetic properties of the + phases. These are contained in the ('heavyweight') phase objects. Multiple + mixture objects may be constructed using the same set of phase objects. + Each one stores its own state information locally, and synchronizes the + phases objects whenever it requires phase properties. + """ + + cdef CxxMultiPhase* mix + cdef list _phases + def __cinit__(self, phases): self.mix = new CxxMultiPhase() self._phases = [] @@ -10,35 +42,300 @@ cdef class Mixture: self.mix.init() if self._phases: - self.pressure = self._phases[0].P - self.temperature = self._phases[0].T + self.P = self._phases[0].P + self.T = self._phases[0].T def __dealloc__(self): del self.mix - def phase(self, n): - return self._phases[n] + def report(self): + """ + Generate a report describing the thermodynamic state of this mixture. To + print the report to the screen, simply call the mixture object. The + following two statements are equivalent:: - def equilibrate(self, XY, solver=1, int estimateEquil=0, double err=1e-9, - int maxsteps=1000, - int maxiter=200, int loglevel=0, printlevel=0): - XY = XY.upper() - vcs_equilibrate(deref(self.mix), stringify(XY).c_str(), - estimateEquil, printlevel, solver, err, - maxsteps, maxiter, loglevel) + >>> mix() + >>> print(mix.report()) + """ + s = [] + for i,phase in enumerate(self._phases): + s.append('************ Phase {0} ************'.format(phase.name)) + s.append('Moles: {0}'.format(self.phaseMoles(i))) + s.append(phase.report()) + + return '\n'.join(s) + + def __call__(self): + print(self.report()) + + property nElements: + """Total number of elements present in the mixture.""" + def __get__(self): + return self.mix.nElements() + + cpdef int elementIndex(self, element) except *: + """Index of element with name 'element':: + + >>> mix.elementIndex('H') + 2 + """ + if isinstance(element, (str, unicode)): + index = self.mix.elementIndex(stringify(element)) + elif isinstance(element, (int, float)): + index = element + else: + raise TypeError("'element' must be a string or a number") + + if not 0 <= index < self.nElements: + raise ValueError('No such element.') + + return index property nSpecies: + """Number of species.""" def __get__(self): return self.mix.nSpecies() - property temperature: + def speciesName(self, k): + """Name of the species with index *k*. Note that index numbers + are assigned in order as phases are added.""" + return pystr(self.mix.speciesName(k)) + + property speciesNames: + def __get__(self): + return [self.speciesName(k) for k in range(self.nSpecies)] + + def speciesIndex(self, phase, species): + """ + :param phase: + Phase object, index or name + :param species: + Species name or index + + Returns the global index of species *species* in phase *phase*. + """ + p = self.phaseIndex(phase) + + if isinstance(species, (str, unicode)): + k = self.phase(p).speciesIndex(species) + elif isinstance(species, (int, float)): + k = species + if not 0 <= k < self.nSpecies: + raise ValueError('Species index out of range') + else: + raise TypeError("'species' must be a string or number") + + return self.mix.speciesIndex(k, p) + + def nAtoms(self, k, m): + """ + Number of atoms of element *m* in the species with global index *k*. + The element may be referenced either by name or by index. + + >>> n = mix.nAtoms(3, 'H') + 4.0 + """ + if not 0 <= k < self.nSpecies: + raise IndexError('Species index ({}) out of range (0 < {})'.format(k, self.nSpecies)) + return self.mix.nAtoms(k, self.elementIndex(m)) + + property nPhases: + """Number of phases""" + def __get__(self): + return len(self._phases) + + def phase(self, n): + return self._phases[n] + + def phaseIndex(self, p): + """Index of the phase named *p*.""" + if isinstance(p, ThermoPhase): + p = p.name + + if isinstance(p, (int, float)): + if p == int(p) and 0 <= p < self.nPhases: + return int(p) + else: + raise IndexError("Phase index '{0}' out of range.".format(p)) + elif isinstance(p, (str, unicode)): + for i, phase in enumerate(self._phases): + if phase.name == p: + return i + raise KeyError("No such phase: '{0}'".format(p)) + + property phaseNames: + """Names of all phases in the order added.""" + def __get__(self): + return [phase.name for phase in self._phases] + + property T: + """ + The Temperature [K] of all phases in the mixture. When set, the + pressure of the mixture is held fixed. + """ def __get__(self): return self.mix.temperature() def __set__(self, T): self.mix.setTemperature(T) - property pressure: + property minTemp: + """ + The minimum temperature for which all species in multi-species + solutions have valid thermo data. Stoichiometric phases are not + considered in determining minTemp. + """ + def __get__(self): + return self.mix.minTemp() + + property maxTemp: + """ + The maximum temperature for which all species in multi-species + solutions have valid thermo data. Stoichiometric phases are not + considered in determining maxTemp. + """ + def __get__(self): + return self.mix.maxTemp() + + property P: + """The Pressure [Pa] of all phases in the mixture. When set, the + temperature of the mixture is held fixed.""" def __get__(self): return self.mix.pressure() def __set__(self, P): self.mix.setPressure(P) + + property charge: + """The total charge in Coulombs, summed over all phases.""" + def __get__(self): + return self.mix.charge() + + def phaseCharge(self, p): + """The charge of phase *p* in Coulumbs.""" + return self.mix.phaseCharge(self.phaseIndex(p)) + + def phaseMoles(self, p=None): + """ + Moles in phase *p*, if *p* is specified, otherwise the number of + moles in all phases. + """ + if p is None: + return [self.mix.phaseMoles(n) for n in range(self.nPhases)] + else: + return self.mix.phaseMoles(self.phaseIndex(p)) + + def setPhaseMoles(self, p, moles): + """ + Set the number of moles of phase *p* to *moles* + """ + self.mix.setPhaseMoles(self.phaseIndex(p), moles) + + def speciesMoles(self, species=None): + """ + Returns the number of moles of species *k* if *k* is specified, + or the number of of moles of each species otherwise. + """ + if species is not None: + return self.mix.speciesMoles(species) + + cdef np.ndarray[np.double_t, ndim=1] data = np.empty(self.nSpecies) + for k in range(self.nSpecies): + data[k] = self.mix.speciesMoles(k) + return data + + def setSpeciesMoles(self, moles): + """ + Set the moles of the species [kmol]. The moles may be specified either + as a string, or as an array. If an array is used, it must be + dimensioned at least as large as the total number of species in the + mixture. Note that the species may belong to any phase, and + unspecified species are set to zero. + + >>> mix.setSpeciesMoles('C(s):1.0, CH4:2.0, O2:0.2') + + """ + if isinstance(moles, (str, unicode)): + self.mix.setMolesByName(stringify(moles)) + return + + if len(moles) != self.nSpecies: + raise ValueError('mole array must be of length nSpecies') + + cdef np.ndarray[np.double_t, ndim=1] data = \ + np.ascontiguousarray(moles, dtype=np.double) + self.mix.setMoles(&data[0]) + + def elementMoles(self, e): + """ + Total number of moles of element *e*, summed over all species. + The element may be referenced either by index number or by name. + """ + return self.mix.elementMoles(self.elementIndex(e)) + + property chem_potentials: + """The chemical potentials of all species [J/kmol].""" + def __get__(self): + cdef np.ndarray[np.double_t, ndim=1] data = np.empty(self.nSpecies) + self.mix.getChemPotentials(&data[0]) + return data + + def equilibrate(self, XY, solver='vcs', rtol=1e-9, maxsteps=1000, + maxiter=100, estimateEquil=0, printlevel=0, loglevel=0): + """ + Set to a state of chemical equilibrium holding property pair *XY* + constant. This method uses a version of the VCS algorithm to find the + composition that minimizes the total Gibbs free energy of the mixture, + subject to element conservation constraints. For a description of the + theory, see Smith and Missen, "Chemical Reaction Equilibrium." + + :param XY: + A two-letter string, which must be one of the set:: + + ['TP', 'HP', 'SP'] + :param solver: + Set to either 'vcs' or 'gibbs' to choose implementation + of the solver to use. 'vcs' uses the solver implemented in the + C++ class 'VCSnonideal', and 'gibbs' uses the one implemented + in class 'MultiPhaseEquil'. + :param rtol: + Error tolerance. Iteration will continue until (Delta mu)/RT is + less than this value for each reaction. Note that this default is + very conservative, and good equilibrium solutions may be obtained + with larger error tolerances. + :param maxsteps: + Maximum number of steps to take while solving the equilibrium + problem for specified *T* and *P*. + :param maxiter: + Maximum number of temperature and/or pressure iterations. + This is only relevant if a property pair other than (T,P) is + specified. + :param estimateEquil: + Flag indicating whether the solver should estimate its own initial + condition. If 0, the initial mole fraction vector in the phase + objects are used as the initial condition. If 1, the initial mole + fraction vector is used if the element abundances are satisfied. + if -1, the initial mole fraction vector is thrown out, and an + estimate is formulated. + :param printlevel: + Determines the amount of output displayed during the solution + process. 0 indicates no output, while larger numbers produce + successively more verbose information. + :param loglevel: + Controls the amount of diagnostic output written to an HTML log + file. If loglevel = 0, no diagnostic output is written. For + values > 0, more detailed information is written to the log file as + loglevel increases. The default log file name is + "equilibrium_log.html", but if this file exists, the log + information will be written to "equilibrium_log{n}.html", + where {n} is an integer chosen to avoid overwriting existing + log files. + """ + if solver == 'vcs': + iSolver = 2 + elif solver == 'gibbs': + iSolver = 1 + else: + raise ValueError('Unrecognized equilibrium solver ' + 'specified: "{}"'.format(solver)) + + vcs_equilibrate(deref(self.mix), stringify(XY).c_str(), estimateEquil, + printlevel, iSolver, rtol, maxsteps, maxiter, loglevel) diff --git a/interfaces/cython/cantera/test/test_mixture.py b/interfaces/cython/cantera/test/test_mixture.py index 893801d9f..13e3266e6 100644 --- a/interfaces/cython/cantera/test/test_mixture.py +++ b/interfaces/cython/cantera/test/test_mixture.py @@ -1,20 +1,182 @@ import unittest - import cantera as ct from . import utilities + class TestMixture(utilities.CanteraTest): @classmethod def setUpClass(cls): cls.phase1 = ct.Solution('h2o2.xml') cls.phase2 = ct.Solution('air.xml') + def setUp(self): + self.mix = ct.Mixture([(self.phase1, 1.0), (self.phase2, 2.0)]) + + def test_sizes(self): + self.assertEqual(self.mix.nPhases, 2) + + self.assertEqual(self.mix.nSpecies, + self.phase1.nSpecies + self.phase2.nSpecies) + + E = set(self.phase1.elementNames) | set(self.phase2.elementNames) + self.assertEqual(len(E), self.mix.nElements) + + def test_elementIndex(self): + m_H = self.mix.elementIndex('H') + self.assertEqual(m_H, self.mix.elementIndex(m_H)) + + with self.assertRaises(ValueError): + self.mix.elementIndex('W') + + with self.assertRaises(ValueError): + self.mix.elementIndex(41) + + with self.assertRaises(TypeError): + self.mix.elementIndex(None) + + def test_speciesIndex(self): + names = self.mix.speciesNames + kOH = names.index('OH') + kN2 = names.index('N2') + self.assertEqual(self.mix.speciesName(kOH), 'OH') + self.assertEqual(self.mix.speciesName(kN2), 'N2') + + self.assertEqual(self.mix.speciesIndex(0, 'OH'), kOH) + self.assertEqual(self.mix.speciesIndex(self.phase1, 'OH'), kOH) + self.assertEqual(self.mix.speciesIndex(self.phase1.name, 'OH'), kOH) + self.assertEqual(self.mix.speciesIndex(0, self.phase1.speciesIndex('OH')), kOH) + self.assertEqual(self.mix.speciesIndex(1, self.phase2.speciesIndex('N2')), kN2) + self.assertEqual(self.mix.speciesIndex(1, 'N2'), kN2) + + with self.assertRaises(IndexError): + self.mix.speciesIndex(3, 'OH') + + with self.assertRaises(ValueError): + self.mix.speciesIndex(1, 'OH') + + with self.assertRaises(ValueError): + self.mix.speciesIndex(0, -2) + + with self.assertRaises(ValueError): + self.mix.speciesIndex(1, 'CO2') + + def test_nAtoms(self): + names = self.mix.speciesNames + kOH = names.index('OH') + kN2 = names.index('N2') + mH = self.mix.elementIndex('H') + mN = self.mix.elementIndex('N') + + self.assertEqual(self.mix.nAtoms(kOH, 'H'), 1) + self.assertEqual(self.mix.nAtoms(kOH, 'O'), 1) + self.assertEqual(self.mix.nAtoms(kOH, mH), 1) + self.assertEqual(self.mix.nAtoms(kOH, mN), 0) + + self.assertEqual(self.mix.nAtoms(kN2, mN), 2) + self.assertEqual(self.mix.nAtoms(kN2, mH), 0) + + def test_phase(self): + self.assertEqual(self.phase1, self.mix.phase(0)) + self.assertEqual(self.phase2, self.mix.phase(1)) + + phaseNames = self.mix.phaseNames + self.assertEqual(len(phaseNames), self.mix.nPhases) + self.assertEqual(phaseNames[0], self.phase1.name) + self.assertEqual(phaseNames[1], self.phase2.name) + + def test_phaseIndex(self): + self.assertEqual(self.mix.phaseIndex(self.phase1), 0) + self.assertEqual(self.mix.phaseIndex(self.phase2), 1) + self.assertEqual(self.mix.phaseIndex(self.phase2.name), 1) + self.assertEqual(self.mix.phaseIndex(1), 1) + + with self.assertRaises(KeyError): + self.mix.phaseIndex('foobar') + + with self.assertRaises(IndexError): + self.mix.phaseIndex(2) + def test_properties(self): - mix = ct.Mixture([(self.phase1, 1.0), (self.phase2, 2.0)]) - self.assertEqual(mix.nSpecies, self.phase1.nSpecies + self.phase2.nSpecies) + self.mix.T = 350 + self.assertEqual(self.mix.T, 350) - mix.temperature = 350 - self.assertEqual(mix.temperature, 350) + self.mix.P = 2e5 + self.assertEqual(self.mix.P, 2e5) + self.assertEqual(self.mix.T, 350) - mix.pressure = 2e5 - self.assertEqual(mix.pressure, 2e5) + self.assertGreater(self.mix.maxTemp, self.mix.minTemp) + + def test_charge(self): + C = sum(self.mix.phaseCharge(i) for i in range(self.mix.nPhases)) + self.assertEqual(self.mix.charge, C) + + def test_phaseMoles(self): + M = self.mix.phaseMoles() + self.assertEqual(M[0], self.mix.phaseMoles(0)) + self.assertEqual(M[1], self.mix.phaseMoles('air')) + + self.mix.setPhaseMoles('air', 4) + self.assertEqual(self.mix.phaseMoles(1), 4) + + def test_speciesMoles(self): + self.mix.setSpeciesMoles('H2:1.0, N2:4.0') + P = self.mix.phaseMoles() + S = self.mix.speciesMoles() + + self.assertEqual(P[0], 1) + self.assertEqual(P[1], 4) + + self.assertEqual(S[self.mix.speciesIndex(0, 'H2')], 1) + self.assertEqual(S[self.mix.speciesIndex(1, 'N2')], 4) + + S[2] = 7 + self.mix.setSpeciesMoles(S) + self.assertNear(self.mix.speciesMoles(2), S[2]) + self.assertNear(self.mix.phaseMoles(0), sum(S[:self.phase1.nSpecies])) + + with self.assertRaises(ValueError): + self.mix.setSpeciesMoles((1,2,3)) + + with self.assertRaises(TypeError): + self.mix.setSpeciesMoles(9) + + def test_elementMoles(self): + self.mix.setSpeciesMoles('H2:1.0, OH:4.0') + + self.assertNear(self.mix.elementMoles('H'), 6) + self.assertNear(self.mix.elementMoles('O'), 4) + self.assertNear(self.mix.elementMoles('N'), 0) + + def test_chem_potentials(self): + C = self.mix.chem_potentials + C1 = self.phase1.chem_potentials + C2 = self.phase2.chem_potentials + + self.assertArrayNear(C[:self.phase1.nSpecies], C1) + self.assertArrayNear(C[self.phase1.nSpecies:], C2) + + def test_equilibrate1(self): + self.mix.setSpeciesMoles('H2:1.0, O2:0.5, N2:1.0') + self.mix.T = 400 + self.mix.P = 2 * ct.OneAtm + + E1 = [self.mix.elementMoles(m) for m in range(self.mix.nElements)] + self.mix.equilibrate('TP') + + E2 = [self.mix.elementMoles(m) for m in range(self.mix.nElements)] + self.assertArrayNear(E1, E2) + self.assertNear(self.mix.T, 400) + self.assertNear(self.mix.P, 2 * ct.OneAtm) + + def test_equilibrate2(self): + self.mix.setSpeciesMoles('H2:1.0, O2:0.5, N2:1.0') + self.mix.T = 400 + self.mix.P = 2 * ct.OneAtm + + E1 = [self.mix.elementMoles(m) for m in range(self.mix.nElements)] + self.mix.equilibrate('TP', solver='gibbs') + + E2 = [self.mix.elementMoles(m) for m in range(self.mix.nElements)] + self.assertArrayNear(E1, E2) + self.assertNear(self.mix.T, 400) + self.assertNear(self.mix.P, 2 * ct.OneAtm)