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

moved files to equil subdirectory

Browse source
This commit is contained in:
Dave Goodwin 2007-05-04 14:12:33 +00:00
parent a04e7309a2
commit 9581e8db9a
10 changed files with 5920 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

823
Cantera/src/equil/BasisOptimize.cpp Normal file
View file

@ -0,0 +1,823 @@
/**
* @file BasisOptimize.cpp
* Functions which calculation optimized basis of the
* stoichiometric coefficient matrix (see /ref equil functions)
*/
/*
* $Author$
* $Date$
* $Revision$
*/
#include "ct_defs.h"
#include "ThermoPhase.h"
#include "MultiPhase.h"
using namespace Cantera;
using namespace std;
#ifdef DEBUG_HKM
namespace Cantera {
int Cantera::BasisOptimize_print_lvl = 0;
static char sbuf[1024];
}
static void print_stringTrunc(const char *str, int space, int alignment);
#endif
static int amax(double *x, int j, int n);
static void switch_pos(vector_int &orderVector, int jr, int kspec);
static int mlequ(double *c, int idem, int n, double *b, int m);
//@{
#ifndef MIN
#define MIN(x,y) (( (x) < (y) ) ? (x) : (y))
#endif
//@}
/*
* Choose the optimum basis for the calculations. This is done by
* choosing the species with the largest mole fraction
* not currently a linear combination of the previous components.
* Then, calculate the stoichiometric coefficient matrix for that
* basis.
*
* Calculates the identity of the component species in the mechanism.
* Rearranges the solution data to put the component data at the
* front of the species list.
*
* Then, calculates SC(J,I) the formation reactions for all noncomponent
* species in the mechanism.
*
* Input
* ---------
* mphase Pointer to the multiphase object. Contains the
* species mole fractions, which are used to pick the
* current optimal species component basis.
* orderVectorElement
* Order vector for the elements. The element rows
* in the formula matrix are
* rearranged according to this vector.
* orderVectorSpecies
* Order vector for the species. The species are
* rearranged according to this formula. The first
* nCompoments of this vector contain the calculated
* species components on exit.
* doFormRxn If true, the routine calculates the formation
* reaction matrix based on the calculated
* component species. If false, this step is skipped.
*
* Output
* ---------
* usedZeroedSpecies = If true, then a species with a zero concentration
* was used as a component. The problem may be
* converged.
* formRxnMatrix
*
* Return
* --------------
* returns the number of components.
*
*
*/
int Cantera::BasisOptimize(int *usedZeroedSpecies, bool doFormRxn,
MultiPhase *mphase, vector_int & orderVectorSpecies,
vector_int & orderVectorElements,
vector_fp & formRxnMatrix) {
int j, jj, k, kk, l, i, jl, ml;
bool lindep;
std::string ename;
std::string sname;
/*
* Get the total number of elements defined in the multiphase object
*/
int ne = mphase->nElements();
/*
* Get the total number of species in the multiphase object
*/
int nspecies = mphase->nSpecies();
doublereal tmp;
doublereal const USEDBEFORE = -1;
/*
* Perhaps, initialize the element ordering
*/
if ((int) orderVectorElements.size() < ne) {
orderVectorElements.resize(ne);
for (j = 0; j < ne; j++) {
orderVectorElements[j] = j;
}
}
/*
* Perhaps, initialize the species ordering
*/
if ((int) orderVectorSpecies.size() != nspecies) {
orderVectorSpecies.resize(nspecies);
for (k = 0; k < nspecies; k++) {
orderVectorSpecies[k] = k;
}
}
#ifdef DEBUG_HKM
double molSave = 0.0;
if (BasisOptimize_print_lvl >= 1) {
writelog(" "); for(i=0; i<77; i++) writelog("-"); writelog("\n");
writelog(" --- Subroutine BASOPT called to ");
writelog("calculate the number of components and ");
writelog("evaluate the formation matrix\n");
if (BasisOptimize_print_lvl > 0) {
writelog(" ---\n");
writelog(" --- Formula Matrix used in BASOPT calculation\n");
writelog(" --- Species | Order | ");
for (j = 0; j < ne; j++) {
jj = orderVectorElements[j];
writelog(" ");
ename = mphase->elementName(jj);
print_stringTrunc(ename.c_str(), 4, 1);
sprintf(sbuf,"(%1d)", j); writelog(sbuf);
}
writelog("\n");
for (k = 0; k < nspecies; k++) {
kk = orderVectorSpecies[k];
writelog(" --- ");
sname = mphase->speciesName(kk);
print_stringTrunc(sname.c_str(), 11, 1);
sprintf(sbuf," | %4d |", k); writelog(sbuf);
for (j = 0; j < ne; j++) {
jj = orderVectorElements[j];
double num = mphase->nAtoms(kk,jj);
sprintf(sbuf,"%6.1g ", num); writelog(sbuf);
}
writelog("\n");
}
writelog(" --- \n");
}
}
#endif
/*
* Calculate the maximum value of the number of components possible
* It's equal to the minimum of the number of elements and the
* number of total species.
*/
int nComponents = MIN(ne, nspecies);
int nNonComponents = nspecies - nComponents;
/*
* Set this return variable to false
*/
*usedZeroedSpecies = false;
/*
* Create an array of mole numbers
*/
vector_fp molNum(nspecies,0.0);
mphase->getMoles(DATA_PTR(molNum));
/*
* Other workspace
*/
vector_fp sm(ne*ne, 0.0);
vector_fp ss(ne, 0.0);
vector_fp sa(ne, 0.0);
if ((int) formRxnMatrix.size() < nspecies*ne) {
formRxnMatrix.resize(nspecies*ne, 0.0);
}
#ifdef DEBUG_HKM
/*
* For debugging purposes keep an unmodified copy of the array.
*/
vector_fp molNumBase(molNum);
#endif
int jr = -1;
/*
* Top of a loop of some sort based on the index JR. JR is the
* current number of component species found.
*/
do {
++jr;
/* - Top of another loop point based on finding a linearly */
/* - independent species */
do {
/*
* Search the remaining part of the mole number vector, molNum
* for the largest remaining species. Return its identity.
* kk is the raw number. k is the orderVectorSpecies index.
*/
kk = amax(DATA_PTR(molNum), 0, nspecies);
for (j = 0; j < nspecies; j++) {
if (orderVectorSpecies[j] == kk) {
k = j;
break;
}
}
if (j == nspecies) {
throw CanteraError("BasisOptimize", "orderVectorSpecies contains an error");
}
if (molNum[kk] == 0.0) *usedZeroedSpecies = true;
/*
* If the largest molNum is negative, then we are done.
*/
if (molNum[kk] == USEDBEFORE) {
nComponents = jr;
nNonComponents = nspecies - nComponents;
goto L_END_LOOP;
}
/*
* Assign a small negative number to the component that we have
* just found, in order to take it out of further consideration.
*/
#ifdef DEBUG_HKM
molSave = molNum[kk];
#endif
molNum[kk] = USEDBEFORE;
/* *********************************************************** */
/* **** CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES ****** */
/* *********************************************************** */
/*
* Modified Gram-Schmidt Method, p. 202 Dalquist
* QR factorization of a matrix without row pivoting.
*/
jl = jr;
for (j = 0; j < ne; ++j) {
jj = orderVectorElements[j];
sm[j + jr*ne] = mphase->nAtoms(kk,jj);
}
if (jl > 0) {
/*
* Compute the coefficients of JA column of the
* the upper triangular R matrix, SS(J) = R_J_JR
* (this is slightly different than Dalquist)
* R_JA_JA = 1
*/
for (j = 0; j < jl; ++j) {
ss[j] = 0.0;
for (i = 0; i < ne; ++i) {
ss[j] += sm[i + jr*ne] * sm[i + j*ne];
}
ss[j] /= sa[j];
}
/*
* Now make the new column, (*,JR), orthogonal to the
* previous columns
*/
for (j = 0; j < jl; ++j) {
for (l = 0; l < ne; ++l) {
sm[l + jr*ne] -= ss[j] * sm[l + j*ne];
}
}
}
/*
* Find the new length of the new column in Q.
* It will be used in the denominator in future row calcs.
*/
sa[jr] = 0.0;
for (ml = 0; ml < ne; ++ml) {
tmp = sm[ml + jr*ne];
sa[jr] += tmp * tmp;
}
/* **************************************************** */
/* **** IF NORM OF NEW ROW .LT. 1E-3 REJECT ********** */
/* **************************************************** */
if (sa[jr] < 1.0e-6) lindep = true;
else lindep = false;
} while(lindep);
/* ****************************************** */
/* **** REARRANGE THE DATA ****************** */
/* ****************************************** */
if (jr != k) {
#ifdef DEBUG_HKM
if (BasisOptimize_print_lvl >= 1) {
kk = orderVectorSpecies[k];
sname = mphase->speciesName(kk);
sprintf(sbuf," --- %-12.12s", sname.c_str()); writelog(sbuf);
jj = orderVectorSpecies[jr];
ename = mphase->speciesName(jj);
sprintf(sbuf,"(%9.2g) replaces %-12.12s", molSave, ename.c_str());
writelog(sbuf);
sprintf(sbuf,"(%9.2g) as component %3d\n", molNum[jj], jr);
writelog(sbuf);
}
#endif
switch_pos(orderVectorSpecies, jr, k);
}
/* - entry point from up above */
L_END_LOOP: ;
/*
* If we haven't found enough components, go back
* and find some more. (nc -1 is used below, because
* jr is counted from 0, via the C convention.
*/
} while (jr < (nComponents-1));
if (! doFormRxn) return nComponents;
/* ****************************************************** */
/* **** EVALUATE THE STOICHIOMETRY ********************** */
/* ****************************************************** */
/*
* Formulate the matrix problem for the stoichiometric
* coefficients. CX + B = 0
* C will be an nc x nc matrix made up of the formula
* vectors for the components. Each component's formular
* vector is a column. The rows are the elements.
* n rhs's will be solved for. Thus, B is an nc x n
* matrix.
*
* BIG PROBLEM 1/21/99:
*
* This algorithm makes the assumption that the
* first nc rows of the formula matrix aren't rank deficient.
* However, this might not be the case. For example, assume
* that the first element in FormulaMatrix[] is argon. Assume that
* no species in the matrix problem actually includes argon.
* Then, the first row in sm[], below will be indentically
* zero. bleh.
* What needs to be done is to perform a rearrangement
* of the ELEMENTS -> i.e. rearrange, FormulaMatrix, sp, and gai, such
* that the first nc elements form in combination with the
* nc components create an invertible sm[]. not a small
* project, but very doable.
* An alternative would be to turn the matrix problem
* below into an ne x nc problem, and do QR elimination instead
* of Gauss-Jordon elimination.
* Note the rearrangement of elements need only be done once
* in the problem. It's actually very similar to the top of
* this program with ne being the species and nc being the
* elements!!
*/
for (k = 0; k < nComponents; ++k) {
kk = orderVectorSpecies[k];
for (j = 0; j < nComponents; ++j) {
jj = orderVectorElements[j];
sm[j + k*ne] = mphase->nAtoms(kk, jj);
}
}
for (i = 0; i < nNonComponents; ++i) {
k = nComponents + i;
kk = orderVectorSpecies[k];
for (j = 0; j < nComponents; ++j) {
jj = orderVectorElements[j];
formRxnMatrix[j + i * ne] = mphase->nAtoms(kk, jj);
}
}
/*
* Use Gauss-Jordon block elimination to calculate
* the reaction matrix
*/
j = mlequ(DATA_PTR(sm), ne, nComponents, DATA_PTR(formRxnMatrix), nNonComponents);
if (j == 1) {
writelog("ERROR: mlequ returned an error condition\n");
throw CanteraError("basopt", "mlequ returned an error condition");
}
#ifdef DEBUG_HKM
if (Cantera::BasisOptimize_print_lvl >= 1) {
writelog(" ---\n");
sprintf(sbuf," --- Number of Components = %d\n", nComponents);
writelog(sbuf);
writelog(" --- Formula Matrix:\n");
writelog(" --- Components: ");
for (k = 0; k < nComponents; k++) {
kk = orderVectorSpecies[k];
sprintf(sbuf," %3d (%3d) ", k, kk); writelog(sbuf);
}
writelog("\n --- Components Moles: ");
for (k = 0; k < nComponents; k++) {
kk = orderVectorSpecies[k];
sprintf(sbuf,"%-11.3g", molNumBase[kk]); writelog(sbuf);
}
writelog("\n --- NonComponent | Moles | ");
for (i = 0; i < nComponents; i++) {
kk = orderVectorSpecies[i];
sname = mphase->speciesName(kk);
sprintf(sbuf,"%-11.10s", sname.c_str()); writelog(sbuf);
}
writelog("\n");
for (i = 0; i < nNonComponents; i++) {
k = i + nComponents;
kk = orderVectorSpecies[k];
sprintf(sbuf," --- %3d (%3d) ", k, kk); writelog(sbuf);
sname = mphase->speciesName(kk);
sprintf(sbuf,"%-10.10s", sname.c_str()); writelog(sbuf);
sprintf(sbuf,"|%10.3g|", molNumBase[kk]); writelog(sbuf);
/*
* Print the negative of formRxnMatrix[]; it's easier to interpret.
*/
for (j = 0; j < nComponents; j++) {
sprintf(sbuf," %6.2f", - formRxnMatrix[j + i * ne]);
writelog(sbuf);
}
writelog("\n");
}
writelog(" "); for (i=0; i<77; i++) writelog("-"); writelog("\n");
}
#endif
return nComponents;
} /* basopt() ************************************************************/
#ifdef DEBUG_HKM
static void print_stringTrunc(const char *str, int space, int alignment)
/***********************************************************************
* vcs_print_stringTrunc():
*
* Print a string within a given space limit. This routine
* limits the amount of the string that will be printed to a
* maximum of "space" characters.
*
* str = String -> must be null terminated.
* space = space limit for the printing.
* alignment = 0 centered
* 1 right aligned
* 2 left aligned
***********************************************************************/
{
int i, ls=0, rs=0;
int len = strlen(str);
if ((len) >= space) {
for (i = 0; i < space; i++) {
sprintf(sbuf,"%c", str[i]); writelog(sbuf);
}
} else {
if (alignment == 1) {
ls = space - len;
} else if (alignment == 2) {
rs = space - len;
} else {
ls = (space - len) / 2;
rs = space - len - ls;
}
if (ls != 0) {
for (i = 0; i < ls; i++) writelog(" ");
}
sprintf(sbuf,"%s", str); writelog(sbuf);
if (rs != 0) {
for (i = 0; i < rs; i++) writelog(" ");
}
}
}
#endif
/*
* Finds the location of the maximum component in a double vector
* INPUT
* x(*) - Vector to search
* j <= i < n : i is the range of indecises to search in X(*)
*
* RETURN
* return index of the greatest value on X(*) searched
*/
static int amax(double *x, int j, int n) {
int i;
int largest = j;
double big = x[j];
for (i = j + 1; i < n; ++i) {
if (x[i] > big) {
largest = i;
big = x[i];
}
}
return largest;
}
static void switch_pos(vector_int &orderVector, int jr, int kspec) {
int kcurr = orderVector[jr];
orderVector[jr] = orderVector[kspec];
orderVector[kspec] = kcurr;
}
/*
* vcs_mlequ:
*
* Invert an nxn matrix and solve m rhs's
*
* Solve C X + B = 0;
*
* This routine uses Gauss elimination and is optimized for the solution
* of lots of rhs's.
* A crude form of row pivoting is used here.
*
*
* c[i+j*idem] = c_i_j = Matrix to be inverted: i = row number
* j = column number
* b[i+j*idem] = b_i_j = vectors of rhs's: i = row number
* j = column number
* (each column is a new rhs)
* n = number of rows and columns in the matrix
* m = number of rhs to be solved for
* idem = first dimension in the calling routine
* idem >= n must be true
*
* Return Value
* 1 : Matrix is singluar
* 0 : solution is OK
*
* The solution is returned in the matrix b.
*/
static int mlequ(double *c, int idem, int n, double *b, int m) {
int i, j, k, l;
double R;
/*
* Loop over the rows
* -> At the end of each loop, the only nonzero entry in the column
* will be on the diagonal. We can therfore just invert the
* diagonal at the end of the program to solve the equation system.
*/
for (i = 0; i < n; ++i) {
if (c[i + i * idem] == 0.0) {
/*
* Do a simple form of row pivoting to find a non-zero pivot
*/
for (k = i + 1; k < n; ++k) {
if (c[k + i * idem] != 0.0) goto FOUND_PIVOT;
}
#ifdef DEBUG_HKM
sprintf(sbuf,"vcs_mlequ ERROR: Encountered a zero column: %d\n", i);
writelog(sbuf);
#endif
return 1;
FOUND_PIVOT: ;
for (j = 0; j < n; ++j) c[i + j * idem] += c[k + j * idem];
for (j = 0; j < m; ++j) b[i + j * idem] += b[k + j * idem];
}
for (l = 0; l < n; ++l) {
if (l != i && c[l + i * idem] != 0.0) {
R = c[l + i * idem] / c[i + i * idem];
c[l + i * idem] = 0.0;
for (j = i+1; j < n; ++j) c[l + j * idem] -= c[i + j * idem] * R;
for (j = 0; j < m; ++j) b[l + j * idem] -= b[i + j * idem] * R;
}
}
}
/*
* The negative in the last expression is due to the form of B upon
* input
*/
for (i = 0; i < n; ++i) {
for (j = 0; j < m; ++j)
b[i + j * idem] = -b[i + j * idem] / c[i + i*idem];
}
return 0;
} /* mlequ() *************************************************************/
/*
*
* ElemRearrange:
*
* This subroutine handles the rearrangement of the constraint
* equations represented by the Formula Matrix. Rearrangement is only
* necessary when the number of components is less than the number of
* elements. For this case, some constraints can never be satisfied
* exactly, because the range space represented by the Formula
* Matrix of the components can't span the extra space. These
* constraints, which are out of the range space of the component
* Formula matrix entries, are migrated to the back of the Formula
* matrix.
*
* A prototypical example is an extra element column in
* FormulaMatrix[],
* which is identically zero. For example, let's say that argon is
* has an element column in FormulaMatrix[], but no species in the
* mechanism
* actually contains argon. Then, nc < ne. Unless the entry for
* desired elementabundance vector for Ar is zero, then this
* element abundance constraint can never be satisfied. The
* constraint vector is not in the range space of the formula
* matrix.
* Also, without perturbation
* of FormulaMatrix[], BasisOptimize[] would produce a zero pivot
* because the matrix
* would be singular (unless the argon element column was already the
* last column of FormulaMatrix[].
* This routine borrows heavily from BasisOptimize algorithm. It
* finds nc constraints which span the range space of the Component
* Formula matrix, and assigns them as the first nc components in the
* formular matrix. This guarrantees that BasisOptimize has a
* nonsingular matrix to invert.
*/
int Cantera::ElemRearrange(int nComponents, const vector_fp & elementAbundances,
MultiPhase *mphase,
vector_int & orderVectorSpecies,
vector_int & orderVectorElements) {
int j, k, l, i, jl, ml, jr, ielem, jj, kk;
bool lindep = false;
int nelements = mphase->nElements();
std::string ename;
/*
* Get the total number of species in the multiphase object
*/
int nspecies = mphase->nSpecies();
double test = -1.0E10;
#ifdef DEBUG_HKM
if (BasisOptimize_print_lvl > 0) {
writelog(" "); for(i=0; i<77; i++) writelog("-"); writelog("\n");
writelog(" --- Subroutine ElemRearrange() called to ");
writelog("check stoich. coefficent matrix\n");
writelog(" --- and to rearrange the element ordering once\n");
}
#endif
/*
* Perhaps, initialize the element ordering
*/
if ((int) orderVectorElements.size() < nelements) {
orderVectorElements.resize(nelements);
for (j = 0; j < nelements; j++) {
orderVectorElements[j] = j;
}
}
/*
* Perhaps, initialize the species ordering. However, this is
* dangerous, as this ordering is assumed to yield the
* component species for the problem
*/
if ((int) orderVectorSpecies.size() != nspecies) {
orderVectorSpecies.resize(nspecies);
for (k = 0; k < nspecies; k++) {
orderVectorSpecies[k] = k;
}
}
/*
* If the elementAbundances aren't input, just create a fake one
* based on summing the column of the stoich matrix.
* This will force elements with zero species to the
* end of the element ordering.
*/
vector_fp eAbund(nelements,0.0);
if ((int) elementAbundances.size() != nelements) {
for (j = 0; j < nelements; j++) {
eAbund[j] = 0.0;
for (k = 0; k < nspecies; k++) {
eAbund[j] += fabs(mphase->nAtoms(k, j));
}
}
} else {
copy(elementAbundances.begin(), elementAbundances.end(),
eAbund.begin());
}
vector_fp sa(nelements,0.0);
vector_fp ss(nelements,0.0);
vector_fp sm(nelements*nelements,0.0);
/*
* Top of a loop of some sort based on the index JR. JR is the
* current number independent elements found.
*/
jr = -1;
do {
++jr;
/*
* Top of another loop point based on finding a linearly
* independent element
*/
do {
/*
* Search the element vector. We first locate elements that
* are present in any amount. Then, we locate elements that
* are not present in any amount.
* Return its identity in K.
*/
k = nelements;
for (ielem = jr; ielem < nelements; ielem++) {
kk = orderVectorElements[ielem];
if (eAbund[kk] != test && eAbund[kk] > 0.0) {
k = ielem;
break;
}
}
for (ielem = jr; ielem < nelements; ielem++) {
kk = orderVectorElements[ielem];
if (eAbund[kk] != test) {
k = ielem;
break;
}
}
if (k == nelements) {
// When we are here, there is an error usually.
// We haven't found the number of elements necessary.
// This is signalled by returning jr != nComponents.
#ifdef DEBUG_HKM
if (BasisOptimize_print_lvl > 0) {
sprintf(sbuf,"Error exit: returning with nComponents = %d\n", jr);
writelog(sbuf);
}
#endif
return jr;
}
/*
* Assign a large negative number to the element that we have
* just found, in order to take it out of further consideration.
*/
eAbund[kk] = test;
/* *********************************************************** */
/* **** CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX */
/* **** LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX ****** */
/* *********************************************************** */
/*
* Modified Gram-Schmidt Method, p. 202 Dalquist
* QR factorization of a matrix without row pivoting.
*/
jl = jr;
/*
* Fill in the row for the current element, k, under consideration
* The row will contain the Formula matrix value for that element
* with respect to the vector of component species.
* (note j and k indecises are flipped compared to the previous routine)
*/
for (j = 0; j < nComponents; ++j) {
jj = orderVectorSpecies[j];
kk = orderVectorElements[k];
sm[j + jr*nComponents] = mphase->nAtoms(jj,kk);
}
if (jl > 0) {
/*
* Compute the coefficients of JA column of the
* the upper triangular R matrix, SS(J) = R_J_JR
* (this is slightly different than Dalquist)
* R_JA_JA = 1
*/
for (j = 0; j < jl; ++j) {
ss[j] = 0.0;
for (i = 0; i < nComponents; ++i) {
ss[j] += sm[i + jr*nComponents] * sm[i + j*nComponents];
}
ss[j] /= sa[j];
}
/*
* Now make the new column, (*,JR), orthogonal to the
* previous columns
*/
for (j = 0; j < jl; ++j) {
for (l = 0; l < nComponents; ++l) {
sm[l + jr*nComponents] -= ss[j] * sm[l + j*nComponents];
}
}
}
/*
* Find the new length of the new column in Q.
* It will be used in the denominator in future row calcs.
*/
sa[jr] = 0.0;
for (ml = 0; ml < nComponents; ++ml) {
double tmp = sm[ml + jr*nComponents];
sa[jr] += tmp * tmp;
}
/* **************************************************** */
/* **** IF NORM OF NEW ROW .LT. 1E-6 REJECT ********** */
/* **************************************************** */
if (sa[jr] < 1.0e-6) lindep = true;
else lindep = false;
} while(lindep);
/* ****************************************** */
/* **** REARRANGE THE DATA ****************** */
/* ****************************************** */
if (jr != k) {
#ifdef DEBUG_HKM
if (BasisOptimize_print_lvl > 0) {
kk = orderVectorElements[k];
ename = mphase->elementName(kk);
writelog(" --- ");
sprintf(sbuf,"%-2.2s", ename.c_str()); writelog(sbuf);
writelog("replaces ");
kk = orderVectorElements[jr];
ename = mphase->elementName(kk);
sprintf(sbuf,"%-2.2s", ename.c_str()); writelog(sbuf);
sprintf(sbuf," as element %3d\n", jr); writelog(sbuf);
}
#endif
switch_pos(orderVectorElements, jr, k);
}
/*
* If we haven't found enough components, go back
* and find some more. (nc -1 is used below, because
* jr is counted from 0, via the C convention.
*/
} while (jr < (nComponents-1));
return nComponents;
} /* vcs_elem_rearrange() ****************************************************/

