From 987e1ddbb0563f3511793b5ee6c87d007f37f358 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Wed, 1 Dec 2004 22:57:23 +0000 Subject: [PATCH] support for multiphase equilibrium --- Cantera/clib/src/Makefile.in | 2 +- Cantera/clib/src/ctmultiphase.cpp | 176 +++++ Cantera/clib/src/ctmultiphase.h | 28 + Cantera/cxx/include/equilibrium.h | 4 +- Cantera/python/Cantera/Mixture.py | 83 +++ Cantera/python/src/ctmultiphase_methods.cpp | 234 +++++++ Cantera/python/src/methods.h | 17 + Cantera/python/src/pycantera.cpp | 2 + Cantera/src/DenseMatrix.cpp | 2 +- Cantera/src/Makefile.in | 2 +- Cantera/src/MultiPhaseEquil.cpp | 669 ++++++++++++++++++++ Cantera/src/MultiPhaseEquil.h | 122 ++++ Cantera/src/phasereport.cpp | 12 +- config.h.in | 14 +- 14 files changed, 1355 insertions(+), 12 deletions(-) create mode 100644 Cantera/clib/src/ctmultiphase.cpp create mode 100644 Cantera/clib/src/ctmultiphase.h create mode 100644 Cantera/python/Cantera/Mixture.py create mode 100644 Cantera/python/src/ctmultiphase_methods.cpp create mode 100644 Cantera/src/MultiPhaseEquil.cpp create mode 100644 Cantera/src/MultiPhaseEquil.h diff --git a/Cantera/clib/src/Makefile.in b/Cantera/clib/src/Makefile.in index ecf7bd14e..d606fd2d3 100755 --- a/Cantera/clib/src/Makefile.in +++ b/Cantera/clib/src/Makefile.in @@ -15,7 +15,7 @@ SUFFIXES= .cpp .d .o CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) OBJS = ct.o Storage.o ctsurf.o ctrpath.o \ - ctreactor.o ctfunc.o ctxml.o ctonedim.o + ctreactor.o ctfunc.o ctxml.o ctonedim.o ctmultiphase.o DEPENDS = $(OBJS:.o=.d) diff --git a/Cantera/clib/src/ctmultiphase.cpp b/Cantera/clib/src/ctmultiphase.cpp new file mode 100644 index 000000000..327d8a847 --- /dev/null +++ b/Cantera/clib/src/ctmultiphase.cpp @@ -0,0 +1,176 @@ + +// Cantera includes +#include "MultiPhase.h" +#include "MultiPhaseEquil.h" + +#include "Cabinet.h" +#include "Storage.h" + +// Build as a DLL under Windows +#ifdef WIN32 +#define DLL_EXPORT __declspec(dllexport) +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#else +#define DLL_EXPORT +#endif + +// Values returned for error conditions +#define ERR -999 +#define DERR -999.999 + +typedef MultiPhase mix_t; + +Cabinet* Cabinet::__storage = 0; + +inline mix_t* _mix(int i) { + return Cabinet::cabinet()->item(i); +} + +inline ThermoPhase* _th(int n) { + return Storage::__storage->__thtable[n]; +} + +static bool checkSpecies(int i, int k) { + try { + if (k < 0 || k >= _mix(i)->nSpecies()) + throw CanteraError("checkSpecies", + "illegal species index ("+int2str(k)+") "); + return true; + } + catch (CanteraError) { + return false; + } +} + +static bool checkElement(int i, int m) { + try { + if (m < 0 || m >= _mix(i)->nElements()) + throw CanteraError("checkElement", + "illegal element index ("+int2str(m)+") "); + return true; + } + catch (CanteraError) { + return false; + } +} + +static bool checkPhase(int i, int n) { + try { + if (n < 0 || n >= _mix(i)->nPhases()) + throw CanteraError("checkPhase", + "illegal phase index ("+int2str(n)+") "); + return true; + } + catch (CanteraError) { + return false; + } +} + +extern "C" { + + int DLL_EXPORT mix_new() { + mix_t* m = new MultiPhase(); + return Cabinet::cabinet()->add(m); + } + + int DLL_EXPORT mix_del(int i) { + Cabinet::cabinet()->del(i); + return 0; + } + + int DLL_EXPORT mix_copy(int i) { + return Cabinet::cabinet()->newCopy(i); + } + + int DLL_EXPORT mix_assign(int i, int j) { + return Cabinet::cabinet()->assign(i,j); + } + + int DLL_EXPORT mix_addPhase(int i, int j, double moles) { + _mix(i)->addPhase(_th(j), moles); + return 0; + } + + int DLL_EXPORT mix_nElements(int i) { + return _mix(i)->nElements(); + } + + int DLL_EXPORT mix_nSpecies(int i) { + return _mix(i)->nSpecies(); + } + + doublereal DLL_EXPORT mix_nAtoms(int i, int k, int m) { + bool ok = (checkSpecies(i,k) && checkElement(i,m)); + if (ok) + return _mix(i)->nAtoms(k,m); + else + return DERR; + } + + doublereal DLL_EXPORT mix_phaseMoles(int i, int n) { + if (!checkPhase(i, n)) return DERR; + return _mix(i)->phaseMoles(n); + } + + int DLL_EXPORT mix_setPhaseMoles(int i, int n, double v) { + if (!checkPhase(i, n)) return ERR; + if (v < 0.0) return -1; + _mix(i)->setPhaseMoles(n, v); + return 0; + } + + int DLL_EXPORT mix_setTemperature(int i, double t) { + if (t < 0.0) return -1; + _mix(i)->setTemperature(t); + return 0; + } + + doublereal DLL_EXPORT mix_temperature(int i) { + return _mix(i)->temperature(); + } + + int DLL_EXPORT mix_setPressure(int i, double p) { + if (p < 0.0) return -1; + _mix(i)->setPressure(p); + return 0; + } + + doublereal DLL_EXPORT mix_pressure(int i) { + return _mix(i)->pressure(); + } + + doublereal DLL_EXPORT mix_speciesMoles(int i, int k) { + if (!checkSpecies(i,k)) return DERR; + return _mix(i)->speciesMoles(k); + } + + doublereal DLL_EXPORT mix_elementMoles(int i, int m) { + if (!checkElement(i,m)) return DERR; + return _mix(i)->elementMoles(m); + } + + + doublereal DLL_EXPORT mix_equilibrate(int i, char* XY, + doublereal err, int maxiter) { + try { + return equilibrate(*_mix(i), XY, err, maxiter); + } + catch (CanteraError) { + return DERR; + } + } + + int DLL_EXPORT mix_getChemPotentials(int i, int lenmu, double* mu) { + try { + if (lenmu < _mix(i)->nSpecies()) + throw CanteraError("getChemPotentials","array too small"); + _mix(i)->getChemPotentials(mu); + return 0; + } + catch (CanteraError) { + return -1; + } + } + +} diff --git a/Cantera/clib/src/ctmultiphase.h b/Cantera/clib/src/ctmultiphase.h new file mode 100644 index 000000000..dce11e320 --- /dev/null +++ b/Cantera/clib/src/ctmultiphase.h @@ -0,0 +1,28 @@ +#ifndef CTC_MULTIPHASE_H +#define CTC_MULTIPHASE_H + +#include "clib_defs.h" + +extern "C" { + + int DLL_IMPORT mix_new(); + int DLL_IMPORT mix_del(int i); + int DLL_IMPORT mix_copy(int i); + int DLL_IMPORT mix_assign(int i, int j); + int DLL_IMPORT mix_addPhase(int i, int j, double moles); + int DLL_IMPORT mix_nElements(int i); + int DLL_IMPORT mix_nSpecies(int i); + int DLL_IMPORT mix_setTemperature(int i, double t); + double DLL_IMPORT mix_temperature(int i); + int DLL_IMPORT mix_setPressure(int i, double p); + double DLL_IMPORT mix_pressure(int i); + double DLL_IMPORT mix_nAtoms(int i, int k, int m); + double DLL_IMPORT mix_phaseMoles(int i, int n); + int DLL_IMPORT mix_setPhaseMoles(int i, int n, double v); + double DLL_IMPORT mix_speciesMoles(int i, int k); + double DLL_IMPORT mix_elementMoles(int i, int m); + double DLL_IMPORT mix_equilibrate(int i, char* XY, + double err, int maxiter); + int DLL_IMPORT mix_getChemPotentials(int i, int lenmu, double* mu); +} +#endif diff --git a/Cantera/cxx/include/equilibrium.h b/Cantera/cxx/include/equilibrium.h index edf6a1d4f..cb77bda99 100755 --- a/Cantera/cxx/include/equilibrium.h +++ b/Cantera/cxx/include/equilibrium.h @@ -4,7 +4,9 @@ #ifndef CT_EQUIL_INCL #define CT_EQUIL_INCL #include "kernel/ChemEquil.h" -#include "kernel/MultiPhaseEquil.h" +//#ifdef DEV_EQUIL +//#include "kernel/MultiPhaseEquil.h" +//#endif #endif diff --git a/Cantera/python/Cantera/Mixture.py b/Cantera/python/Cantera/Mixture.py new file mode 100644 index 000000000..b441c85bd --- /dev/null +++ b/Cantera/python/Cantera/Mixture.py @@ -0,0 +1,83 @@ +import _cantera +import types +from Numeric import zeros + +class Mixture: + """Class Mixture represents mixtures of one or more phases of matter.""" + + def __init__(self, phases=[]): + self.__mixid = _cantera.mix_new() + self._spnames = [] + self._phases = [] + if phases: + for p in phases: + try: + ph = p[0] + moles = p[1] + except: + ph = p + moles = 0 + self.addPhase(ph, moles) + self._phases.append(ph) + + + def __del__(self): + _cantera.mix_del(self.__mixid) + + def __repr__(self): + s = '' + for p in range(len(self._phases)): + s += '\n******************* Phase '+`p`+' ******************************\n' + s += '\n Moles: '+`self.phaseMoles(p)`+'\n' + s += self._phases[p].__repr__()+'\n\n' + return s + + def addPhase(self, phase = None, moles = 0.0): + for k in range(phase.nSpecies()): + self._spnames.append(phase.speciesName(k)) + _cantera.mix_addPhase(self.__mixid, phase.thermo_hndl(), moles) + + def nElements(self): + """Total number of elements present in the mixture.""" + return _cantera.mix_nElements(self.__mixid) + def nSpecies(self): + """Total number of species present in the mixture. This is the + sum of the numbers of species in each phase.""" + return _cantera.mix_nSpecies(self.__mixid) + def speciesName(self, k): + return self._spnames[k] + def speciesIndex(self, species): + if type(species) == types.StringType: + return self._spnames.index(species) + else: + return species + def nAtoms(self, k, m): + """Number of atoms of element m in species k.""" + return _cantera.mix_nAtoms(self.__mixid, k, m) + def setTemperature(self, t): + return _cantera.mix_setTemperature(self.__mixid, t) + def temperature(self): + return _cantera.mix_temperature(self.__mixid) + def setPressure(self, p): + return _cantera.mix_setPressure(self.__mixid, p) + def pressure(self): + return _cantera.mix_pressure(self.__mixid) + def phaseMoles(self, n): + """Moles of phase n.""" + return _cantera.mix_phaseMoles(self.__mixid, n) + def setPhaseMoles(self, n, moles): + """Set the moles of phase n.""" + return _cantera.mix_setPhaseMoles(self.__mixid, n, moles) + def speciesMoles(self, species): + """Moles of species k.""" + k = self.speciesIndex(species) + return _cantera.mix_speciesMoles(self.__mixid, k) + def elementMoles(self, m): + return _cantera.mix_elementMoles(self.__mixid, m) + def chemPotentials(self): + mu = zeros(self.nSpecies(),'d') + _cantera.mix_getChemPotentials(self.__mixid, mu) + return mu + def equilibrate(self, XY = "TP", err = 1.0e-9, maxiter = 1000): + return _cantera.mix_equilibrate(self.__mixid, XY, err, maxiter) + diff --git a/Cantera/python/src/ctmultiphase_methods.cpp b/Cantera/python/src/ctmultiphase_methods.cpp new file mode 100644 index 000000000..9a42d7d3e --- /dev/null +++ b/Cantera/python/src/ctmultiphase_methods.cpp @@ -0,0 +1,234 @@ + +static PyObject * +py_mix_new(PyObject *self, PyObject *args) +{ + int _val; + _val = mix_new(); + return Py_BuildValue("i",_val); +} + + +static PyObject * +py_mix_del(PyObject *self, PyObject *args) +{ + int _val; + int i; + if (!PyArg_ParseTuple(args, "i:mix_del", &i)) + return NULL; + + _val = mix_del(i); + if (int(_val) < 0) return reportCanteraError(); + return Py_BuildValue("i",_val); +} + + +static PyObject * +py_mix_addPhase(PyObject *self, PyObject *args) +{ + int _val; + int i; + int j; + double moles; + if (!PyArg_ParseTuple(args, "iid:mix_addPhase", &i, &j, &moles)) + return NULL; + + _val = mix_addPhase(i,j,moles); + if (int(_val) < 0) return reportCanteraError(); + return Py_BuildValue("i",_val); +} + + +static PyObject * +py_mix_nElements(PyObject *self, PyObject *args) +{ + int _val; + int i; + if (!PyArg_ParseTuple(args, "i:mix_nElements", &i)) + return NULL; + + _val = mix_nElements(i); + if (int(_val) < -900) return reportCanteraError(); + return Py_BuildValue("i",_val); +} + +static PyObject * +py_mix_nSpecies(PyObject *self, PyObject *args) +{ + int _val; + int i; + if (!PyArg_ParseTuple(args, "i:mix_nSpecies", &i)) + return NULL; + + _val = mix_nSpecies(i); + if (int(_val) < -900) return reportCanteraError(); + return Py_BuildValue("i",_val); +} + + +static PyObject * +py_mix_nAtoms(PyObject *self, PyObject *args) +{ + double _val; + int i; + int k; + int m; + if (!PyArg_ParseTuple(args, "iii:mix_nAtoms", &i, &k, &m)) + return NULL; + + _val = mix_nAtoms(i,k,m); + if (int(_val) < -900) return reportCanteraError(); + return Py_BuildValue("d",_val); +} + +static PyObject * +py_mix_setTemperature(PyObject *self, PyObject *args) +{ + int _val; + int i; + double t; + if (!PyArg_ParseTuple(args, "id:mix_setTemperature", &i, &t)) + return NULL; + + _val = mix_setTemperature(i,t); + if (int(_val) == -1) return reportCanteraError(); + return Py_BuildValue("i",_val); +} + + +static PyObject * +py_mix_temperature(PyObject *self, PyObject *args) +{ + double _val; + int i; + if (!PyArg_ParseTuple(args, "i:mix_temperature", &i)) + return NULL; + + _val = mix_temperature(i); + if (int(_val) == -1) return reportCanteraError(); + return Py_BuildValue("d",_val); +} + + +static PyObject * +py_mix_setPressure(PyObject *self, PyObject *args) +{ + int _val; + int i; + double p; + if (!PyArg_ParseTuple(args, "id:mix_setPressure", &i, &p)) + return NULL; + + _val = mix_setPressure(i,p); + if (int(_val) == -1) return reportCanteraError(); + return Py_BuildValue("i",_val); +} + + +static PyObject * +py_mix_pressure(PyObject *self, PyObject *args) +{ + double _val; + int i; + if (!PyArg_ParseTuple(args, "i:mix_pressure", &i)) + return NULL; + + _val = mix_pressure(i); + if (int(_val) == -1) return reportCanteraError(); + return Py_BuildValue("d",_val); +} + +static PyObject * +py_mix_phaseMoles(PyObject *self, PyObject *args) +{ + double _val; + int i; + int n; + if (!PyArg_ParseTuple(args, "ii:mix_phaseMoles", &i, &n)) + return NULL; + + _val = mix_phaseMoles(i,n); + if (int(_val) < -900) return reportCanteraError(); + return Py_BuildValue("d",_val); +} + + +static PyObject * +py_mix_setPhaseMoles(PyObject *self, PyObject *args) +{ + int _val; + int i; + int n; + double v; + if (!PyArg_ParseTuple(args, "iid:mix_setPhaseMoles", &i, &n, &v)) + return NULL; + + _val = mix_setPhaseMoles(i,n,v); + if (int(_val) < 0) return reportCanteraError(); + return Py_BuildValue("i",_val); +} + + +static PyObject * +py_mix_speciesMoles(PyObject *self, PyObject *args) +{ + double _val; + int i; + int k; + if (!PyArg_ParseTuple(args, "ii:mix_speciesMoles", &i, &k)) + return NULL; + + _val = mix_speciesMoles(i,k); + if (int(_val) < -900) return reportCanteraError(); + return Py_BuildValue("d",_val); +} + + +static PyObject * +py_mix_elementMoles(PyObject *self, PyObject *args) +{ + double _val; + int i; + int m; + if (!PyArg_ParseTuple(args, "ii:mix_elementMoles", &i, &m)) + return NULL; + + _val = mix_elementMoles(i,m); + //if (int(_val) < -900) return reportCanteraError(); + return Py_BuildValue("d",_val); +} + +static PyObject * +py_mix_equilibrate(PyObject *self, PyObject *args) +{ + double _val; + int i; + char* XY; + double err; + int maxiter; + if (!PyArg_ParseTuple(args, "isdi:mix_equilibrate", &i, &XY, &err, &maxiter)) + return NULL; + + _val = mix_equilibrate(i,XY,err,maxiter); + if (int(_val) < -900) return reportCanteraError(); + return Py_BuildValue("d",_val); +} + + +static PyObject * +py_mix_getChemPotentials(PyObject *self, PyObject *args) +{ + int i; + int _val; + PyObject* mu; + if (!PyArg_ParseTuple(args, "iO:mix_getChemPotentials", &i, &mu)) + return NULL; + + PyArrayObject* mu_array = (PyArrayObject*)mu; + double* mu_data = (double*)mu_array->data; + int mu_len = mu_array->dimensions[0]; + + _val = mix_getChemPotentials(i, mu_len, mu_data); + if (int(_val) < 0) return reportCanteraError(); + return Py_BuildValue("i",_val); +} + diff --git a/Cantera/python/src/methods.h b/Cantera/python/src/methods.h index 8ce3b131e..b0fc10da5 100644 --- a/Cantera/python/src/methods.h +++ b/Cantera/python/src/methods.h @@ -251,6 +251,23 @@ static PyMethodDef ct_methods[] = { {"func_del", py_func_del, METH_VARARGS}, {"func_value", py_func_value, METH_VARARGS}, + {"mix_new", py_mix_new, METH_VARARGS}, + {"mix_del", py_mix_del, METH_VARARGS}, + {"mix_addPhase", py_mix_addPhase, METH_VARARGS}, + {"mix_nElements", py_mix_nElements, METH_VARARGS}, + {"mix_nSpecies", py_mix_nSpecies, METH_VARARGS}, + {"mix_nAtoms", py_mix_nAtoms, METH_VARARGS}, + {"mix_setTemperature", py_mix_setTemperature, METH_VARARGS}, + {"mix_temperature", py_mix_temperature, METH_VARARGS}, + {"mix_setPressure", py_mix_setPressure, METH_VARARGS}, + {"mix_pressure", py_mix_pressure, METH_VARARGS}, + {"mix_phaseMoles", py_mix_phaseMoles, METH_VARARGS}, + {"mix_setPhaseMoles", py_mix_setPhaseMoles, METH_VARARGS}, + {"mix_speciesMoles", py_mix_speciesMoles, METH_VARARGS}, + {"mix_elementMoles", py_mix_elementMoles, METH_VARARGS}, + {"mix_equilibrate", py_mix_equilibrate, METH_VARARGS}, + {"mix_getChemPotentials", py_mix_getChemPotentials, METH_VARARGS}, + #ifdef INCL_USER_PYTHON #include "usermethods.h" #endif diff --git a/Cantera/python/src/pycantera.cpp b/Cantera/python/src/pycantera.cpp index 5d10e6bd1..8819bdeef 100644 --- a/Cantera/python/src/pycantera.cpp +++ b/Cantera/python/src/pycantera.cpp @@ -24,6 +24,7 @@ #include "ctreactor.h" #include "ctfunc.h" #include "ctonedim.h" +#include "ctmultiphase.h" #include using namespace std; @@ -46,6 +47,7 @@ static PyObject *ErrorObject; #include "ctreactor_methods.cpp" #include "ctfunc_methods.cpp" #include "ctonedim_methods.cpp" +#include "ctmultiphase_methods.cpp" #ifdef INCL_USER_PYTHON #include "ctuser.h" diff --git a/Cantera/src/DenseMatrix.cpp b/Cantera/src/DenseMatrix.cpp index 681fb69f6..04e6d9795 100755 --- a/Cantera/src/DenseMatrix.cpp +++ b/Cantera/src/DenseMatrix.cpp @@ -118,7 +118,7 @@ namespace Cantera { */ void multiply(const DenseMatrix& A, const double* b, double* prod) { ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, - static_cast(A.nRows()), static_cast(A.nRows()), 1.0, + static_cast(A.nRows()), static_cast(A.nColumns()), 1.0, A.begin(), static_cast(A.nRows()), b, 1, 0.0, prod, 1); } diff --git a/Cantera/src/Makefile.in b/Cantera/src/Makefile.in index 6e0dd6af9..bf83a8994 100755 --- a/Cantera/src/Makefile.in +++ b/Cantera/src/Makefile.in @@ -40,7 +40,7 @@ HETEROKIN = InterfaceKinetics.o ImplicitSurfChem.o SurfPhase.o EdgeKinetics.o $( CK = $(KINETICS) # chemical equilibrium -EQUIL = ChemEquil.o sort.o $(THERMO) +EQUIL = ChemEquil.o MultiPhaseEquil.o sort.o $(THERMO) # reaction path analysis RPATH = Group.o ReactionPath.o diff --git a/Cantera/src/MultiPhaseEquil.cpp b/Cantera/src/MultiPhaseEquil.cpp new file mode 100644 index 000000000..962006130 --- /dev/null +++ b/Cantera/src/MultiPhaseEquil.cpp @@ -0,0 +1,669 @@ +#include "MultiPhaseEquil.h" +#include "MultiPhase.h" +#include "sort.h" +#include "recipes.h" + +#include +#include +using namespace std; + + +namespace Cantera { + + const doublereal TINY = 1.0e-20; + + /// Used to print reaction equations. Given a stoichiometric + /// coefficient 'nu' and a chemical symbol 'sym', return a string + /// for this species in the reaction. + /// @param first if this is false, then a " + " string will be + /// added to the beginning of the string. + /// @param nu Stoichiometric coefficient. May be positive or negative. + /// @param sym Species chemical symbol. + /// + static string coeffString(bool first, doublereal nu, string sym) { + if (nu == 0.0) return ""; + string strt = " + "; + if (first) strt = ""; + if (nu == 1.0 || nu == -1.0) + return strt + sym; + string s = fp2str(fabs(nu)); + return strt + s + " " + sym; + } + + + /// Constructor. Construct a multiphase equilibrium manager for + /// a multiphase mixture. + /// @param mix Pointer to a multiphase mixture object. + MultiPhaseEquil::MultiPhaseEquil(mix_t* mix) : m_mix(mix) + { + // the multi-phase mixture + m_mix = mix; + + // store some mixture parameters locally + m_nel_mix = mix->nElements(); + m_nsp_mix = mix->nSpecies(); + m_np = mix->nPhases(); + m_press = mix->pressure(); + m_temp = mix->temperature(); + + index_t m, k; + m_nel = 0; + m_nsp = 0; + m_incl_species.resize(m_nsp_mix,1); + m_incl_element.resize(m_nel_mix,1); + for (m = 0; m < m_nel_mix; m++) { + if (m_mix->elementMoles(m) <= 0.0) { + m_incl_element[m] = 0; + for (k = 0; k < m_nsp_mix; k++) { + if (m_mix->nAtoms(k,m) != 0.0) { + m_incl_species[k] = 0; + } + } + } + } + for (m = 0; m < m_nel_mix; m++) { + if (m_incl_element[m] == 1) { + m_nel++; + m_element.push_back(m); + } + } + for (k = 0; k < m_nsp_mix; k++) { + if (m_incl_species[k] ==1) { + m_nsp++; + m_species.push_back(k); + } + } + //cout << "nsp = " << m_nsp << endl; + //cout << m_element << endl << m_species << endl; + + // some work arrays for internal use + m_work.resize(m_nsp); + m_work2.resize(m_nsp); + m_mu.resize(m_nsp_mix); + + // number of moles of each species + m_moles.resize(m_nsp); + m_lastmoles.resize(m_nsp); + m_dxi.resize(m_nsp - m_nel); + + index_t ik; + for (ik = 0; ik < m_nsp; ik++) { + m_moles[ik] = m_mix->speciesMoles(m_species[ik]); + } + + // Delta G / RT for each reaction + m_deltaG_RT.resize(m_nsp - m_nel, 0.0); + m_majorsp.resize(m_nsp); + m_sortindex.resize(m_nsp,0); + m_lastsort.resize(m_nel); + m_solnrxn.resize(m_nsp - m_nel); + m_A.resize(m_nel, m_nsp, 0.0); + m_N.resize(m_nsp, m_nsp - m_nel); + m_order.resize(m_nsp, 0); + + setInitialMoles(); + computeN(); + + // make sure the components are non-zero + for (k = 0; k < m_nel; k++) { + if (m_moles[m_order[k]] <= 0.0) { + m_moles[m_order[k]] = 1.0e-17; + } + } + vector_fp dxi(m_nsp - m_nel, 1.0e-20); + multiply(m_N, dxi.begin(), m_work.begin()); + unsort(m_work); + + for (k = 0; k < m_nsp; k++) { + m_moles[k] += m_work[k]; + m_lastmoles[k] = m_moles[k]; + if (m_mix->solutionSpecies(m_species[k])) + m_dsoln.push_back(1); + else + m_dsoln.push_back(0); + } + setMoles(); + } + + void MultiPhaseEquil::setMoles() { + vector_fp n(m_nsp_mix, 0.0); + index_t k; + for (k = 0; k < m_nsp; k++) { + n[m_species[k]] = m_moles[k]; + } + m_mix->setMoles(n.begin()); + } + + /** + * Estimate the initial mole fractions. Uses the Simplex method + * to estimate the initial number of moles of each species. The + * linear Gibbs minimization problem is solved, neglecting the + * free energy of mixing terms. This procedure produces a good + * estimate of the low-temperature equilibrium composition. + * + * @param s phase object + * @param elementMoles vector of elemental moles + */ + int MultiPhaseEquil::setInitialMoles() { + int m, n; + double lp = log(m_press/OneAtm); + + DenseMatrix aa(m_nel+2, m_nsp+1, 0.0); + + // first column contains fixed element moles + for (m = 0; m < m_nel; m++) { + aa(m+1,0) = m_mix->elementMoles(m_element[m]); + } + + // get the array of non-dimensional Gibbs functions for the pure + // species + m_mix->getStandardChemPotentials(m_mu.begin()); + + int kpp = 0; + doublereal rt = GasConstant * m_temp; + for (int k = 0; k < m_nsp; k++) { + kpp++; + aa(0, kpp) = -m_mu[m_species[k]]/rt; + aa(0, kpp) -= m_dsoln[k]*lp; // ideal gas + for (int q = 0; q < m_nel; q++) + aa(q+1, kpp) = -m_mix->nAtoms(m_species[k], m_element[q]); + } + + integer mp = m_nel+2; // parameters for SIMPLX + integer np = m_nsp+1; + integer m1 = 0; + integer m2 = 0; + integer m3 = m_nel; + integer icase=0; + + vector_int iposv(m_nel); + vector_int izrov(m_nsp); + + // solve the linear programming problem + + simplx_(&aa(0,0), &m_nel, &m_nsp, &mp, &np, &m1, &m2, &m3, + &icase, izrov.begin(), iposv.begin()); + + fill(m_moles.begin(), m_moles.end(), 0.0); + for (n = 0; n < m_nel; n++) { + int ksp = 0; + int ip = iposv[n] - 1; + for (int k = 0; k < m_nsp; k++) { + if (ip == ksp) { + m_moles[k] = aa(n+1, 0); + } + ksp++; + } + } + setMoles(); + return icase; + } + + + /// This method finds a set of constituent species and a complete + /// set of formation reactions for the non-constituents in terms + /// of the constituents. Note that in most cases, many different + /// constituent sets are possible, and therefore neither the + /// constituents returned by this method nor the formation + /// reactions are unique. The algorithm used here is described in + /// Smith and Missen, Chemical Reaction Equilibrium Analysis. + /// + /// The constituent species are taken to be the first M species + /// in array 'species' that have linearly-independent compositions. + /// + /// Arguments: + /// + /// On entry, vector species shold contain species index numbers + /// in the order of decreasing desirability as a constituent. For + /// example, if it is desired to choose the constituents from + /// among the major species, this array might list species index + /// numbers in decreasing order of mole fraction. If array + /// 'species' does not have length = nSpecies(), then the species + /// will be considered as candidates to be constituents in + /// declaration order, beginning with the first phase added. + /// + /// On return, the first M entries of array 'species' contain the index + /// numbers of the constituent species. + /// + /// Matrix nu is an output array that contains the stoichiometric + /// coefficents for a set of K - M formation reactions for the + /// non-constituent species, such that nu(k,i) is the net + /// stoichiometric coefficent of species k in reaction i. Matrix + /// nu will be resized to (K, K-M) and its initial values, if + /// any, will be erased. + + void MultiPhaseEquil::getComponents(const vector_int& order) { + int m, n, k, j; + + // if the input species array has the wrong size, ignore it + // and consider the species for constituents in declarationi order. + if (order.size() != m_nsp) { + for (k = 0; k < m_nsp; k++) m_order[k] = k; + } + else { + for (k = 0; k < m_nsp; k++) m_order[k] = order[k]; + } + // cout << m_order << endl; + doublereal tmp; + index_t itmp; + + index_t nRows = m_nel; + index_t nColumns = m_nsp; + doublereal fctr; + + // set up the atomic composition matrix + for (m = 0; m < nRows; m++) { + for (k = 0; k < nColumns; k++) { + m_A(m, k) = m_mix->nAtoms(m_species[m_order[k]], m_element[m]); + } + } + + // Do Gauss elimination + for (m = 0; m < nRows; m++) { + // if a pivot is zero, exchange columns + if (m_A(m,m) == 0.0) { + for (k = m+1; k < nColumns; k++) { + if (m_A(m,k) != 0.0) { + for (n = 0; n < nRows; n++) { + tmp = m_A(n,m); + m_A(n, m) = m_A(n, k); + m_A(n, k) = tmp; + } + // exchange the species labels on the columns + itmp = m_order[m]; + m_order[m] = m_order[k]; + m_order[k] = itmp; + break; + } + } + // throw an exception if the entire row is zero + if (k >= m_nsp) + throw CanteraError("getComponents","all zeros!"); + } + + // scale row m so that the diagonal element is unity + fctr = 1.0/m_A(m,m); + for (k = 0; k < nColumns; k++) { + m_A(m,k) *= fctr; + } + + // subtract A(n,m)/A(m,m) * (row m) from row n, so that + // A(n,m) = 0. + for (n = m+1; n < m_nel; n++) { + fctr = m_A(n,m)/m_A(m,m); + for (k = 0; k < m_nsp; k++) { + m_A(n,k) -= m_A(m,k)*fctr; + } + } + } + + + // The left m_nel columns of A are now upper-diagonal. + // Now reduce it to diagonal form by back-solving + for (m = nRows-1; m > 0; m--) { + for (n = m-1; n>= 0; n--) { + if (m_A(n,m) != 0.0) { + fctr = m_A(n,m); + for (k = m; k < m_nsp; k++) { + m_A(n,k) -= fctr*m_A(m,k); + } + } + } + } + + // create stoichometric coefficient matrix. + for (n = 0; n < m_nsp; n++) { + if (n < m_nel) + for (k = 0; k < m_nsp - m_nel; k++) + m_N(n, k) = -m_A(n, k + m_nel); + else { + for (k = 0; k < m_nsp - m_nel; k++) m_N(n, k) = 0.0; + m_N(n, n - m_nel) = 1.0; + } + } + + // find reactions involving solution phase species + for (j = 0; j < m_nsp - m_nel; j++) { + m_solnrxn[j] = false; + for (k = 0; k < m_nsp; k++) { + if (m_N(k, j) != 0) + if (m_mix->solutionSpecies(m_species[m_order[k]])) + m_solnrxn[j] = true; + } + } + //cout << "exit: " << m_order << endl; + //for (j = 0; j < m_nsp - m_nel; j++) { + // cout << reactionString(j) << endl; + //} + } + + + /// Re-arrange a vector of species properties in sequential form + /// into sorted (components first) form. + void MultiPhaseEquil::sort(vector_fp& x) { + copy(x.begin(), x.end(), m_work2.begin()); + index_t k; + for (k = 0; k < m_nsp; k++) { + x[k] = m_work2[m_order[k]]; + } + } + + /// Re-arrange a vector of species properties in sorted form + /// (components first) into unsorted, sequential form. + void MultiPhaseEquil::unsort(vector_fp& x) { + copy(x.begin(), x.end(), m_work2.begin()); + index_t k; + for (k = 0; k < m_nsp; k++) { + x[m_order[k]] = m_work2[k]; + } + } + + /// Return a string specifying the jth reaction. + string MultiPhaseEquil::reactionString(index_t j) { + string sr = "", sp = ""; + index_t i, k; + bool rstrt = true; + bool pstrt = true; + doublereal nu; + for (i = 0; i < m_nsp; i++) { + nu = m_N(i, j); + k = m_species[m_order[i]]; + if (nu < 0.0) { + sr += coeffString(rstrt, nu, m_mix->speciesName(k)); + rstrt = false; + } + if (nu > 0.0) { + sp += coeffString(pstrt, nu, m_mix->speciesName(k)); + pstrt = false; + } + } + return sr + " <=> " + sp; + } + + doublereal MultiPhaseEquil::step(doublereal omega, vector_fp& deltaN) { + index_t k, ik; + if (omega < 0.0) + throw CanteraError("step","negative omega"); + //cout << "entering step " << m_moles << endl << deltaN << endl; + for (ik = 0; ik < m_nel; ik++) { + k = m_order[ik]; + m_lastmoles[k] = m_moles[k]; + m_moles[k] += omega * deltaN[k]; + } + + for (ik = m_nel; ik < m_nsp; ik++) { + k = m_order[ik]; + m_lastmoles[k] = m_moles[k]; + if (m_majorsp[k]) { + m_moles[k] += omega * deltaN[k]; + } + else { + m_moles[k] = fabs(m_moles[k])*fmin(10.0, exp(-m_deltaG_RT[ik - m_nel])); + } + } + setMoles(); + } + + /// Take one step in composition, given the gradient of G at the + /// starting point, and a vector of reaction steps dxi. + doublereal MultiPhaseEquil:: + stepComposition() { + + m_iter++; + index_t m, ip, ik, nsp, j, k = 0; + doublereal grad0 = computeReactionSteps(m_dxi); + + if (grad0 > 0.0) + throw CanteraError("stepComposition", "positive gradient!"); + + + // compute mole the fraction changes. + //multiply(m_N, dxi.begin(), m_work.begin()); + for (ik = 0; ik < m_nsp; ik++) { + m_work[ik] = 0.0; + k = m_order[ik]; + for (j = 0; j < m_nsp - m_nel; j++) { + m_work[ik] += m_N(ik, j) * m_dxi[j]; + } + } + + // change to sequential form + unsort(m_work); + + // scale omega to keep the major species non-negative + const doublereal FCTR = 0.99; + const doublereal MAJOR_THRESHOLD = 1.0e-12; + + doublereal omega = 1.0, omax, omegamax = 1.0; + for (ik = 0; ik < m_nsp; ik++) { + k = m_order[ik]; + + // if species k is in a multi-species solution phase, then its + // mole number must remain positive, unless the entire phase + // goes away. First we'll determine an upper bound on omega, + // such that all + if (m_dsoln[k] == 1) { + + if ((m_moles[k] > MAJOR_THRESHOLD) || (ik < m_nel)) { + omax = m_moles[k]*FCTR/(fabs(m_work[k]) + TINY); + if (m_work[k] < 0.0 && omax < omegamax) { + omegamax = omax; + if (omegamax < 1.0e-5) { + cout << m_mix->speciesName(m_species[k]) << " results in " + << " omega = " << omegamax << endl; + //cout << m_moles[k] << " " << m_work[k] << endl; + if (ik < m_nel) cout << "component" << endl; + } + } + m_majorsp[k] = true; + } + else { + m_majorsp[k] = false; + } + } + else { + if (m_work[k] < 0.0 && m_moles[k] > 0.0) { + omax = -m_moles[k]/m_work[k]; + if (omax < omegamax) { + omegamax = omax*1.000001; + if (omegamax < 1.0e-5) { + cout << m_mix->speciesName(m_species[k]) << " results in " + << " omega = " << omegamax << endl; + //cout << m_moles[k] << " " << m_work[k] << endl; + if (ik < m_nel) cout << "component" << endl; + } + } + } + m_majorsp[k] = true; + } + } + + // now take a step with this scaled omega + step(omegamax, m_work); + + // compute the gradient of G at this new position in the + // current direction. If it is positive, then we have overshot + // the minimum. In this case, interpolate back. + m_mix->getChemPotentials(m_mu.begin()); + doublereal grad1 = 0.0; + for (k = 0; k < m_nsp; k++) { + grad1 += m_work[k] * m_mu[m_species[k]]; + } + // doublereal grad1 = dot(m_work.begin(), m_work.end(), m_work2.begin()); + + omega = omegamax; + if (grad1 > 0.0) { + omega *= -grad0 / (grad1 - grad0); + for (k = 0; k < m_nsp; k++) m_moles[k] = m_lastmoles[k]; + step(omega, m_work); + } + //cout << m_moles << endl; + //cout << "omega: " << omega << " " << m_mix->gibbs() << " " << error() << endl; + return omega; + } + + + /// Compute the change in extent of reaction for each reaction. + + doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi) { + + index_t i, j, k, ik, kc, ip; + int inu; + doublereal stoich, nmoles, csum, term1, fctr, dg_rt; + vector_fp nu; + const doublereal TINY = 1.0e-20; + doublereal grad = 0.0; + + dxi.resize(m_nsp - m_nel); + computeN(); + + m_mix->getChemPotentials(m_mu.begin()); + + for (j = 0; j < m_nsp - m_nel; j++) { + + // get stoichiometric vector + getStoichVector(j, nu); + + // compute Delta G + doublereal dg_rt = 0.0; + for (k = 0; k < m_nsp; k++) { + dg_rt += m_mu[m_species[k]] * nu[k]; + } + dg_rt /= (m_temp * GasConstant); + + m_deltaG_RT[j] = dg_rt; + fctr = 1.0; + // if this is a formation reaction for a single-component phase, + // check whether reaction should be included + ik = j + m_nel; + k = m_order[ik]; + if (!m_dsoln[k]) { + if (m_moles[k] <= 0.0 && dg_rt > 0.0) { + fctr = 0.0; + } + else { + fctr = 0.05; + } + } + else if (!m_solnrxn[j]) { + fctr = 1.0; + } + else { + + // component sum + csum = 0.0; + for (k = 0; k < m_nel; k++) { + kc = m_order[k]; + stoich = nu[kc]; + nmoles = fabs(m_mix->speciesMoles(m_species[kc])) + TINY; + csum += stoich*stoich*m_dsoln[kc]/nmoles; + } + + // noncomponent term + kc = m_order[j + m_nel]; + nmoles = fabs(m_mix->speciesMoles(m_species[kc])) + TINY; + term1 = m_dsoln[kc]/nmoles; + + // sum over solution phases + doublereal sum = 0.0, psum; + for (ip = 0; ip < m_np; ip++) { + phase_t& p = m_mix->phase(ip); + if (p.nSpecies() > 1) { + psum = 0.0; + for (k = 0; k < m_nsp; k++) { + kc = m_species[k]; + if (m_mix->speciesPhaseIndex(kc) == ip) { + stoich = nu[kc]; + psum += stoich * stoich; + } + } + sum -= psum / (fabs(m_mix->phaseMoles(ip)) + TINY); + } + } + fctr = 1.0/(term1 + csum + sum); + //cout << "fctr terms = " << term1 << " " << csum << " " << sum << endl; + } + dxi[j] = -fctr*dg_rt; + index_t m; + for (m = 0; m < m_nel; m++) { + if (m_moles[m_order[m]] <= 0.0 && (m_N(m, j)*dxi[j] < 0.0)) + dxi[j] = 0.0; + } + //cout << reactionString(j) << " " << dxi[j] << " " << fctr << " " << dg_rt << endl; + grad += dxi[j]*dg_rt; + + } + return grad*GasConstant*m_temp; + } + + void MultiPhaseEquil::computeN() { + index_t m, k, isp; + + const doublereal THRESHOLD = 0.01; + + // get the species moles + + // sort mole fractions + doublereal molesum = 0.0; + for (k = 0; k < m_nsp; k++) { + m_work[k] = m_mix->speciesMoles(m_species[k]); + m_sortindex[k] = k; + molesum += m_work[k]; + } + heapsort(m_work, m_sortindex); + + // reverse order in sort index + index_t itmp; + for (k = 0; k < m_nsp/2; k++) { + itmp = m_sortindex[m_nsp-k-1]; + m_sortindex[m_nsp-k-1] = m_sortindex[k]; + m_sortindex[k] = itmp; + } + index_t ik, ij; + bool ok; + for (m = 0; m < m_nel; m++) { + for (ik = 0; ik < m_nsp; ik++) { + k = m_sortindex[ik]; + if (m_mix->nAtoms(m_species[k],m_element[m]) > 0) break; + } + ok = false; + for (ij = 0; ij < m_nel; ij++) { + if (k == m_order[ij]) ok = true; + } + if (!ok) { + //for (ij = 0; ij < m_nel; ij++) { + // cout << m_mix->speciesName(m_order[ij]) << endl; + //} + //cout << "mismatch: " << m << " " << m_mix->elementName(m) << " " << m_mix->speciesName(k) << endl; + //cout << "sortindex = " << m_sortindex << endl; + getComponents(m_sortindex); + break; + } + // for (ij = 0; ij < m_nel; ij++) + // m_lastsort[ij] = m_sortindex[ij]; + // break; + //} + } + } + + doublereal MultiPhaseEquil::error() { + index_t j, ik, k, maxj; + doublereal err, maxerr = 0.0; + for (j = 0; j < m_nsp - m_nel; j++) { + ik = j + m_nel; + k = m_order[ik]; + if (m_dsoln[k] == 0 && m_moles[k] <= 0.0) { + if (m_deltaG_RT[j] >= 0.0) err = 0.0; + else err = 1.0; + } + else { + err = fabs(m_deltaG_RT[j]); + } + if (err > maxerr) { + maxerr = err; + maxj = j; + } + } + return maxerr; + } +} diff --git a/Cantera/src/MultiPhaseEquil.h b/Cantera/src/MultiPhaseEquil.h new file mode 100644 index 000000000..75537bc63 --- /dev/null +++ b/Cantera/src/MultiPhaseEquil.h @@ -0,0 +1,122 @@ +#ifndef CT_MULTIPHASE_EQUIL +#define CT_MULTIPHASE_EQUIL + +#include "ct_defs.h" +#include "MultiPhase.h" + +namespace Cantera { + + int _equilflag(const char* xy); + + class MultiPhaseEquil { + + public: + + typedef MultiPhase mix_t; + typedef size_t index_t; + typedef DenseMatrix matrix_t; + + MultiPhaseEquil(mix_t* mix); + + virtual ~MultiPhaseEquil() {} + + int constituent(index_t m) { + if (m < m_nel) return m_order[m]; + else return -1; + } + + void getStoichVector(index_t rxn, vector_fp& nu) { + index_t k; + nu.resize(m_nsp, 0.0); + if (rxn > m_nsp - m_nel) return; + for (k = 0; k < m_nsp; k++) { + nu[m_order[k]] = m_N(k, rxn); + } + } + + int iterations() { return m_iter; } + + doublereal equilibrate(int XY, doublereal err = 1.0e-9, + int maxsteps = 1000) { + int i; + m_iter = 0; + for (i = 0; i < maxsteps; i++) { + stepComposition(); + if (error() < err) break; + } + if (i >= maxsteps) { + throw CanteraError("MultiPhaseEquil::equilibrate", + "no convergence in " + int2str(maxsteps) + + " iterations. Error = " + fp2str(error())); + } + return error(); + } + + string reactionString(index_t j); + doublereal error(); + + protected: + + void getComponents(const vector_int& order); + int setInitialMoles(); + void computeN(); + doublereal stepComposition(); + void sort(vector_fp& x); + void unsort(vector_fp& x); + doublereal step(doublereal omega, vector_fp& deltaN); + doublereal computeReactionSteps(vector_fp& dxi); + void setMoles(); + + index_t m_nel_mix, m_nsp_mix, m_np; + index_t m_nel, m_nsp; + int m_iter; + mix_t* m_mix; + doublereal m_press, m_temp; + vector_int m_order; + matrix_t m_N, m_A; + vector_fp m_work, m_work2; + vector_fp m_moles, m_lastmoles, m_dxi; + vector_fp m_deltaG_RT, m_mu; + vector m_majorsp; + vector_int m_sortindex; + vector_int m_lastsort; + vector_int m_dsoln; + vector_int m_incl_element, m_incl_species; + vector_int m_species, m_element; + vector m_solnrxn; + }; + + //----------------------------------------------------------- + // convenience functions + //----------------------------------------------------------- + + /** + * Set a mixture to a state of chemical equilibrium. The flag 'XY' + * determines the two properties that will be held fixed in the + * calculation. + */ + inline doublereal equilibrate(MultiPhase& s, int XY, + doublereal tol = 1.0e-9, int maxsteps = 1000) { + s.init(); + MultiPhaseEquil e(&s); + if (XY == TP) + return e.equilibrate(XY, tol, maxsteps); + else { + throw CanteraError("equilibrate","only fixed T, P supported"); + return -1.0; + } + } + + /** + * Set a mixture to a state of chemical equilibrium. The flag 'XY' + * determines the two properties that will be held fixed in the + * calculation. + */ + inline doublereal equilibrate(MultiPhase& s, const char* XY, + doublereal tol = 1.0e-9, int maxsteps = 1000) { + return equilibrate(s,_equilflag(XY), tol, maxsteps); + } +} + + +#endif diff --git a/Cantera/src/phasereport.cpp b/Cantera/src/phasereport.cpp index ecb9d361c..2f23cacbe 100644 --- a/Cantera/src/phasereport.cpp +++ b/Cantera/src/phasereport.cpp @@ -65,18 +65,20 @@ namespace Cantera { int kk = th.nSpecies(); array_fp x(kk); array_fp y(kk); + array_fp mu(kk); th.getMoleFractions(x.begin()); th.getMassFractions(y.begin()); - + th.getChemPotentials(mu.begin()); + doublereal rt = GasConstant * th.temperature(); int k; - sprintf(p, "\n X Y \n"); + sprintf(p, "\n X Y Chem. Pot. / RT \n"); s += p; - sprintf(p, " ------------- ------------\n"); + sprintf(p, " ------------- ------------ ------------\n"); s += p; for (k = 0; k < kk; k++) { - sprintf(p, "%18s %12.6e %12.6e\n", - th.speciesName(k).c_str(), x[k], y[k]); + sprintf(p, "%18s %12.6g %12.6g %12.6g\n", + th.speciesName(k).c_str(), x[k], y[k], mu[k]/rt); s += p; } return s; diff --git a/config.h.in b/config.h.in index 992d5856c..345f9f94f 100755 --- a/config.h.in +++ b/config.h.in @@ -5,6 +5,14 @@ #define CT_CONFIG_H +//------------------------ Development flags ------------------// +// +// These flags turn on or off features that are still in +// development and are not yet stable. + +#define DEV_EQUIL + + //------------------------ Fortran settings -------------------// @@ -12,9 +20,9 @@ // corresponding Fortran data types on your system. The defaults // are OK for most systems -typedef double doublereal; // Fortran double precision -typedef int integer; // Fortran integer -typedef int ftnlen; // Fortran hidden string length type +typedef double doublereal; // Fortran double precision +typedef int integer; // Fortran integer +typedef int ftnlen; // Fortran hidden string length type