research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

support for multiphase equilibrium

Browse source
This commit is contained in:
Dave Goodwin 2004-12-01 22:57:23 +00:00
parent bd7c2087f5
commit 987e1ddbb0
14 changed files with 1355 additions and 12 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
Cantera/clib/src/Makefile.in
View file

@ -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)

176
Cantera/clib/src/ctmultiphase.cpp Normal file
View file

@ -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<mix_t>* Cabinet<mix_t>::__storage = 0;
inline mix_t* _mix(int i) {
return Cabinet<mix_t>::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<mix_t>::cabinet()->add(m);
}
int DLL_EXPORT mix_del(int i) {
Cabinet<mix_t>::cabinet()->del(i);
return 0;
}
int DLL_EXPORT mix_copy(int i) {
return Cabinet<mix_t>::cabinet()->newCopy(i);
}
int DLL_EXPORT mix_assign(int i, int j) {
return Cabinet<mix_t>::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;
}
}
}

28
Cantera/clib/src/ctmultiphase.h Normal file
View file

@ -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

4
Cantera/cxx/include/equilibrium.h
View file

@ -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

83
Cantera/python/Cantera/Mixture.py Normal file
View file

@ -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)

234
Cantera/python/src/ctmultiphase_methods.cpp Normal file
View file

@ -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);
}

17
Cantera/python/src/methods.h
View file

@ -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

2
Cantera/python/src/pycantera.cpp
View file

@ -24,6 +24,7 @@
#include "ctreactor.h"
#include "ctfunc.h"
#include "ctonedim.h"
#include "ctmultiphase.h"
#include <iostream>
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"

2
Cantera/src/DenseMatrix.cpp
View file

@ -118,7 +118,7 @@ namespace Cantera {
*/
void multiply(const DenseMatrix& A, const double* b, double* prod) {
ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose,
static_cast<int>(A.nRows()), static_cast<int>(A.nRows()), 1.0,
static_cast<int>(A.nRows()), static_cast<int>(A.nColumns()), 1.0,
A.begin(), static_cast<int>(A.nRows()), b, 1, 0.0, prod, 1);
}

2
Cantera/src/Makefile.in
View file

@ -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

669
Cantera/src/MultiPhaseEquil.cpp Normal file
View file

@ -0,0 +1,669 @@
#include "MultiPhaseEquil.h"
#include "MultiPhase.h"
#include "sort.h"
#include "recipes.h"
#include <math.h>
#include <iostream>
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;
}
}

122
Cantera/src/MultiPhaseEquil.h Normal file
View file

@ -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<bool> 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<bool> 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

12
Cantera/src/phasereport.cpp
View file

@ -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;

14
config.h.in
View file

@ -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