1906
Cantera/src/equil/ChemEquil.cpp Executable file
View file

File diff suppressed because it is too large Load diff

257
Cantera/src/equil/ChemEquil.h Executable file
View file

@ -0,0 +1,257 @@
/**
* @file ChemEquil.h
*
* Chemical equilibrium.
*
* $Author$
* $Date$
* $Revision$
*
* Copyright 2001 California Institute of Technology
*
*/
#ifndef CT_CHEM_EQUIL_H
#define CT_CHEM_EQUIL_H
// Cantera includes
#include "ct_defs.h"
#include "vec_functions.h"
#include "ctexceptions.h"
#include "ThermoPhase.h"
#include "DenseMatrix.h"
#include "MultiPhaseEquil.h"
namespace Cantera {
int _equilflag(const char* xy);
/**
* Chemical equilibrium options. Used internally by class ChemEquil.
*/
class EquilOpt {
public:
EquilOpt() : relTolerance(1.e-8), absElemTol(1.0E-70),maxIterations(1000),
iterations(0),
maxStepSize(10.0), propertyPair(TP), contin(false) {}
doublereal relTolerance; ///< Relative tolerance
doublereal absElemTol; ///< Abs Tol in element number
int maxIterations; ///< Maximum number of iterations
int iterations; ///< Iteration counter
/**
* Maximum step size. Largest change in any element potential or
* in log(T) allowed in one Newton step. Default: 10.0
*/
doublereal maxStepSize;
/**
* Property pair flag. Determines which two thermodynamic properties
* are fixed.
*/
int propertyPair;
/**
* Continuation flag. Set true if the calculation should be
* initialized from the last calculation. Otherwise, the
* calculation will be started from scratch and the initial
* composition and element potentials estimated.
*/
bool contin;
};
template<class M>
class PropertyCalculator;
/**
* @defgroup equil Chemical Equilibrium
*
*/
/**
* Class ChemEquil implements a chemical equilibrium solver for
* single-phase solutions. It is a "non-stoichiometric" solver in
* the terminology of Smith and Missen, meaning that every
* intermediate state is a valid chemical equilibrium state, but
* does not necessarily satisfy the element constraints. In
* contrast, the solver implemented in class MultiPhaseEquil uses
* a "stoichiometric" algorithm, in which each intermediate state
* satisfies the element constraints but is not a state of
* chemical equilibrium. Non-stoichiometric methods are faster
* when they converge, but stoichiometric ones tend to be more
* robust and can be used also for problems with multiple
* condensed phases. As expected, the ChemEquil solver is faster
* than MultiPhaseEquil for many single-phase equilibrium
* problems (particularly if there are only a few elements but
* vvery many species), but can be less stable. Problem
* situations include low temperatures where only a few species
* have non-zero mole fractions, precisely stoichiometric
* compositions (e.g. 2 H2 + O2). In general, if speed is
* important, this solver should be tried first, and if it fails
* then use MultiPhaseEquil.
* @ingroup equil
*/
class ChemEquil {
public:
//! Default Constructor
ChemEquil();
//! Constructor combined with the initialization function
/*!
* This constructor initializes the ChemEquil object with everything it
* needs to start solving equilibrium problems.
* @param s ThermoPhase object that will be used in the equilibrium calls.
*/
ChemEquil(thermo_t& s);
virtual ~ChemEquil();
int equilibrate(thermo_t& s, const char* XY,
bool useThermoPhaseElementPotentials = false);
int equilibrate(thermo_t& s, const char* XY, vector_fp& elMoles,
bool useThermoPhaseElementPotentials = false);
const vector_fp& elementPotentials() const { return m_lambda; }
/**
* Options controlling how the calculation is carried out.
* @see EquilOptions
*/
EquilOpt options;
protected:
//! Pointer to the %ThermoPhase object used to initialize this object.
/*!
* This %ThermoPhase object must be compatible with the %ThermoPhase
* objects input from the equilibrate function. Currently, this
* means that the 2 %ThermoPhases have to have consist of the same
* species and elements.
*/
thermo_t* m_phase;
/// number of atoms of element m in species k.
doublereal nAtoms(int k, int m) const { return m_comp[k*m_mm + m]; }
void initialize(thermo_t& s);
void setToEquilState(thermo_t& s,
const vector_fp& x, doublereal t);
int setInitialMoles(thermo_t& s, vector_fp& elMoleGoal);
int estimateElementPotentials(thermo_t& s, vector_fp& lambda,
vector_fp& elMolesGoal);
int estimateEP_Brinkley(thermo_t&s, vector_fp& lambda, vector_fp& elMoles);
int dampStep(thermo_t& s, vector_fp& oldx,
double oldf, vector_fp& grad, vector_fp& step, vector_fp& x,
double& f, vector_fp& elmols, double xval, double yval );
void equilResidual(thermo_t& s, const vector_fp& x,
const vector_fp& elmtotal, vector_fp& resid,
double xval, double yval);
void equilJacobian(thermo_t& s, vector_fp& x,
const vector_fp& elmols, DenseMatrix& jac,
double xval, double yval);
void adjustEloc(thermo_t& s, vector_fp & elMolesGoal);
void update(const thermo_t& s);
double calcEmoles(thermo_t& s, vector_fp& x,
const double & n_t, const vector_fp & Xmol_i_calc,
vector_fp& eMolesCalc, vector_fp& n_i_calc,
double pressureConst);
int m_mm;
int m_kk;
int m_skip;
/**
* This is equal to the rank of the stoichiometric coefficient
* matrix when it is computed. It's initialized to m_mm.
*/
int m_nComponents;
PropertyCalculator<thermo_t> *m_p1, *m_p2;
/**
* Current value of the mole fractions in the single phase.
* -> length = m_kk.
*/
vector_fp m_molefractions;
/**
* Current value of the dimensional element potentials
* -> length = m_mm
*/
vector_fp m_lambda;
/*
* Current value of the sum of the element abundances given the
* current element potentials.
*/
doublereal m_elementTotalSum;
/*
* Current value of the element mole fractions. Note these aren't
* the goal element mole fractions.
*/
vector_fp m_elementmolefracs;
vector_fp m_reswork;
vector_fp m_jwork1;
vector_fp m_jwork2;
/*
* Storage of the element compositions
* natom(k,m) = m_comp[k*m_mm+ m];
*/
vector_fp m_comp;
doublereal m_temp, m_dens;
doublereal m_p0;
/**
* Index of the element id corresponding to the electric charge of each
* species. Equal to -1 if there is no such element id.
*/
int m_eloc;
doublereal m_startTemp, m_startDens;
vector_fp m_startSoln;
vector_fp m_grt;
vector_fp m_mu_RT;
/**
* Dimensionless values of the gibbs free energy for the
* standard state of each species, at the temperature and
* pressure of the solution (the star standard state).
*/
vector_fp m_muSS_RT;
vector_int m_component;
/*
* element fractional cutoff, below which the element will be
* zeroed.
*/
double m_elemFracCutoff;
bool m_doResPerturb;
vector_int m_orderVectorElements;
vector_int m_orderVectorSpecies;
};
#ifdef DEBUG_HKM
extern int ChemEquil_print_lvl;
#endif
}
#endif

879
Cantera/src/equil/MultiPhase.cpp Normal file
View file

@ -0,0 +1,879 @@
/**
* @file MultiPhase.cpp
* Definitions for the \link Cantera::MultiPhase MultiPhase\endlink
* object that is used to set up multiphase equilibrium problems (see \ref equilfunctions).
*/
/*
* $Author$
* $Date$
* $Revision$
*/
#include "MultiPhase.h"
#include "MultiPhaseEquil.h"
#include "ThermoPhase.h"
#include "DenseMatrix.h"
#include "stringUtils.h"
using namespace std;
namespace Cantera {
/// Constructor.
MultiPhase::MultiPhase() : m_temp(0.0), m_press(0.0),
m_nel(0), m_nsp(0), m_init(false), m_eloc(-1),
m_Tmin(1.0), m_Tmax(100000.0) {
}
void MultiPhase::
addPhases(MultiPhase& mix) {
index_t n;
for (n = 0; n < mix.m_np; n++) {
addPhase(mix.m_phase[n], mix.m_moles[n]);
}
}
void MultiPhase::
addPhases(phase_list& phases, const vector_fp& phaseMoles) {
index_t np = phases.size();
index_t n;
for (n = 0; n < np; n++) {
addPhase(phases[n], phaseMoles[n]);
}
init();
}
void MultiPhase::
addPhase(phase_t* p, doublereal moles) {
if (m_init) {
throw CanteraError("addPhase",
"phases cannot be added after init() has been called.");
}
// save the pointer to the phase object
m_phase.push_back(p);
// store its number of moles
m_moles.push_back(moles);
m_temp_OK.push_back(true);
// update the number of phases and the total number of
// species
m_np = m_phase.size();
m_nsp += p->nSpecies();
// determine if this phase has new elements
// for each new element, add an entry in the map
// from names to index number + 1:
string ename;
// iterate over the elements in this phase
index_t m, nel = p->nElements();
for (m = 0; m < nel; m++) {
ename = p->elementName(m);
// if no entry is found for this element name, then
// it is a new element. In this case, add the name
// to the list of names, increment the element count,
// and add an entry to the name->(index+1) map.
if (m_enamemap.find(ename) == m_enamemap.end()) {
m_enamemap[ename] = m_nel + 1;
m_enames.push_back(ename);
m_atomicNumber.push_back(p->atomicNumber(m));
// Element 'E' (or 'e') is special. Note its location.
if (ename == "E" || ename == "e") m_eloc = m_nel;
m_nel++;
}
}
// If the mixture temperature hasn't been set, then set the
// temperature and pressure to the values for the phase being
// added.
if (m_temp == 0.0 && p->temperature() > 0.0) {
m_temp = p->temperature();
m_press = p->pressure();
}
// If this is a solution phase, update the minimum and maximum
// mixture temperatures. Stoichiometric phases are excluded,
// since a mixture may define multiple stoichiometric phases,
// each of which has thermo data valid only over a limited
// range. For example, a mixture might be defined to contain a
// phase representing water ice and one representing liquid
// water, only one of which should be present if the mixture
// represents an equilibrium state.
if (p->nSpecies() > 1) {
double t = p->minTemp();
if (t > m_Tmin) m_Tmin = t;
t = p->maxTemp();
if (t < m_Tmax) m_Tmax = t;
}
}
// Process phases and build atomic composition array. This method
// must be called after all phases are added, before doing
// anything else with the mixture. After init() has been called,
// no more phases may be added.
void MultiPhase::init() {
if (m_init) return;
index_t ip, kp, k = 0, nsp, m;
int mlocal;
string sym;
// allocate space for the atomic composition matrix
m_atoms.resize(m_nel, m_nsp, 0.0);
m_moleFractions.resize(m_nsp, 0.0);
m_elemAbundances.resize(m_nel, 0.0);
// iterate over the elements
// -> fill in m_atoms(m,k), m_snames(k), m_spphase(k),
// m_sptart(ip)
for (m = 0; m < m_nel; m++) {
sym = m_enames[m];
k = 0;
// iterate over the phases
for (ip = 0; ip < m_np; ip++) {
phase_t* p = m_phase[ip];
nsp = p->nSpecies();
mlocal = p->elementIndex(sym);
for (kp = 0; kp < nsp; kp++) {
if (mlocal >= 0) {
m_atoms(m, k) = p->nAtoms(kp, mlocal);
}
if (m == 0) {
m_snames.push_back(p->speciesName(kp));
if (kp == 0) {
m_spstart.push_back(m_spphase.size());
}
m_spphase.push_back(ip);
}
k++;
}
}
}
if (m_eloc >= 0) {
doublereal esum;
for (k = 0; k < m_nsp; k++) {
esum = 0.0;
for (m = 0; m < m_nel; m++) {
if (int(m) != m_eloc)
esum += m_atoms(m,k) * m_atomicNumber[m];
}
//m_atoms(m_eloc, k) += esum;
}
}
/// set the initial composition within each phase to the
/// mole fractions stored in the phase objects
m_init = true;
updateMoleFractions();
}
// Return a reference to phase n. The state of phase n is
// also updated to match the state stored locally in the
// mixture object.
MultiPhase::phase_t& MultiPhase::phase(index_t n) {
if (!m_init) init();
m_phase[n]->setState_TPX(m_temp, m_press,
DATA_PTR(m_moleFractions) + m_spstart[n]);
return *m_phase[n];
}
/// Moles of species \c k.
doublereal MultiPhase::speciesMoles(index_t k) const {
index_t ip = m_spphase[k];
return m_moles[ip]*m_moleFractions[k];
}
/// Total moles of element m, summed over all
/// phases
doublereal MultiPhase::elementMoles(index_t m) const {
doublereal sum = 0.0, phasesum;
index_t i, k = 0, ik, nsp;
for (i = 0; i < m_np; i++) {
phasesum = 0.0;
nsp = m_phase[i]->nSpecies();
for (ik = 0; ik < nsp; ik++) {
k = speciesIndex(ik, i);
phasesum += m_atoms(m,k)*m_moleFractions[k];
}
sum += phasesum * m_moles[i];
}
return sum;
}
/// Total charge, summed over all phases
doublereal MultiPhase::charge() const {
doublereal sum = 0.0;
index_t i;
for (i = 0; i < m_np; i++) {
sum += phaseCharge(i);
}
return sum;
}
/// Net charge of one phase (Coulombs). The net charge is computed as
/// \f[ Q_p = N_p \sum_k F z_k X_k \f]
/// where the sum runs only over species in phase \a p.
/// @param p index of the phase for which the charge is desired.
doublereal MultiPhase::phaseCharge(index_t p) const {
doublereal phasesum = 0.0;
int ik, k, nsp = m_phase[p]->nSpecies();
for (ik = 0; ik < nsp; ik++) {
k = speciesIndex(ik, p);
phasesum += m_phase[p]->charge(ik)*m_moleFractions[k];
}
return Faraday*phasesum*m_moles[p];
}
/// Get the chemical potentials of all species in all phases.
void MultiPhase::getChemPotentials(doublereal* mu) const {
index_t i, loc = 0;
updatePhases();
for (i = 0; i < m_np; i++) {
m_phase[i]->getChemPotentials(mu + loc);
loc += m_phase[i]->nSpecies();
}
}
// Get chemical potentials of species with valid thermo
// data. This method is designed for use in computing chemical
// equilibrium by Gibbs minimization. For solution phases (more
// than one species), this does the same thing as
// getChemPotentials. But for stoichiometric phases, this writes
// into array \a mu the user-specified value \a not_mu instead of
// the chemical potential if the temperature is outside the range
// for which the thermo data for the one species in the phase are
// valid. The need for this arises since many condensed phases
// have thermo data fit only for the temperature range for which
// they are stable. For example, in the NASA database, the fits
// for H2O(s) are only done up to 0 C, the fits for H2O(L) are
// only done from 0 C to 100 C, etc. Using the polynomial fits outside
// the range for which the fits were done can result in spurious
// chemical potentials, and can lead to condensed phases
// appearing when in fact they should be absent.
//
// By setting \a not_mu to a large positive value, it is possible
// to force routines which seek to minimize the Gibbs free energy
// of the mixture to zero out any phases outside the temperature
// range for which their thermo data are valid.
//
// If this method is called with \a standard set to true, then
// the composition-independent standard chemical potentials are
// returned instead of the composition-dependent chemical
// potentials.
//
void MultiPhase::getValidChemPotentials(doublereal not_mu,
doublereal* mu, bool standard) const {
index_t i, loc = 0;
updatePhases();
// iterate over the phases
for (i = 0; i < m_np; i++) {
if (tempOK(i) || m_phase[i]->nSpecies() > 1) {
if (!standard)
m_phase[i]->getChemPotentials(mu + loc);
else
m_phase[i]->getStandardChemPotentials(mu + loc);
}
else
fill(mu + loc, mu + loc + m_phase[i]->nSpecies(), not_mu);
loc += m_phase[i]->nSpecies();
}
}
/// True if species \a k belongs to a solution phase.
bool MultiPhase::solutionSpecies(index_t k) const {
if (m_phase[m_spphase[k]]->nSpecies() > 1)
return true;
else
return false;
}
/// The Gibbs free energy of the mixture (J).
doublereal MultiPhase::gibbs() const {
index_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < m_np; i++)
sum += m_phase[i]->gibbs_mole() * m_moles[i];
return sum;
}
/// The enthalpy of the mixture (J).
doublereal MultiPhase::enthalpy() const {
index_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < m_np; i++)
sum += m_phase[i]->enthalpy_mole() * m_moles[i];
return sum;
}
/// The entropy of the mixture (J/K).
doublereal MultiPhase::entropy() const {
index_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < m_np; i++)
sum += m_phase[i]->entropy_mole() * m_moles[i];
return sum;
}
/// The specific heat at constant pressure and composition (J/K).
/// Note that this does not account for changes in composition of
/// the mixture with temperature.
doublereal MultiPhase::cp() const {
index_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < m_np; i++)
sum += m_phase[i]->cp_mole() * m_moles[i];
return sum;
}
/// Set the mole fractions of phase \a n to the values in
/// array \a x.
void MultiPhase::setPhaseMoleFractions(index_t n, doublereal* x) {
phase_t* p = m_phase[n];
p->setState_TPX(m_temp, m_press, x);
}
// Set the species moles using a map. The map \a xMap maps
// species name strings to mole numbers. Mole numbers that are
// less than or equal to zero will be set to zero.
void MultiPhase::setMolesByName(compositionMap& xMap) {
int kk = nSpecies();
doublereal x;
vector_fp moles(kk, 0.0);
for (int k = 0; k < kk; k++) {
x = xMap[speciesName(k)];
if (x > 0.0) moles[k] = x;
}
setMoles(DATA_PTR(moles));
}
// Set the species moles using a string. Unspecified species are
// set to zero.
void MultiPhase::setMolesByName(const std::string& x) {
compositionMap xx;
// add an entry in the map for every species, with value -1.0.
// Function parseCompString (stringUtils.cpp) uses the names
// in the map to specify the allowed species.
int kk = nSpecies();
for (int k = 0; k < kk; k++) {
xx[speciesName(k)] = -1.0;
}
// build the composition map from the string, and then set the
// moles.
parseCompString(x, xx);
setMolesByName(xx);
}
// Get the mole numbers of all species in the multiphase
// object
void MultiPhase::getMoles(doublereal * molNum) const {
/*
* First copy in the mole fractions
*/
copy(m_moleFractions.begin(), m_moleFractions.end(), molNum);
index_t ik;
doublereal *dtmp = molNum;
for (index_t ip = 0; ip < m_np; ip++) {
doublereal phasemoles = m_moles[ip];
phase_t* p = m_phase[ip];
index_t nsp = p->nSpecies();
for (ik = 0; ik < nsp; ik++) {
*(dtmp++) *= phasemoles;
}
}
}
/// Set the species moles to the values in array \a n. The state
/// of each phase object is also updated to have the specified
/// composition and the mixture temperature and pressure.
void MultiPhase::setMoles(doublereal* n) {
if (!m_init) init();
index_t ip, loc = 0;
index_t ik, k = 0, nsp;
doublereal phasemoles;
for (ip = 0; ip < m_np; ip++) {
phase_t* p = m_phase[ip];
nsp = p->nSpecies();
phasemoles = 0.0;
for (ik = 0; ik < nsp; ik++) {
phasemoles += n[k];
k++;
}
m_moles[ip] = phasemoles;
if (nsp > 1) {
p->setState_TPX(m_temp, m_press, n + loc);
p->getMoleFractions(DATA_PTR(m_moleFractions) + loc);
}
else {
m_moleFractions[loc] = 1.0;
}
loc += nsp;
}
}
void MultiPhase::getElemAbundances(doublereal *elemAbundances) const {
index_t eGlobal;
calcElemAbundances();
for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
elemAbundances[eGlobal] = m_elemAbundances[eGlobal];
}
}
// Internal routine to calculate the element abundance vector
void MultiPhase::calcElemAbundances() const {
index_t loc = 0;
index_t eGlobal;
int ik, kGlobal;
doublereal spMoles;
for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
m_elemAbundances[eGlobal] = 0.0;
}
for (index_t ip = 0; ip < m_np; ip++) {
phase_t* p = m_phase[ip];
int nspPhase = p->nSpecies();
doublereal phasemoles = m_moles[ip];
for (ik = 0; ik < nspPhase; ik++) {
kGlobal = loc + ik;
spMoles = m_moleFractions[kGlobal] * phasemoles;
for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
m_elemAbundances[eGlobal] += m_atoms(eGlobal, kGlobal) * spMoles;
}
}
loc += nspPhase;
}
}
/// The total mixture volume [m^3].
doublereal MultiPhase::volume() const {
int i;
doublereal sum = 0;
for (i = 0; i < int(m_np); i++) {
sum += m_moles[i]/m_phase[i]->molarDensity();
}
return sum;
}
doublereal MultiPhase::equilibrate(int XY, doublereal err,
int maxsteps, int maxiter, int loglevel) {
doublereal error;
bool strt = false;
doublereal dt;
doublereal h0;
int n;
bool start;
doublereal ferr, hnow, herr = 1.0;
doublereal snow, serr = 1.0, s0;
doublereal Tlow = -1.0, Thigh = -1.0;
doublereal Hlow = Undef, Hhigh = Undef, tnew;
doublereal dta=0.0, dtmax, cpb;
MultiPhaseEquil* e = 0;
if (!m_init) init();
beginLogGroup("MultiPhase::equilibrate", loglevel);
if (XY == TP) {
addLogEntry("problem type","fixed T,P");
addLogEntry("Temperature",temperature());
addLogEntry("Pressure", pressure());
// create an equilibrium manager
e = new MultiPhaseEquil(this);
try {
error = e->equilibrate(XY, err, maxsteps);
}
catch (CanteraError err) {
endLogGroup();
delete e;
e = 0;
throw err;
}
goto done;
}
else if (XY == HP) {
h0 = enthalpy();
Tlow = 0.5*m_Tmin; // lower bound on T
Thigh = 2.0*m_Tmax; // upper bound on T
addLogEntry("problem type","fixed H,P");
addLogEntry("H target",fp2str(h0));
for (n = 0; n < maxiter; n++) {
// if 'strt' is false, the current composition will be used as
// the starting estimate; otherwise it will be estimated
// if (e) {
// cout << "e should be zero, but it is not!" << endl;
// delete e;
// }
e = new MultiPhaseEquil(this, strt);
// start with a loose error tolerance, but tighten it as we get
// close to the final temperature
beginLogGroup("iteration "+int2str(n));
try {
error = e->equilibrate(TP, err, maxsteps);
hnow = enthalpy();
// the equilibrium enthalpy monotonically increases with T;
// if the current value is below the target, the we know the
// current temperature is too low. Set
if (hnow < h0) {
if (m_temp > Tlow) {
Tlow = m_temp;
Hlow = hnow;
}
}
// the current enthalpy is greater than the target; therefore the
// current temperature is too high.
else {
if (m_temp < Thigh) {
Thigh = m_temp;
Hhigh = hnow;
}
}
if (Hlow != Undef && Hhigh != Undef) {
cpb = (Hhigh - Hlow)/(Thigh - Tlow);
dt = (h0 - hnow)/cpb;
dta = fabs(dt);
dtmax = 0.5*fabs(Thigh - Tlow);
if (dta > dtmax) dt *= dtmax/dta;
}
else {
tnew = sqrt(Tlow*Thigh);
dt = tnew - m_temp;
//cpb = cp();
}
herr = fabs((h0 - hnow)/h0);
addLogEntry("T",fp2str(temperature()));
addLogEntry("H",fp2str(hnow));
addLogEntry("H rel error",fp2str(herr));
addLogEntry("lower T bound",fp2str(Tlow));
addLogEntry("upper T bound",fp2str(Thigh));
endLogGroup(); // iteration
if (herr < err) { // || dta < 1.0e-4) {
addLogEntry("T iterations",int2str(n));
addLogEntry("Final T",fp2str(temperature()));
addLogEntry("H rel error",fp2str(herr));
goto done;
}
tnew = m_temp + dt;
if (tnew < 0.0) tnew = 0.5*m_temp;
//dta = fabs(tnew - m_temp);
setTemperature(tnew);
// if the size of Delta T is not too large, use
// the current composition as the starting estimate
if (dta < 100.0) strt = false;
}
catch (CanteraError err) {
if (!strt) {
addLogEntry("no convergence",
"try estimating starting composition");
strt = true;
}
else {
tnew = 0.5*(m_temp + Thigh);
if (fabs(tnew - m_temp) < 1.0) tnew = m_temp + 1.0;
setTemperature(tnew);
addLogEntry("no convergence",
"trying T = "+fp2str(m_temp));
}
endLogGroup();
}
delete e;
e = 0;
}
addLogEntry("reached max number of T iterations",int2str(maxiter));
endLogGroup();
throw CanteraError("MultiPhase::equilibrate",
"No convergence for T");
}
else if (XY == SP) {
s0 = entropy();
start = true;
Tlow = 1.0; // m_Tmin; // lower bound on T
Thigh = 1.0e6; // m_Tmax; // upper bound on T
addLogEntry("problem type","fixed S,P");
addLogEntry("S target",fp2str(s0));
addLogEntry("min T",fp2str(Tlow));
addLogEntry("max T",fp2str(Thigh));
for (n = 0; n < maxiter; n++) {
if (e) delete e;
e = new MultiPhaseEquil(this, strt);
ferr = 0.1;
if (fabs(dt) < 1.0) ferr = err;
//start = false;
beginLogGroup("iteration "+int2str(n));
try {
error = e->equilibrate(TP, err, maxsteps);
snow = entropy();
if (snow < s0) {
if (m_temp > Tlow) Tlow = m_temp;
}
else {
if (m_temp < Thigh) Thigh = m_temp;
}
serr = fabs((s0 - snow)/s0);
addLogEntry("T",fp2str(temperature()));
addLogEntry("S",fp2str(snow));
addLogEntry("S rel error",fp2str(serr));
endLogGroup();
dt = (s0 - snow)*m_temp/cp();
dtmax = 0.5*fabs(Thigh - Tlow);
dtmax = (dtmax > 500.0 ? 500.0 : dtmax);
dta = fabs(dt);
if (dta > dtmax) dt *= dtmax/dta;
if (herr < err || dta < 1.0e-4) {
addLogEntry("T iterations",int2str(n));
addLogEntry("Final T",fp2str(temperature()));
addLogEntry("S rel error",fp2str(serr));
goto done;
}
tnew = m_temp + dt;
setTemperature(tnew);
// if the size of Delta T is not too large, use
// the current composition as the starting estimate
if (dta < 100.0) strt = false;
}
catch (CanteraError err) {
if (!strt) {
addLogEntry("no convergence",
"setting strt to True");
strt = true;
}
else {
tnew = 0.5*(m_temp + Thigh);
setTemperature(tnew);
addLogEntry("no convergence",
"trying T = "+fp2str(m_temp));
}
endLogGroup();
}
delete e;
e = 0;
}
addLogEntry("reached max number of T iterations",int2str(maxiter));
endLogGroup();
throw CanteraError("MultiPhase::equilibrate",
"No convergence for T");
}
// else if (XY == SP) {
// if (loglevel > 0) {
// addLogEntry("problem type","fixed S,P");
// }
// doublereal dt = 1.0e3;
// doublereal s0 = entropy();
// int n;
// bool start = true;
// doublereal ferr, snow, serr, tnew;
// for (n = 0; n < maxiter; n++) {
// e = new MultiPhaseEquil(this, start);
// ferr = 0.1;
// start = false;
// if (fabs(dt) < 1.0) ferr = err;
// if (loglevel > 1) {
// beginLogGroup("iteration "+int2str(n));
// }
// try {
// error = e->equilibrate(TP, ferr, maxsteps, loglevel-1);
// snow = entropy();
// tnew = exp(0.5*(s0 - snow)/cp())*temperature();
// serr = fabs((s0 - snow)/s0);
// if (loglevel > 1) {
// addLogEntry("T",fp2str(temperature()));
// addLogEntry("S rel error",fp2str(serr));
// endLogGroup();
// }
// if (serr < err) {
// if (loglevel > 0) {
// addLogEntry("T iterations",int2str(n));
// addLogEntry("Final T",fp2str(temperature()));
// addLogEntry("S rel error",fp2str(serr));
// }
// goto done;
// }
// setTemperature(tnew);
// }
// catch (CanteraError err) {
// delete e;
// if (!strt) {
// if (loglevel > 0)
// addLogEntry("no convergence",
// "setting strt to True");
// strt = true;
// }
// else {
// tnew = 0.5*(m_temp + Thigh);
// setTemperature(tnew);
// if (loglevel > 0)
// addLogEntry("no convergence",
// "trying T = "+fp2str(m_temp));
// }
// }
// endLogGroup();
// }
// if (loglevel > 0) write_logfile("equil_err.html");
// throw CanteraError("MultiPhase::equilibrate",
// "No convergence for T");
// }
else if (XY == TV) {
addLogEntry("problem type","fixed T, V");
// doublereal dt = 1.0e3;
doublereal v0 = volume();
doublereal dVdP;
int n;
bool start = true;
doublereal error, vnow, pnow, verr;
for (n = 0; n < maxiter; n++) {
pnow = pressure();
MultiPhaseEquil e(this, start);
start = false;
beginLogGroup("iteration "+int2str(n));
error = e.equilibrate(TP, err, maxsteps);
vnow = volume();
verr = fabs((v0 - vnow)/v0);
addLogEntry("P",fp2str(pressure()));
addLogEntry("V rel error",fp2str(verr));
endLogGroup();
if (verr < err) {
addLogEntry("P iterations",int2str(n));
addLogEntry("Final P",fp2str(pressure()));
addLogEntry("V rel error",fp2str(verr));
goto done;
}
// find dV/dP
setPressure(pnow*1.01);
dVdP = (volume() - vnow)/(0.01*pnow);
setPressure(pnow + 0.5*(v0 - vnow)/dVdP);
}
}
else {
endLogGroup();
throw CanteraError("MultiPhase::equilibrate","unknown option");
}
return -1.0;
done:
delete e;
e = 0;
endLogGroup();
return err;
}
#ifdef MULTIPHASE_DEVEL
void importFromXML(string infile, string id) {
XML_Node* root = get_XML_File(infile);
if (id == "-") id = "";
XML_Node* x = get_XML_Node(string("#")+id, root);
if (x.name() != "multiphase")
throw CanteraError("MultiPhase::importFromXML",
"Current XML_Node is not a multiphase element.");
vector<XML_Node*> phases;
x.getChildren("phase",phases);
int np = phases.size();
int n;
ThermoPhase* p;
for (n = 0; n < np; n++) {
XML_Node& ph = *phases[n];
srcfile = infile;
if (ph.hasAttrib("src")) srcfile = ph["src"];
idstr = ph["id"];
p = newPhase(srcfile, idstr);
if (p) {
addPhase(p, ph.value());
}
}
}
#endif
// Name of element \a m.
std::string MultiPhase::elementName(int m) const {
return m_enames[m];
}
// Index of element with name \a name.
int MultiPhase::elementIndex(std::string name) const {
for (size_t e = 0; e < m_nel; e++) {
if (m_enames[e] == name) {
return (int) e;
}
}
return -1;
}
// Name of species with global index \a k.
std::string MultiPhase::speciesName(int k) const {
return m_snames[k];
}
//-------------------------------------------------------------
//
// protected methods
//
//-------------------------------------------------------------
/// Update the locally-stored species mole fractions.
void MultiPhase::updateMoleFractions() {
index_t ip, loc = 0;
for (ip = 0; ip < m_np; ip++) {
phase_t* p = m_phase[ip];
p->getMoleFractions(DATA_PTR(m_moleFractions) + loc);
loc += p->nSpecies();
}
calcElemAbundances();
}
/// synchronize the phase objects with the mixture state. This
/// method sets each phase to the mixture temperature and
/// pressure, and sets the phase mole fractions based on the
/// mixture mole numbers.
void MultiPhase::updatePhases() const {
index_t p, nsp, loc = 0;
for (p = 0; p < m_np; p++) {
nsp = m_phase[p]->nSpecies();
const doublereal* x = DATA_PTR(m_moleFractions) + loc;
loc += nsp;
m_phase[p]->setState_TPX(m_temp, m_press, x);
m_temp_OK[p] = true;
if (m_temp < m_phase[p]->minTemp()
|| m_temp > m_phase[p]->maxTemp()) m_temp_OK[p] = false;
}
}
}

721
Cantera/src/equil/MultiPhase.h Normal file
View file

@ -0,0 +1,721 @@
/**
* @file MultiPhase.h
* Headers for the \link Cantera::MultiPhase MultiPhase\endlink
* object that is used to set up multiphase equilibrium problems (see \ref equilfunctions).
*
*/
/*
* $Author$
* $Date$
* $Revision$
*/
#ifndef CT_MULTIPHASE_H
#define CT_MULTIPHASE_H
#include "ct_defs.h"
#include "DenseMatrix.h"
#include "ThermoPhase.h"
namespace Cantera {
//! A class for multiphase mixtures. The mixture can contain any
//! number of phases of any type.
/*!
* All phases have the same
* temperature and pressure, and a specified number of moles.
* The phases do not need to have the same elements. For example,
* a mixture might consist of a gaseous phase with elements (H,
* C, O, N), a solid carbon phase containing only element C,
* etc. A master element set will be constructed for the mixture
* that is the union of the elements of each phase.
*
* Below, reference is made to global species and global elements.
* These refer to the collective species and elements encompassing
* all of the phases tracked by the object.
*
* @ingroup equilfunctions
*/
class MultiPhase {
public:
//! Shorthand for an index variable that can't be negative
typedef size_t index_t;
//! Shorthand for a ThermoPhase
typedef ThermoPhase phase_t;
//! shorthand for a 2D matrix
typedef DenseMatrix array_t;
//! Shorthand for a vector of pointers to ThermoPhase's
typedef std::vector<phase_t*> phase_list;
/// Constructor. The constructor takes no arguments, since
/// phases are added using method addPhase.
MultiPhase();
/// Destructor. Does nothing. Class MultiPhase does not take
/// "ownership" (i.e. responsibility for destroying) the
/// phase objects.
virtual ~MultiPhase() {}
//! Add a vector of phases to the mixture
/*!
* See the single addPhases command. This just does a bunch of phases
* at one time
* @param phases Vector of pointers to phases
* @param phaseMoles Vector of mole numbers in each phase (kmol)
*/
void addPhases(phase_list& phases, const vector_fp& phaseMoles);
//! Add all phases present in 'mix' to this mixture.
/*!
* @param mix Add all of the phases in another MultiPhase
* object to the current object.
*/
void addPhases(MultiPhase& mix);
//! Add a phase to the mixture.
/*!
* This function must be called befure the init() function is called,
* which serves to freeze the MultiPhase.
*
* @param p pointer to the phase object
* @param moles total number of moles of all species in this phase
*/
void addPhase(phase_t* p, doublereal moles);
/// Number of elements.
int nElements() const { return int(m_nel); }
//! Returns the string name of the global element \a m.
/*!
* @param m index of the global element
*/
std::string elementName(int m) const;
//! Returns the index of the element with name \a name.
/*!
* @param name String name of the global element
*/
int elementIndex(std::string name) const;
//! Number of species, summed over all phases.
int nSpecies() const { return int(m_nsp); }
//! Name of species with global index \a kGlob
/*!
* @param kGlob global species index
*/
std::string speciesName(int kGlob) const;
//! Returns the Number of atoms of global element \a mGlob in
//! global species \a kGlob.
/*!
* @param kGlob global species index
* @param mGlob global element index
* @return returns the number of atoms.
*/
doublereal nAtoms(int kGlob, int mGlob) {
if (!m_init) init();
return m_atoms(mGlob, kGlob);
}
/// Returns the global Species mole fractions.
/*!
* Write the array of species mole
* fractions into array \c x. The mole fractions are
* normalized to sum to one in each phase.
*
* @param x vector of mole fractions.
* Length = number of global species.
*/
void getMoleFractions(doublereal* x) const {
std::copy(m_moleFractions.begin(), m_moleFractions.end(), x);
}
//! Process phases and build atomic composition array.
/*!This method
* must be called after all phases are added, before doing
* anything else with the mixture. After init() has been called,
* no more phases may be added.
*/
void init();
//! Return the number of moles in phase n.
/*!
* @param n Index of the phase.
*/
doublereal phaseMoles(index_t n) const {
return m_moles[n];
}
//! Set the number of moles of phase with index n.
/*!
* @param n Index of the phase
* @param moles Number of moles in the phase (kmol)
*/
void setPhaseMoles(index_t n, doublereal moles) {
m_moles[n] = moles;
}
/// Return a %ThermoPhase reference to phase n.
/*! The state of phase n is
* also updated to match the state stored locally in the
* mixture object.
*
* @param n Phase Index
*
* @return Reference to the %ThermoPhase object for the phase
*/
phase_t& phase(index_t n);
//! Returns the moles of global species \c k.
/*!
* Returns the moles of global species k.
* units = kmol
*
* @param kGlob Global species index k
*/
doublereal speciesMoles(index_t kGlob) const;
//! Index of the species belonging to phase number \c p
//! with local index \c k within the phase.
/*!
* Returns the index of the global species
*
* @param k local index of the species within the phase
* @param p index of the phase
*/
int speciesIndex(index_t k, index_t p) const {
return m_spstart[p] + k;
}
/// Minimum temperature for which all solution phases have
/// valid thermo data. Stoichiometric phases are not
/// considered, since they may have thermo data only valid for
/// conditions for which they are stable.
doublereal minTemp() const { return m_Tmin; }
/// Maximum temperature for which all solution phases have
/// valid thermo data. Stoichiometric phases are not
/// considered, since they may have thermo data only valid for
/// conditions for which they are stable.
doublereal maxTemp() const { return m_Tmax; }
/// Total charge (Coulombs).
doublereal charge() const;
/// Charge (Coulombs) of phase with index \a p.
/*!
* @param p Phase Index
*/
doublereal phaseCharge(index_t p) const;
/// Total moles of global element \a m, summed over all phases.
/*!
* @param m Index of the global element
*/
doublereal elementMoles(index_t m) const;
//! Returns a vector of Chemical potentials.
/*!
* Write into array \a mu the chemical
* potentials of all species [J/kmol]. The chemical
* potentials are related to the activities by
*
* \f$
* \mu_k = \mu_k^0(T, P) + RT \ln a_k.
* \f$.
*
* @param mu Chemical potential vector.
* Length = num global species.
* Units = J/kmol.
*/
void getChemPotentials(doublereal* mu) const;
/// Returns a vector of Valid chemical potentials.
/*!
* Write into array \a mu the
* chemical potentials of all species with thermo data valid
* for the current temperature [J/kmol]. For other species,
* set the chemical potential to the value \a not_mu. If \a
* standard is set to true, then the values returned are
* standard chemical potentials.
*
* This method is designed for use in computing chemical
* equilibrium by Gibbs minimization. For solution phases (more
* than one species), this does the same thing as
* getChemPotentials. But for stoichiometric phases, this writes
* into array \a mu the user-specified value \a not_mu instead of
* the chemical potential if the temperature is outside the range
* for which the thermo data for the one species in the phase are
* valid. The need for this arises since many condensed phases
* have thermo data fit only for the temperature range for which
* they are stable. For example, in the NASA database, the fits
* for H2O(s) are only done up to 0 C, the fits for H2O(L) are
* only done from 0 C to 100 C, etc. Using the polynomial fits outside
* the range for which the fits were done can result in spurious
* chemical potentials, and can lead to condensed phases
* appearing when in fact they should be absent.
*
* By setting \a not_mu to a large positive value, it is possible
* to force routines which seek to minimize the Gibbs free energy
* of the mixture to zero out any phases outside the temperature
* range for which their thermo data are valid.
*
* @param not_mu Value of the chemical potential to set
* species in phases, for which the thermo data
* is not valid
*
* @param mu Vector of chemical potentials
* length = Global species, units = J kmol-1
*
* @param standard If this method is called with \a standard set to true, then
* the composition-independent standard chemical potentials are
* returned instead of the composition-dependent chemical
* potentials.
*/
void getValidChemPotentials(doublereal not_mu, doublereal* mu,
bool standard = false) const;
//! Temperature [K].
doublereal temperature() const { return m_temp; }
//! Set the mixture to a state of chemical equilibrium.
/*!
* @param XY Integer flag specifying properties to hold fixed.
* @param err Error tolerance for \f$\Delta \mu/RT \f$ for
* all reactions. Also used as the relative error tolerance
* for the outer loop.
* @param maxsteps Maximum number of steps to take in solving
* the fixed TP problem.
* @param maxiter Maximum number of "outer" iterations for
* problems holding fixed something other than (T,P).
* @param loglevel Level of diagnostic output, written to a
* file in HTML format.
*/
doublereal equilibrate(int XY, doublereal err = 1.0e-9,
int maxsteps = 1000, int maxiter = 200, int loglevel = -99);
/// Set the temperature [K].
/*!
* @param T value of the temperature (Kelvin)
*/
void setTemperature(doublereal T) {
m_temp = T;
updatePhases();
}
/// Pressure [Pa].
doublereal pressure() const {
return m_press;
}
/// Volume [m^3].
/*!
* Returns the cummulative sum of the volumes of all the
* phases in the %MultiPhase.
*/
doublereal volume() const;
//! Set the pressure [Pa].
/*!
* @param P Set the pressure in the %MultiPhase object (Pa)
*/
void setPressure(doublereal P) {
m_press = P;
updatePhases();
}
/// Enthalpy [J].
doublereal enthalpy() const;
/// Entropy [J/K].
doublereal entropy() const;
/// Gibbs function [J].
doublereal gibbs() const;
/// Heat capacity at constant pressure [J/K].
doublereal cp() const;
/// Number of phases.
index_t nPhases() const {
return m_np;
}
//! Return true is species \a kGlob is a species in a
//! multicomponent solution phase.
/*!
* @param kGlob index of the global species
*/
bool solutionSpecies(index_t kGlob) const;
//! Returns the phase index of the Kth "global" species
/*!
* @param kGlob Global species index.
*
* @return
* Returns the index of the owning phase.
*/
index_t speciesPhaseIndex(index_t kGlob) const {
return m_spphase[kGlob];
}
//! Returns the mole fraction of global species k
/*!
* @param kGlob Index of the global species.
*/
doublereal moleFraction(index_t kGlob) const{
return m_moleFractions[kGlob];
}
//! Set the Mole fractions of the nth phase
/*!
* This function sets the mole fractions of the
* nth phase. Note, the mole number of the phase
* stays constant
*
* @param n ID of the phase
* @param x Vector of input mole fractions.
*/
void setPhaseMoleFractions(index_t n, doublereal* x);
//! Set the number numbers of species in the MultiPhase
/*!
* @param xMap CompositionMap of the species with
* nonzero mole numbers
* units = kmol.
*/
void setMolesByName(compositionMap& xMap);
//! Set the Moles via a string containing their names.
/*!
* The string x is in the form of a composition map
* Species which are not listed by name in the composition
* map are set to zero.
*
* @param x string x in the form of a composition map
* where values are the moles of the species.
*/
void setMolesByName(const std::string& x);
//! Return a vector of global species mole numbers
/*!
* Returns a vector of the number of moles of each species
* in the multiphase object.
*
* @param molNum Vector of doubles of length nSpecies
* containing the global mole numbers
* (kmol).
*/
void getMoles(doublereal * molNum) const;
//! Sets all of the global species mole numbers
/*!
* Sets the number of moles of each species
* in the multiphase object.
*
* @param n Vector of doubles of length nSpecies
* containing the global mole numbers
* (kmol).
*/
void setMoles(doublereal* n);
//! Retrieves a vector of element abundances
/*!
* @param elemAbundances Vector of element abundances
* Length = number of elements in the MultiPhase object.
* Index is the global element index
* units is in kmol.
*/
void getElemAbundances(doublereal * elemAbundances) const;
//! Return true if the phase \a p has valid thermo data for
//! the current temperature.
/*!
* @param p Index of the phase.
*/
bool tempOK(index_t p) const {
return m_temp_OK[p];
}
// These methods are meant for internal use.
/// update the locally-stored composition to match the current
/// compositions of the phase objects.
void updateMoleFractions();
protected:
/// Set the states of the phase objects to the locally-stored
/// state. Note that if individual phases have T and P different
/// than that stored locally, the phase T and P will be modified.
void updatePhases() const;
//! Calculate the element abundance vector
void calcElemAbundances() const;
/**
* Vector of the number of moles in each phase.
* Length = m_np, number of phases.
*/
vector_fp m_moles;
/**
* Vector of the ThermoPhase Pointers.
*/
std::vector<phase_t*> m_phase;
//! Global Stoichiometric Coefficient array
/*!
* This is a two dimensional array m_atoms(m, k). The first
* index is the global element index. The second index, k, is the
* global species index.
* The value is the number of atoms of type m in species k.
*/
array_t m_atoms;
/**
* Locally storred vector of mole fractions of all species
* comprising the MultiPhase object.
*/
vector_fp m_moleFractions;
//! Mapping between the global species number and the phase ID
/*!
* m_spphase[kGlobal] = iPhase
* Length = number of global species
*/
vector_int m_spphase;
//! Vector of ints containing of first species index in the global list of species
//! for each phase
/*!
* kfirst = m_spstart[ip], kfirst is the index of the first species in the ip'th
* phase.
*/
vector_int m_spstart;
//! String names of the global elements
/*!
* This has a length equal to the number of global elements.
*/
std::vector<std::string> m_enames;
//! Atomic number of each element
/*!
* This is the atomic number of each global element.
*/
vector_int m_atomicNumber;
//! Vector of species names in the problem
/*!
* Vector is over all species defined in the object,
* the global species index.
*/
std::vector<std::string> m_snames;
//! Returns the global element index, given the element string name
/*!
* -> used in the construction. However, wonder if it needs to be global.
*/
std::map<std::string, int> m_enamemap;
/**
* Number of phases in the MultiPhase object
*/
index_t m_np;
//! Current value of the temperature (kelvin)
doublereal m_temp;
//! Current value of the pressure (Pa)
doublereal m_press;
/**
* Number of distinct elements in all of the phases
*/
index_t m_nel;
/**
* Number of distinct species in all of the phases
*/
index_t m_nsp;
//! True if the init() routine has been called, and the MultiPhase frozen
bool m_init;
//! Global ID of the element corresponding to the electronic charge.
/*!
* If there is none, then this is equal to -1
*/
int m_eloc;
//! Vector of bools indicating whether temperatures are ok for phases.
/*!
* If the current temperature is outside the range of valid temperatures
* for the phase thermodynamics, the phase flag is set to false.
*/
mutable std::vector<bool> m_temp_OK;
//! Minimum temperature for which thermo parameterizations are valid
/*!
* Stoichiometric phases are ignored in this determination.
* units Kelvin
*/
doublereal m_Tmin;
//! Minimum temperature for which thermo parameterizations are valid
/*!
* Stoichiometric phases are ignored in this determination.
* units Kelvin
*/
doublereal m_Tmax;
//! Vector of element abundances
/*!
* m_elemAbundances[mGlobal] = kmol of element mGlobal summed over all
* species in all phases.
*/
mutable vector_fp m_elemAbundances;
};
//! Function to output a MultiPhase description to a stream
/*!
* Writes out a description of the contents of each phase of the
* MultiPhase using the report function.
*
* @param s ostream
* @param x Reference to a MultiPhase
* @return returns a reference to the ostream
*/
inline std::ostream& operator<<(std::ostream& s, Cantera::MultiPhase& x) {
size_t ip;
for (ip = 0; ip < x.nPhases(); ip++) {
if (x.phase(ip).name() != "") {
s << "*************** " << x.phase(ip).name() << " *****************" << std::endl;
}
else {
s << "*************** Phase " << ip << " *****************" << std::endl;
}
s << "Moles: " << x.phaseMoles(ip) << std::endl;
s << report(x.phase(ip)) << std::endl;
}
return s;
}
//! Choose the optimum basis of species for the equilibrium calculations.
/*!
* This is done by
* choosing the species with the largest mole fraction
* not currently a linear combination of the previous components.
* Then, calculate the stoichiometric coefficient matrix for that
* basis.
*
* Calculates the identity of the component species in the mechanism.
* Rearranges the solution data to put the component data at the
* front of the species list.
*
* Then, calculates SC(J,I) the formation reactions for all noncomponent
* species in the mechanism.
*
* Input
* ---------
* @param mphase Pointer to the multiphase object. Contains the
* species mole fractions, which are used to pick the
* current optimal species component basis.
* @param orderVectorElements
* Order vector for the elements. The element rows
* in the formula matrix are
* rearranged according to this vector.
* @param orderVectorSpecies
* Order vector for the species. The species are
* rearranged according to this formula. The first
* nCompoments of this vector contain the calculated
* species components on exit.
* @param doFormRxn If true, the routine calculates the formation
* reaction matrix based on the calculated
* component species. If false, this step is skipped.
*
* Output
* ---------
* @param usedZeroedSpecies = If true, then a species with a zero concentration
* was used as a component. The problem may be
* converged.
* @param formRxnMatrix
*
* @return Returns the number of components.
*
* @ingroup equilfunctions
*/
int BasisOptimize( int *usedZeroedSpecies, bool doFormRxn,
MultiPhase *mphase, vector_int & orderVectorSpecies,
vector_int & orderVectorElements,
vector_fp & formRxnMatrix);
//! This subroutine handles the potential rearrangement of the constraint
//! equations represented by the Formula Matrix.
/*!
* Rearrangement is only
* necessary when the number of components is less than the number of
* elements. For this case, some constraints can never be satisfied
* exactly, because the range space represented by the Formula
* Matrix of the components can't span the extra space. These
* constraints, which are out of the range space of the component
* Formula matrix entries, are migrated to the back of the Formula
* matrix.
*
* A prototypical example is an extra element column in
* FormulaMatrix[],
* which is identically zero. For example, let's say that argon is
* has an element column in FormulaMatrix[], but no species in the
* mechanism
* actually contains argon. Then, nc < ne. Unless the entry for
* desired element abundance vector for Ar is zero, then this
* element abundance constraint can never be satisfied. The
* constraint vector is not in the range space of the formula
* matrix.
* Also, without perturbation
* of FormulaMatrix[], BasisOptimize[] would produce a zero pivot
* because the matrix
* would be singular (unless the argon element column was already the
* last column of FormulaMatrix[].
* This routine borrows heavily from BasisOptimize algorithm. It
* finds nc constraints which span the range space of the Component
* Formula matrix, and assigns them as the first nc components in the
* formular matrix. This guarrantees that BasisOptimize has a
* nonsingular matrix to invert.
* input
* @param nComponents Number of components calculated previously.
*
* @param elementAbundances Current value of the element abundances
*
* @param mphase Input pointer to a MultiPhase object
*
* @param orderVectorSpecies input vector containing the ordering
* of the global species in mphase. This is used
* to extract the component basis of the mphase object.
*
* output
* @param orderVectorElements Ouput vector containing the order
* of the elements that is necessary for
* calculation of the formula matrix.
*
* @ingroup equilfunctions
*/
int ElemRearrange(int nComponents, const vector_fp & elementAbundances,
MultiPhase *mphase,
vector_int & orderVectorSpecies,
vector_int & orderVectorElements);
#ifdef DEBUG_HKM
extern int BasisOptimize_print_lvl;
#endif
}
#endif

832
Cantera/src/equil/MultiPhaseEquil.cpp Normal file
View file

@ -0,0 +1,832 @@
#include "MultiPhaseEquil.h"
#include "MultiPhase.h"
#include "sort.h"
#include "global.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. The
/// absolute value will be used in the string.
/// @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.
/// @param start If true, the initial composition will be
/// determined by a linear Gibbs minimization, otherwise the
/// initial mixture composition will be used.
MultiPhaseEquil::MultiPhaseEquil(mix_t* mix, bool start) : 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_eloc = 1000;
m_incl_species.resize(m_nsp_mix,1);
m_incl_element.resize(m_nel_mix,1);
for (m = 0; m < m_nel_mix; m++) {
string enm = mix->elementName(m);
// element 'E' or 'e' represents an electron; this
// requires special handling, so save its index
// for later use
if (enm == "E" || enm == "e") {
m_eloc = m;
}
// if an element other than electrons is not present in
// the mixture, then exclude it and all species containing
// it from the calculation. Electrons are a special case,
// since a species can have a negative number of 'atoms'
// of electrons (positive ions).
if (m_mix->elementMoles(m) <= 0.0) {
if (m != m_eloc) {
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;
}
}
}
}
}
// Now build the list of elements to be included, starting with
// electrons, if they are present.
if (m_eloc < m_nel_mix) {
m_element.push_back(m_eloc);
m_nel++;
}
// add the included elements other than electrons
for (m = 0; m < m_nel_mix; m++) {
if (m_incl_element[m] == 1 && m != m_eloc) {
m_nel++;
m_element.push_back(m);
}
}
// include pure single-constituent phases only if their thermo
// data are valid for this temperature. This is necessary,
// since some thermo polynomial fits are done only for a
// limited temperature range. For example, using the NASA
// polynomial fits for solid ice and liquid water, if this
// were not done the calculation would predict solid ice to be
// present far above its melting point, since the thermo
// polynomial fits only extend to 273.15 K, and give
// unphysical results above this temperature, leading
// (incorrectly) to Gibbs free energies at high temperature
// lower than for liquid water.
index_t ip;
for (k = 0; k < m_nsp_mix; k++) {
ip = m_mix->speciesPhaseIndex(k);
if (!m_mix->solutionSpecies(k) &&
!m_mix->tempOK(ip)) {
m_incl_species[k] = 0;
if (m_mix->speciesMoles(k) > 0.0) {
throw CanteraError("MultiPhaseEquil",
"condensed-phase species"+ m_mix->speciesName(k)
+ " is excluded since its thermo properties are \n"
"not valid at this temperature, but it has "
"non-zero moles in the initial state.");
}
}
}
// Now build the list of all species to be included in the
// calculation.
for (k = 0; k < m_nsp_mix; k++) {
if (m_incl_species[k] ==1) {
m_nsp++;
m_species.push_back(k);
}
}
// some work arrays for internal use
m_work.resize(m_nsp);
m_work2.resize(m_nsp);
m_work3.resize(m_nsp_mix);
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);
// initialize the mole numbers to the mixture composition
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);
// if the 'start' flag is set, estimate the initial mole
// numbers by doing a linear Gibbs minimization. In this case,
// only the elemental composition of the initial mixture state
// matters.
if (start) {
setInitialMoles();
}
computeN();
// Take a very small step in composition space, so that no
// species has precisely zero moles.
vector_fp dxi(m_nsp - m_nel, 1.0e-20);
multiply(m_N, DATA_PTR(dxi), DATA_PTR(m_work));
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);
}
m_force = false;
updateMixMoles();
// At this point, the instance has been created, the species
// to be included have been determined, and an initial
// composition has been selected that has all non-zero mole
// numbers for the included species.
}
doublereal MultiPhaseEquil::equilibrate(int XY, doublereal err,
int maxsteps, int loglevel) {
int i;
m_iter = 0;
string iterstr;
beginLogGroup("MultiPhaseEquil::equilibrate", loglevel);
for (i = 0; i < maxsteps; i++) {
iterstr = "iteration "+int2str(i);
beginLogGroup(iterstr);
stepComposition();
addLogEntry("error",fp2str(error()));
endLogGroup(iterstr);
if (error() < err) break;
}
if (i >= maxsteps) {
addLogEntry("Error","no convergence in "+int2str(maxsteps)
+" iterations");
endLogGroup("MultiPhaseEquil::equilibrate");
throw CanteraError("MultiPhaseEquil::equilibrate",
"no convergence in " + int2str(maxsteps) +
" iterations. Error = " + fp2str(error()));
}
addLogEntry("iterations",int2str(iterations()));
addLogEntry("error tolerance",fp2str(err));
addLogEntry("error",fp2str(error()));
endLogGroup("MultiPhaseEquil::equilibrate");
finish();
return error();
}
void MultiPhaseEquil::updateMixMoles() {
fill(m_work3.begin(), m_work3.end(), 0.0);
index_t k;
for (k = 0; k < m_nsp; k++) {
m_work3[m_species[k]] = m_moles[k];
}
m_mix->setMoles(DATA_PTR(m_work3));
}
/// Clean up the composition. The solution algorithm can leave
/// some species in stoichiometric condensed phases with very
/// small negative mole numbers. This method simply sets these to
/// zero.
void MultiPhaseEquil::finish() {
fill(m_work3.begin(), m_work3.end(), 0.0);
index_t k;
for (k = 0; k < m_nsp; k++) {
m_work3[m_species[k]] = (m_moles[k] > 0.0 ? m_moles[k] : 0.0);
}
m_mix->setMoles(DATA_PTR(m_work3));
}
/// Extimate the initial mole numbers. This is done by running
/// each reaction as far forward or backward as possible, subject
/// to the constraint that all mole numbers remain
/// non-negative. Reactions for which \f$ \Delta \mu^0 \f$ are
/// positive are run in reverse, and ones for which it is negative
/// are run in the forward direction. The end result is equivalent
/// to solving the linear programming problem of minimizing the
/// linear Gibbs function subject to the element and
/// non-negativity constraints.
int MultiPhaseEquil::setInitialMoles() {
index_t ik, j;
double not_mu = 1.0e12;
beginLogGroup("MultiPhaseEquil::setInitialMoles");
m_mix->getValidChemPotentials(not_mu, DATA_PTR(m_mu), true);
doublereal dg_rt;
int idir;
double nu;
double delta_xi, dxi_min = 1.0e10;
bool redo = true;
int iter = 0;
while (redo) {
// choose a set of components based on the current
// composition
computeN();
addLogEntry("iteration",iter);
redo = false;
iter++;
if (iter > 4) break;
// loop over all reactions
for (j = 0; j < m_nsp - m_nel; j++) {
dg_rt = 0.0;
dxi_min = 1.0e10;
for (ik = 0; ik < m_nsp; ik++) {
dg_rt += mu(ik) * m_N(ik,j);
}
// fwd or rev direction
idir = (dg_rt < 0.0 ? 1 : -1);
for (ik = 0; ik < m_nsp; ik++) {
nu = m_N(ik, j);
// set max change in progress variable by
// non-negativity requirement
if (nu*idir < 0) {
delta_xi = fabs(moles(ik)/nu);
// if a component has nearly zero moles, redo
// with a new set of components
if (!redo && delta_xi < 1.0e-10 && ik < m_nel) {
addLogEntry("component too small",speciesName(ik));
redo = true;
}
if (delta_xi < dxi_min) dxi_min = delta_xi;
}
}
// step the composition by dxi_min
for (ik = 0; ik < m_nsp; ik++) {
moles(ik) += m_N(ik, j) * idir*dxi_min;
}
}
// set the moles of the phase objects to match
updateMixMoles();
}
for (ik = 0; ik < m_nsp; ik++)
if (moles(ik) != 0.0) addLogEntry(speciesName(ik), moles(ik));
endLogGroup("MultiPhaseEquil::setInitialMoles");
return 0;
}
/// This method finds a set of component species and a complete
/// set of formation reactions for the non-components in terms of
/// the components. Note that in most cases, many different
/// component sets are possible, and therefore neither the
/// components 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 component species are taken to be the first M species
/// in array 'species' that have linearly-independent compositions.
///
/// @param order On entry, vector \a order should contain species
/// index numbers in the order of decreasing desirability as a
/// component. For example, if it is desired to choose the
/// components 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 components in declaration order, beginning with the
/// first phase added.
///
void MultiPhaseEquil::getComponents(const vector_int& order) {
index_t m, k, j;
int n;
// if the input species array has the wrong size, ignore it
// and consider the species for components in declaration 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];
}
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. This occurs when
// a species has an elemental composition that is not
// linearly independent of the component species that have
// already been assigned
if (m_A(m,m) == 0.0) {
// First, we need to find a good candidate for a
// component species to swap in for the one that has
// zero pivot. It must contain element m, be linearly
// independent of the components processed so far
// (m_A(m,k) != 0), and should be a major species if
// possible. We'll choose the species with greatest
// mole fraction that satisfies these criteria.
doublereal maxmoles = -999.0;
index_t kmax = 0;
for (k = m+1; k < nColumns; k++) {
if (m_A(m,k) != 0.0) {
if (fabs(m_moles[m_order[k]]) > maxmoles) {
kmax = k;
maxmoles = fabs(m_moles[m_order[k]]);
}
}
}
// Now exchange the column with zero pivot with the
// column for this major species
for (n = 0; n < int(nRows); n++) {
tmp = m_A(n,m);
m_A(n, m) = m_A(n, kmax);
m_A(n, kmax) = tmp;
}
// exchange the species labels on the columns
itmp = m_order[m];
m_order[m] = m_order[kmax];
m_order[kmax] = itmp;
}
// 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;
}
// For all rows below the diagonal, subtract A(n,m)/A(m,m)
// * (row m) from row n, so that A(n,m) = 0.
for (n = int(m+1); n < int(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 the m_nel columns 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 < int(m_nsp); n++) {
if (n < int(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;
}
}
}
/// 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];
}
}
void MultiPhaseEquil::printInfo() {
index_t m, ik, k;
beginLogGroup("info");
beginLogGroup("components");
for (m = 0; m < m_nel; m++) {
ik = m_order[m];
k = m_species[ik];
addLogEntry(m_mix->speciesName(k), fp2str(m_moles[ik]));
}
endLogGroup("components");
beginLogGroup("non-components");
for (m = m_nel; m < m_nsp; m++) {
ik = m_order[m];
k = m_species[ik];
addLogEntry(m_mix->speciesName(k), fp2str(m_moles[ik]));
}
endLogGroup("non-components");
addLogEntry("Error",fp2str(error()));
beginLogGroup("Delta G / RT");
for (k = 0; k < m_nsp - m_nel; k++) {
addLogEntry(reactionString(k), fp2str(m_deltaG_RT[k]));
}
endLogGroup("Delta G / RT");
endLogGroup("info");
}
/// 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;
}
void MultiPhaseEquil::step(doublereal omega, vector_fp& deltaN) {
index_t k, ik;
beginLogGroup("MultiPhaseEquil::step");
if (omega < 0.0)
throw CanteraError("step","negative omega");
for (ik = 0; ik < m_nel; ik++) {
k = m_order[ik];
m_lastmoles[k] = m_moles[k];
addLogEntry("component "+m_mix->speciesName(m_species[k])+" moles",
m_moles[k]);
addLogEntry("component "+m_mix->speciesName(m_species[k])+" step",
omega*deltaN[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])*fminn(10.0,
exp(-m_deltaG_RT[ik - m_nel]));
}
}
updateMixMoles();
endLogGroup("MultiPhaseEquil::step");
}
/// Take one step in composition, given the gradient of G at the
/// starting point, and a vector of reaction steps dxi.
doublereal MultiPhaseEquil::
stepComposition() {
beginLogGroup("MultiPhaseEquil::stepComposition");
m_iter++;
index_t ik, k = 0;
doublereal grad0 = computeReactionSteps(m_dxi);
// compute the mole fraction changes.
multiply(m_N, DATA_PTR(m_dxi), DATA_PTR(m_work));
// change to sequential form
unsort(m_work);
// scale omega to keep the major species non-negative
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 (ik < m_nel) {
FCTR = 0.99;
if (m_moles[k] < MAJOR_THRESHOLD) m_force = true;
}
else FCTR = 0.9;
// 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)) {
if (m_moles[k] < MAJOR_THRESHOLD) m_force = true;
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) {
m_force = true;
}
}
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) {
m_force = true;
}
}
}
if (m_moles[k] < -Tiny) {
addLogEntry("Negative moles for "
+m_mix->speciesName(m_species[k]), fp2str(m_moles[k]));
}
m_majorsp[k] = true;
}
}
// now take a step with this scaled omega
addLogEntry("Stepping by ", fp2str(omegamax));
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.
doublereal not_mu = 1.0e12;
m_mix->getValidChemPotentials(not_mu, DATA_PTR(m_mu));
doublereal grad1 = 0.0;
for (k = 0; k < m_nsp; k++) {
grad1 += m_work[k] * m_mu[m_species[k]];
}
omega = omegamax;
if (grad1 > 0.0) {
omega *= fabs(grad0) / (grad1 + fabs(grad0));
for (k = 0; k < m_nsp; k++) m_moles[k] = m_lastmoles[k];
addLogEntry("Stepped over minimum. Take smaller step ", fp2str(omega));
step(omega, m_work);
}
printInfo();
endLogGroup("MultiPhaseEquil::stepComposition");
return omega;
}
/// Compute the change in extent of reaction for each reaction.
doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi) {
index_t j, k, ik, kc, ip;
doublereal stoich, nmoles, csum, term1, fctr, rfctr;
vector_fp nu;
const doublereal TINY = 1.0e-20;
doublereal grad = 0.0;
dxi.resize(m_nsp - m_nel);
computeN();
doublereal not_mu = 1.0e12;
m_mix->getValidChemPotentials(not_mu, DATA_PTR(m_mu));
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.5;
}
}
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) {
// bug fixed 7/12/06 DGG
stoich = nu[k]; // nu[kc];
psum += stoich * stoich;
}
}
sum -= psum / (fabs(m_mix->phaseMoles(ip)) + TINY);
// if (ISNAN(sum)) {
// cout << " sum is nan. " << endl;
// cout << psum << " " << m_mix->phaseMoles(ip) << endl;
// }
}
}
rfctr = term1 + csum + sum;
if (fabs(rfctr) < TINY)
fctr = 1.0;
else
fctr = 1.0/(term1 + csum + sum);
// if (ISNAN(fctr)) {
// cout << "fctr is nan." << endl;
// cout << term1 << " " << csum << " " << sum << " " << TINY << endl;
//}
}
dxi[j] = -fctr*dg_rt;
//if (ISNAN(dxi[j])) {
// cout << "nan detected. " << endl;
// cout << fctr << " " << dg_rt << endl;
//}
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;
}
grad += dxi[j]*dg_rt;
}
return grad*GasConstant*m_temp;
}
void MultiPhaseEquil::computeN() {
index_t m, k;
// 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 (int(k) == m_order[ij]) ok = true;
}
if (!ok || m_force) {
getComponents(m_sortindex);
m_force = true;
break;
}
}
}
doublereal MultiPhaseEquil::error() {
index_t j, ik, k;
doublereal err, maxerr = 0.0;
// examine every reaction
for (j = 0; j < m_nsp - m_nel; j++) {
ik = j + m_nel;
k = m_order[ik];
// don't require formation reactions for solution species
// present in trace amounts to be equilibrated
if (!isStoichPhase(ik) && fabs(moles(ik)) <= SmallNumber) {
err = 0.0;
}
// for stoichiometric phase species, no error if not present and
// delta G for the formation reaction is positive
if (isStoichPhase(ik) && moles(ik) <= 0.0 &&
m_deltaG_RT[j] >= 0.0) {
err = 0.0;
}
else {
err = fabs(m_deltaG_RT[j]);
}
if (err > maxerr) {
maxerr = err;
}
}
return maxerr;
}
}

124
Cantera/src/equil/MultiPhaseEquil.h Normal file
View file

@ -0,0 +1,124 @@
#ifndef CT_MULTIPHASE_EQUIL
#define CT_MULTIPHASE_EQUIL
#include "ct_defs.h"
#include "MultiPhase.h"
namespace Cantera {
/**
* Multiphase chemical equilibrium solver. Class MultiPhaseEquil
* is designed to be used to set a mixture containing one or more
* phases to a state of chemical equilibrium. It implements the
* VCS algorithm, described in Smith and Missen, "Chemical
* Reaction Equilibrium."
*
* This class only handles chemical equilibrium at a specified
* temperature and pressure. To compute equilibrium holding other
* properties fixed, it is necessary to iterate on T and P in an
* "outer" loop, until the specified properties have the desired
* values. This is done, for example, in method equilibrate of
* class MultiPhase.
*
* This class is primarily meant to be used internally by the
* equilibrate method of class MultiPhase, although there is no
* reason it cannot be used directly in application programs if
* desired.
*
* @ingroup equil
*/
class MultiPhaseEquil {
public:
typedef MultiPhase mix_t;
typedef size_t index_t;
typedef DenseMatrix matrix_t;
MultiPhaseEquil(mix_t* mix, bool start=true);
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 loglevel=-99);
std::string reactionString(index_t j);
doublereal error();
void printInfo();
void setInitialMixMoles() {
setInitialMoles();
finish();
}
index_t componentIndex(index_t n) { return m_species[m_order[n]]; }
protected:
void getComponents(const vector_int& order);
int setInitialMoles();
void computeN();
doublereal stepComposition();
//void sort(vector_fp& x);
void unsort(vector_fp& x);
void step(doublereal omega, vector_fp& deltaN);
doublereal computeReactionSteps(vector_fp& dxi);
void updateMixMoles();
void finish();
// moles of the species with sorted index ns
double moles(int ns) const { return m_moles[m_order[ns]]; }
double& moles(int ns) { return m_moles[m_order[ns]]; }
int solutionSpecies(int n) const { return m_dsoln[m_order[n]]; }
bool isStoichPhase(int n) const { return (m_dsoln[m_order[n]] == 0); }
doublereal mu(int n) const { return m_mu[m_species[m_order[n]]]; }
std::string speciesName(int n) const { return
m_mix->speciesName(m_species[m_order[n]]); }
index_t m_nel_mix, m_nsp_mix, m_np;
index_t m_nel, m_nsp;
index_t m_eloc;
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, m_work3;
vector_fp m_moles, m_lastmoles, m_dxi;
vector_fp m_deltaG_RT, m_mu;
std::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 of indices for species that are included in the
// calculation. This is used to exclude pure-phase species
// with invalid thermo data
vector_int m_species;
vector_int m_element;
std::vector<bool> m_solnrxn;
bool m_force;
};
}
#endif

86
Cantera/src/equil/PropertyCalculator.h Executable file
View file

@ -0,0 +1,86 @@
/**
* @file PropertyCalculator.h
*
* $Author$
* $Revision$
* $Date$
*/
// Copyright 2001 California Institute of Technology
#ifndef CT_PROP_CALC_H
#define CT_PROP_CALC_H
#include "ct_defs.h"
namespace Cantera {
/// Classes used by ChemEquil. These classes are used only by the
/// ChemEquil equilibrium solver. Each one returns a particular
/// property of the object supplied as the argument.
///
template<class M>
class PropertyCalculator {
public:
virtual ~PropertyCalculator(){}
virtual doublereal value(const M& s) =0;
virtual string symbol() =0;
};
template<class M>
class EnthalpyCalculator : public PropertyCalculator<M> {
public:
virtual doublereal value(const M& s) {
return s.enthalpy_mass();
}
virtual string symbol() { return "H"; }
};
template<class M>
class EntropyCalculator : public PropertyCalculator<M> {
public:
virtual doublereal value(const M& s) {
return s.entropy_mass();
}
virtual string symbol() { return "S"; }
};
template<class M>
class TemperatureCalculator : public PropertyCalculator<M> {
public:
virtual doublereal value(const M& s) {
return s.temperature();
}
virtual string symbol() { return "T"; }
};
template<class M>
class PressureCalculator : public PropertyCalculator<M> {
public:
virtual doublereal value(const M& s) {
return s.pressure();
}
virtual string symbol() { return "P"; }
};
template<class M>
class DensityCalculator : public PropertyCalculator<M> {
public:
virtual doublereal value(const M& s) {
return s.density();
}
virtual string symbol() { return "V"; }
};
template<class M>
class IntEnergyCalculator : public PropertyCalculator<M> {
public:
virtual doublereal value(const M& s) {
return s.intEnergy_mass();
}
virtual string symbol() { return "U"; }
};
}
#endif

105
Cantera/src/equil/equil.h Normal file
View file

@ -0,0 +1,105 @@
/***********************************************************************
* $RCSfile$
* $Author$
* $Date$
* $Revision$
***********************************************************************/
// Copyright 2001 California Institute of Technology
/**
* @file equil.h
* This file contains the definition of some high level general equilibration
* routines and the text for the module \ref equilfunctions.
*
* It also contains the Module doxygen text for the Equilibration Solver
* capability within %Cantera. see \ref equilfunctions
*/
#ifndef CT_KERNEL_EQUIL_H
#define CT_KERNEL_EQUIL_H
//#include "ChemEquil.h"
#include "MultiPhase.h"
namespace Cantera {
/*!
* @defgroup equilfunctions Equilibrium Solver Capability
*
* Cantera has several different equilibrium routines.
*/
//-----------------------------------------------------------
// convenience functions
//-----------------------------------------------------------
//! Equilibrate a ThermoPhase object
/*!
* Set a single-phase chemical solution to chemical equilibrium.
* This is a convenience function that uses one or the other of
* the two chemical equilibrium solvers. The XY parameter indicates what two
* thermodynamic quantities, other than element composition, are to be held
* constant during the equilibration process.
*
* @param s ThermoPhase object that will be equilibrated.
* @param XY String representation of what two properties
* are being held constant
* @param solver ID of the solver to be used to equlibrate the phase.
* If solver = 0, the ChemEquil solver will be used,
* and if solver = 1, the
* MultiPhaseEquil solver will be used (slower than ChemEquil,
* but more stable). If solver < 0 (default, then ChemEquil will
* be tried first, and if it fails MultiPhaseEquil will be tried.
* @param rtol Relative tolerance
* @param maxsteps Maximum number of steps to take to find the solution
* @param maxiter For the MultiPhaseEquil solver only, this is
* the maximum number of outer temperature or pressure iterations
* to take when T and/or P is not held fixed.
* @param loglevel loglevel Controls amount of diagnostic output. loglevel
* = 0 suppresses diagnostics, and increasingly-verbose messages
* are written as loglevel increases. The messages are written to
* a file in HTML format for viewing in a web browser.
* @see HTML_logs
*
* @return
* Return variable is equal to the number of subroutine attempts
* it took to equilibrate the system.
*
*
* @ingroup equilfunctions
* @ingroup equil
*/
int equilibrate(thermo_t& s, const char* XY,
int solver = -1, doublereal rtol = 1.0e-9, int maxsteps = 1000,
int maxiter = 100, int loglevel = -99);
//! Equilibrate a MultiPhase object
/*!
* Equilibrate a MultiPhase object. The XY parameter indicates what two
* thermodynamic quantities, other than element composition, are to be held
* constant during the equilibration process.
*
* This is the top-level driver for multiphase equilibrium. It
* doesn't do much more than call the equilibrate method of class
* MultiPhase, except that it adds some messages to the logfile,
* if loglevel is set > 0.
*
* @param s MultiPhase object that will be equilibrated.
* @param XY String representation of what is being held constant
* @param rtol Relative tolerance
* @param maxsteps Maximum number of steps
* @param maxiter Maximum iterations
* @param loglevel loglevel
*
* @return
* Return variable is equal to the number of subroutine attempts
* it took to equilibrate the system.
*
* @ingroup equilfunctions
* @ingroup equil
*/
doublereal equilibrate(MultiPhase& s, const char* XY,
doublereal rtol = 1.0e-9, int maxsteps = 1000, int maxiter = 100,
int loglevel = -99);
}
#endif

187
Cantera/src/equil/equilibrate.cpp Normal file
View file

@ -0,0 +1,187 @@
/**
* @file equilibrate.cpp
* Driver routines for the chemical equilibrium solvers.
*
*/
#include "equil.h"
#include "ChemEquil.h"
#include "MultiPhaseEquil.h"
namespace Cantera {
/*
* Set a multiphase mixture to a state of chemical equilibrium.
* This is the top-level driver for multiphase equilibrium. It
* doesn't do much more than call the equilibrate method of class
* MultiPhase, except that it adds some messages to the logfile,
* if loglevel is set > 0.
*
* @ingroup equil
*/
doublereal equilibrate(MultiPhase& s, const char* XY,
doublereal tol, int maxsteps, int maxiter,
int loglevel) {
beginLogGroup("equilibrate",loglevel);
addLogEntry("multiphase equilibrate function");
beginLogGroup("arguments");
addLogEntry("XY",XY);
addLogEntry("tol",tol);
addLogEntry("maxsteps",maxsteps);
addLogEntry("maxiter",maxiter);
addLogEntry("loglevel",loglevel);
endLogGroup("arguments");
s.init();
int ixy = _equilflag(XY);
if (ixy == TP || ixy == HP || ixy == SP || ixy == TV) {
try {
double err = s.equilibrate(ixy, tol, maxsteps, maxiter);
addLogEntry("Success. Error",err);
endLogGroup("equilibrate");
return err;
}
catch (CanteraError e) {
addLogEntry("Failure.",lastErrorMessage());
endLogGroup("equilibrate");
throw e;
}
}
else {
addLogEntry("multiphase equilibrium can be done only for TP, HP, SP, or TV");
endLogGroup("equilibrate");
throw CanteraError("equilibrate","unsupported option");
return -1.0;
}
}
/*
* Set a single-phase chemical solution to chemical equilibrium.
* This is a convenience function that uses one or the other of
* the two chemical equilibrium solvers.
*
* @param s The object to set to an equilibrium state
*
* @param XY An integer specifying the two properties to be held
* constant.
*
* @param solver The equilibrium solver to use. If solver = 0,
* the ChemEquil solver will be used, and if solver = 1, the
* MultiPhaseEquil solver will be used (slower than ChemEquil,
* but more stable). If solver < 0 (default, then ChemEquil will
* be tried first, and if it fails MultiPhaseEquil will be tried.
*
* @param maxsteps The maximum number of steps to take to find
* the solution.
*
* @param maxiter For the MultiPhaseEquil solver only, this is
* the maximum number of outer temperature or pressure iterations
* to take when T and/or P is not held fixed.
*
* @param loglevel Controls amount of diagnostic output. loglevel
* = 0 suppresses diagnostics, and increasingly-verbose messages
* are written as loglevel increases. The messages are written to
* a file in HTML format for viewing in a web browser.
* @see HTML_logs
*
* @ingroup equil
*/
int equilibrate(thermo_t& s, const char* XY, int solver,
doublereal rtol, int maxsteps, int maxiter, int loglevel) {
MultiPhase* m = 0;
ChemEquil* e = 0;
bool redo = true;
int retn = -1;
int nAttempts = 0;
int retnSub = 0;
beginLogGroup("equilibrate", loglevel);
addLogEntry("Single-phase equilibrate function");
{
beginLogGroup("arguments");
addLogEntry("phase",s.id());
addLogEntry("XY",XY);
addLogEntry("solver",solver);
addLogEntry("rtol",rtol);
addLogEntry("maxsteps",maxsteps);
addLogEntry("maxiter",maxiter);
addLogEntry("loglevel",loglevel);
endLogGroup("arguments");
}
while (redo) {
if (solver > 0) {
m = new MultiPhase;
try {
m->addPhase(&s, 1.0);
m->init();
nAttempts++;
(void) equilibrate(*m, XY, rtol, maxsteps, maxiter, loglevel);
redo = false;
addLogEntry("MultiPhaseEquil solver succeeded.");
delete m;
retn = nAttempts;
}
catch (CanteraError err) {
addLogEntry("MultiPhaseEquil solver failed.");
delete m;
if (nAttempts < 2) {
addLogEntry("Trying single phase ChemEquil solver.");
solver = -1;
} else {
endLogGroup("equilibrate");
throw err;
}
}
}
else { // solver <= 0
/*
* Call the element potential solver
*/
e = new ChemEquil;
try {
e->options.maxIterations = maxsteps;
e->options.relTolerance = rtol;
nAttempts++;
retnSub = e->equilibrate(s,XY);
if (retnSub < 0) {
addLogEntry("ChemEquil solver failed.");
if (nAttempts < 2) {
addLogEntry("Trying MultiPhaseEquil solver.");
solver = 1;
} else {
throw CanteraError("equilibrate",
"Both equilibrium solvers failed");
}
}
retn = nAttempts;
s.setElementPotentials(e->elementPotentials());
redo = false;
delete e;
addLogEntry("ChemEquil solver succeeded.");
}
catch (CanteraError err) {
delete e;
addLogEntry("ChemEquil solver failed.");
// If ChemEquil fails, try the MultiPhase solver
if (solver < 0) {
addLogEntry("Trying MultiPhaseEquil solver.");
solver = 1;
}
else {
redo = false;
endLogGroup("equilibrate");
throw err;
}
}
}
} // while (redo)
/*
* We are here only for a success
*/
endLogGroup("equilibrate");
return retn;
}
}