diff --git a/Cantera/src/equil/BasisOptimize.cpp b/Cantera/src/equil/BasisOptimize.cpp new file mode 100644 index 000000000..02cc48cff --- /dev/null +++ b/Cantera/src/equil/BasisOptimize.cpp @@ -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() ****************************************************/ diff --git a/Cantera/src/equil/ChemEquil.cpp b/Cantera/src/equil/ChemEquil.cpp new file mode 100755 index 000000000..86ee84243 --- /dev/null +++ b/Cantera/src/equil/ChemEquil.cpp @@ -0,0 +1,1906 @@ +/** + * + * @file ChemEquil.cpp + * + * Chemical equilibrium. Implementation file for class + * ChemEquil. + * + * + * $Id$ + * + * Copyright 2001 California Institute of Technology + * + */ + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include +using namespace std; + +#include "ChemEquil.h" +#include "DenseMatrix.h" + +#include "sort.h" +#include "PropertyCalculator.h" +#include "ctexceptions.h" +#include "vec_functions.h" +#include "stringUtils.h" +#include "MultiPhase.h" + +#ifdef DEBUG_HKM +#include "stdio.h" +int Cantera::ChemEquil_print_lvl = 0; +static char sbuf[1024]; +#endif +#ifndef MIN +#define MIN(x,y) (( (x) < (y) ) ? (x) : (y)) +#endif +namespace Cantera { + + /// map property strings to integers + int _equilflag(const char* xy) { + string flag = string(xy); + if (flag == "TP") return TP; + else if (flag == "TV") return TV; + else if (flag == "HP") return HP; + else if (flag == "UV") return UV; + else if (flag == "SP") return SP; + else if (flag == "SV") return SV; + else if (flag == "UP") return UP; + else throw CanteraError("_equilflag","unknown property pair "+flag); + } + + + //----------------------------------------------------------- + // construction / destruction + //----------------------------------------------------------- + + + /// Default Constructor. + ChemEquil::ChemEquil() : m_skip(-1), m_p1(0), m_p2(0), m_elementTotalSum(1.0), + m_p0(OneAtm), m_eloc(-1), + m_elemFracCutoff(1.0E-100), + m_doResPerturb(false) + {} + + //! 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::ChemEquil(thermo_t& s) : + m_skip(-1), m_p1(0), m_p2(0), + m_elementTotalSum(1.0), + m_p0(OneAtm), m_eloc(-1), + m_elemFracCutoff(1.0E-100), + m_doResPerturb(false) + { + initialize(s); + } + + /// Destructor + ChemEquil::~ChemEquil(){ + if (m_p1) { + delete m_p1; + } + if (m_p2) { + delete m_p2; + } + } + + /** + * Prepare for equilibrium calculations. + * @param s object representing the solution phase. + */ + void ChemEquil::initialize(thermo_t& s) + { + // store a pointer to s and some of its properties locally. + m_phase = &s; + + m_p0 = s.refPressure(); + m_kk = s.nSpecies(); + m_mm = s.nElements(); + m_nComponents = m_mm; + //if (m_kk < m_mm) { + //throw CanteraError("ChemEquil::initialize", + // "number of species cannot be less than the number of elements."); + //} + + // allocate space in internal work arrays within the ChemEquil object + m_molefractions.resize(m_kk); + m_lambda.resize(m_mm, -100.0); + m_elementmolefracs.resize(m_mm); + m_comp.resize(m_mm * m_kk); + m_jwork1.resize(m_mm+2); + m_jwork2.resize(m_mm+2); + m_startSoln.resize(m_mm+1); + m_grt.resize(m_kk); + m_mu_RT.resize(m_kk); + m_muSS_RT.resize(m_kk); + m_component.resize(m_mm,-2); + m_orderVectorElements.resize(m_mm); + int m, k; + for (m = 0; m < m_mm; m++) { + m_orderVectorElements[m] = m; + } + m_orderVectorSpecies.resize(m_kk); + for (k = 0; k < m_kk; k++) { + m_orderVectorSpecies[k] = k; + } + + // set up elemental composition matrix + int mneg = -1; + doublereal na, ewt; + for (m = 0; m < m_mm; m++) { + for (k = 0; k < m_kk; k++) { + na = s.nAtoms(k,m); + + // handle the case of negative atom numbers (used to + // represent positive ions, where the 'element' is an + // electron + if (na < 0.0) { + + // if negative atom numbers have already been specified + // for some element other than this one, throw + // an exception + if (mneg >= 0 && mneg != m) + throw CanteraError("ChemEquil::initialize", + "negative atom numbers allowed for only one element"); + mneg = m; + ewt = s.atomicWeight(m); + + // the element should be an electron... if it isn't + // print a warning. + if (ewt > 1.0e-3) + writelog(string("WARNING: species " + +s.speciesName(k) + +" has "+fp2str(s.nAtoms(k,m)) + +" atoms of element " + +s.elementName(m)+ + ", but this element is not an electron.\n")); + } + } + } + m_eloc = mneg; + + // set up the elemental composition matrix + for (k = 0; k < m_kk; k++) { + for (m = 0; m < m_mm; m++) { + m_comp[k*m_mm + m] = s.nAtoms(k,m); + } + } + } + + + /** + * Set mixture to an equilibrium state consistent with specified + * element potentials and temperature. + * + * @param lambda_RT vector of non-dimensional element potentials + * \f[ \lambda_m/RT \f]. + * @param t temperature in K. + * + */ + void ChemEquil::setToEquilState(thermo_t& s, + const vector_fp& lambda_RT, doublereal t) + { + // Construct the chemical potentials by summing element potentials + fill(m_mu_RT.begin(), m_mu_RT.end(), 0.0); + for (int k = 0; k < m_kk; k++) + for (int m = 0; m < m_mm; m++) + m_mu_RT[k] += lambda_RT[m]*nAtoms(k,m); + + // Set the temperature + s.setTemperature(t); + + // Call the phase-specific method to set the phase to the + // equilibrium state with the specified species chemical + // potentials. + s.setToEquilState(DATA_PTR(m_mu_RT)); + update(s); + } + + + /** + * update internally stored state information. + * @todo argument not used. + */ + void ChemEquil::update(const thermo_t& s) { + + // get the mole fractions, temperature, and density + s.getMoleFractions(DATA_PTR(m_molefractions)); + m_temp = s.temperature(); + m_dens = s.density(); + + // compute the elemental mole fractions + double sum = 0.0; + int m, k; + for (m = 0; m < m_mm; m++) { + m_elementmolefracs[m] = 0.0; + for (k = 0; k < m_kk; k++) { + m_elementmolefracs[m] += nAtoms(k,m) * m_molefractions[k]; + if (m_molefractions[k] < 0.0) { + throw CanteraError("update", + "negative mole fraction for "+s.speciesName(k)+ + ": "+fp2str(m_molefractions[k])); + } + } + sum += m_elementmolefracs[m]; + } + // Store the sum for later use + m_elementTotalSum = sum; + // normalize the element mole fractions + for (m = 0; m < m_mm; m++) m_elementmolefracs[m] /= sum; + } + + /// Estimate the initial mole numbers. This version borrows from the + /// MultiPhaseEquil solver. + int ChemEquil::setInitialMoles(thermo_t& s, vector_fp & elMoleGoal) { + MultiPhase* mp = 0; + MultiPhaseEquil* e = 0; + int iok = 0; + beginLogGroup("ChemEquil::setInitialMoles"); + try { + mp = new MultiPhase; + mp->addPhase(&s, 1.0); + mp->init(); + e = new MultiPhaseEquil(mp, true); + e->setInitialMixMoles(); + + // store component indices + if (m_nComponents > m_kk) { + m_nComponents = m_kk; + } + for (int m = 0; m < m_nComponents; m++) { + m_component[m] = e->componentIndex(m); + } + for (int k = 0; k < m_kk; k++) { + if (s.moleFraction(k) > 0.0) { + addLogEntry(s.speciesName(k), + s.moleFraction(k)); + } + } + /* + * Update the current values of the temp, density, and + * mole fraction, and element abundance vectors kept + * within the ChemEquil object. + */ + update(s); + +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + writelog("setInitialMoles: Estimated Mole Fractions\n"); + sprintf(sbuf," Temperature = %g\n", s.temperature()); writelog(sbuf); + sprintf(sbuf," Pressure = %g\n", s.pressure()); writelog(sbuf); + for (int k = 0; k < m_kk; k++) { + string nnn = s.speciesName(k); + double mf = s.moleFraction(k); + sprintf(sbuf," %-12s % -10.5g\n", nnn.c_str(), mf); writelog(sbuf); + } + writelog(" Element_Name ElementGoal ElementMF\n"); + for (int m = 0; m < m_mm; m++) { + string nnn = s.elementName(m); + sprintf(sbuf," %-12s % -10.5g% -10.5g\n", + nnn.c_str(), elMoleGoal[m], m_elementmolefracs[m]); writelog(sbuf); + } + } +#endif + + delete e; + delete mp; + iok = 0; + } + catch (CanteraError) { + delete e; + delete mp; + iok = -1; + } + endLogGroup(); + return iok; + } + + + /** + * Generate a starting estimate for the element potentials. + */ + int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, + vector_fp& elMolesGoal) + { + int m, n; + beginLogGroup("estimateElementPotentials"); + //for (k = 0; k < m_kk; k++) { + // if (m_molefractions[k] > 0.0) { + // m_molefractions[k] = fmaxx(m_molefractions[k], 0.05); + // } + //} + //s.setState_PX(s.pressure(), m_molefractions.begin()); + + + vector_fp b(m_mm, -999.0); + vector_fp mu_RT(m_kk, 0.0); + vector_fp xMF_est(m_kk, 0.0); + + s.getMoleFractions(DATA_PTR(xMF_est)); + for (n = 0; n < s.nSpecies(); n++) { + if (xMF_est[n] < 1.0E-20) { + xMF_est[n] = 1.0E-20; + } + } + s.setMoleFractions(DATA_PTR(xMF_est)); + s.getMoleFractions(DATA_PTR(xMF_est)); + + MultiPhase *mp = new MultiPhase; + mp->addPhase(&s, 1.0); + mp->init(); + int usedZeroedSpecies = 0; + vector_fp formRxnMatrix; + m_nComponents = BasisOptimize(&usedZeroedSpecies, false, + mp, m_orderVectorSpecies, + m_orderVectorElements, formRxnMatrix); + + for (m = 0; m < m_nComponents; m++) { + int k = m_orderVectorSpecies[m]; + m_component[m] = k; + if (xMF_est[k] < 1.0E-8) { + xMF_est[k] = 1.0E-8; + } + } + s.setMoleFractions(DATA_PTR(xMF_est)); + s.getMoleFractions(DATA_PTR(xMF_est)); + + int nct = Cantera::ElemRearrange(m_nComponents, elMolesGoal, mp, + m_orderVectorSpecies, m_orderVectorElements); + if (nct != m_nComponents) { + throw CanteraError("ChemEquil::estimateElementPotentials", + "confused"); + } + + delete mp; + + + s.getChemPotentials(DATA_PTR(mu_RT)); + doublereal rrt = 1.0/(GasConstant* s.temperature()); + scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(), rrt); + +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + for (m = 0; m < m_nComponents; m++) { + int isp = m_component[m]; + string nnn = s.speciesName(isp); + sprintf(sbuf,"isp = %d, %s\n", isp, nnn.c_str()); + writelog(sbuf); + } + double pres = s.pressure(); + double temp = s.temperature(); + sprintf(sbuf,"Pressure = %g\n", pres); writelog(sbuf); + sprintf(sbuf,"Temperature = %g\n", temp); writelog(sbuf); + writelog(" id Name MF mu/RT \n"); + for (n = 0; n < s.nSpecies(); n++) { + string nnn = s.speciesName(n); + sprintf(sbuf,"%10d %15s %10.5g %10.5g\n", + n, nnn.c_str(), xMF_est[n], mu_RT[n]); + writelog(sbuf); + } + } +#endif + DenseMatrix aa(m_nComponents, m_nComponents, 0.0); + for (m = 0; m < m_nComponents; m++) { + for (n = 0; n < m_nComponents; n++) { + aa(m,n) = nAtoms(m_component[m], m_orderVectorElements[n]); + } + b[m] = mu_RT[m_component[m]]; + } + + int info; + try { + info = solve(aa, DATA_PTR(b)); + } + catch (CanteraError) { + addLogEntry("failed to estimate initial element potentials."); + info = -2; + } + for (m = 0; m < m_nComponents; m++) { + lambda_RT[m_orderVectorElements[m]] = b[m]; + } + for (m = m_nComponents; m < m_mm; m++) { + lambda_RT[m_orderVectorElements[m]] = 0.0; + } + if (info == 0) { + for (m = 0; m < m_mm; m++) { + addLogEntry(s.elementName(m),lambda_RT[m]); + } + } + +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + writelog(" id CompSpecies ChemPot EstChemPot Diff\n"); + for (m = 0; m < m_nComponents; m++) { + int isp = m_component[m]; + double tmp = 0.0; + string sname = s.speciesName(isp); + for (n = 0; n < m_mm; n++) { + tmp += nAtoms(isp, n) * lambda_RT[n]; + } + sprintf(sbuf,"%3d %16s %10.5g %10.5g %10.5g\n", + m, sname.c_str(), mu_RT[isp], tmp, tmp - mu_RT[isp]); + writelog(sbuf); + } + + writelog(" id ElName Lambda_RT\n"); + for (m = 0; m < m_mm; m++) { + string ename = s.elementName(m); + sprintf(sbuf," %3d %6s %10.5g\n", m, ename.c_str(), lambda_RT[m]); + writelog(sbuf); + } + } +#endif + endLogGroup(); + return info; + } + + + /** + * Equilibrate a phase, holding the elemental composition fixed + * at the initial value found within the ThermoPhase object. + * + * The value of 2 specified properties are obtained by querying the + * ThermoPhase object. The properties must be already contained + * within the current thermodynamic state of the system. + */ + int ChemEquil::equilibrate(thermo_t& s, const char* XY, + bool useThermoPhaseElementPotentials) { + vector_fp elMolesGoal(s.nElements()); + initialize(s); + update(s); + copy(m_elementmolefracs.begin(), m_elementmolefracs.end(), + elMolesGoal.begin()); + return equilibrate(s, XY, elMolesGoal, useThermoPhaseElementPotentials); + } + + + /** + * Compute the equilibrium composition for 2 specified + * properties and the specified element moles. + * + * elMoles = specified vector of element abundances. + * + * The 2 specified properties are obtained by querying the + * ThermoPhase object. The properties must be already contained + * within the current thermodynamic state of the system. + * + * Return variable: + * Successful returns are indicated by a return value of 0. + * Unsuccessful returns are indicated by a return value of -1 for + * lack of convergence or -3 for a singular jacobian. + */ + int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, vector_fp& elMolesGoal, + bool useThermoPhaseElementPotentials) + { + doublereal xval, yval, tmp; + int fail = 0; + int m, im; + + if (m_p1) delete m_p1; + if (m_p2) delete m_p2; + bool tempFixed = true; + int XY = _equilflag(XYstr); + + vector_fp state; + s.saveState(state); + + /* + * Check Compatibility + */ + if (m_mm != s.nElements() || m_kk != s.nSpecies()) { + throw CanteraError("ChemEquil::equilibrate ERROR", + "Input ThermoPhase is incompatible with initialization"); + } + +#ifdef DEBUG_HKM + int n; + const vector& eNames = s.elementNames(); +#endif + beginLogGroup("ChemEquil::equilibrate"); + initialize(s); + update(s); + switch (XY) { + case TP: case PT: + m_p1 = new TemperatureCalculator; + m_p2 = new PressureCalculator; + break; + case HP: case PH: + tempFixed = false; + m_p1 = new EnthalpyCalculator; + m_p2 = new PressureCalculator; + break; + case SP: case PS: + tempFixed = false; + m_p1 = new EntropyCalculator; + m_p2 = new PressureCalculator; + break; + case SV: case VS: + tempFixed = false; + m_p1 = new EntropyCalculator; + m_p2 = new DensityCalculator; + break; + case TV: case VT: + m_p1 = new TemperatureCalculator; + m_p2 = new DensityCalculator; + break; + case UV: case VU: + tempFixed = false; + m_p1 = new IntEnergyCalculator; + m_p2 = new DensityCalculator; + break; + default: + endLogGroup("ChemEquil::equilibrate"); + throw CanteraError("equilibrate","illegal property pair."); + } + + addLogEntry("Problem type","fixed "+m_p1->symbol()+", "+m_p2->symbol()); + addLogEntry(m_p1->symbol(), m_p1->value(s)); + addLogEntry(m_p2->symbol(), m_p2->value(s)); + + // If the temperature is one of the specified variables, and + // it is outside the valid range, throw an exception. + if (tempFixed) { + double tfixed = s.temperature(); + if (tfixed > s.maxTemp() + 1.0 || tfixed < s.minTemp() - 1.0) { + endLogGroup("ChemEquil::equilibrate"); + throw CanteraError("ChemEquil","Specified temperature (" + +fp2str(s.temperature())+" K) outside " + "valid range of "+fp2str(s.minTemp())+" K to " + +fp2str(s.maxTemp())+" K\n"); + } + } + + /* + * Before we do anything to change the ThermoPhase object, + * we calculate and store the two specified thermodynamic + * properties that we are after. + */ + xval = m_p1->value(s); + yval = m_p2->value(s); + + int mm = m_mm; + int nvar = mm + 1; + DenseMatrix jac(nvar, nvar); // jacobian + vector_fp x(nvar, -102.0); // solution vector + vector_fp res_trial(nvar, 0.0); // residual + + /* + * Replace one of the element abundance fraction equations + * with the specified property calculation. + * + * We choose the equation of the element with the highest element + * abundance. + */ + tmp = -1.0; + for (im = 0; im < m_nComponents; im++) { + m = m_orderVectorElements[im]; + if (elMolesGoal[m] > tmp ) { + m_skip = m; + tmp = elMolesGoal[m]; + } + } + if (tmp <= 0.0) { + throw CanteraError("ChemEquil", + "Element Abundance Vector is zeroed"); + } + + // start with a composition with everything non-zero. Note + // that since we have already save the target element moles, + // changing the composition at this point only affects the + // starting point, not the final solution. + vector_fp xmm(m_kk,0.0); + for (int k = 0; k < m_kk; k++) { + xmm[k] = s.moleFraction(k) + 1.0E-32; + } + s.setMoleFractions(DATA_PTR(xmm)); + + /* + * Update the internally storred values of m_temp, + * m_dens, and the element mole fractions. + */ + update(s); + + // loop to estimate T + if (!tempFixed) { + + beginLogGroup("Initial T Estimate"); + + doublereal tmax = s.maxTemp(); + doublereal tmin = s.minTemp(); + doublereal slope, phigh, plow, pval, dt; + + // first get the property values at the upper and lower + // temperature limits. Since p1 (h, s, or u) is monotonic + // in T, these values determine the upper and lower + // bounnds (phigh, plow) for p1. + + s.setTemperature(tmax); + setInitialMoles(s, elMolesGoal); + phigh = m_p1->value(s); + + s.setTemperature(tmin); + setInitialMoles(s, elMolesGoal); + plow = m_p1->value(s); + + // start with T at the midpoint of the range + doublereal t0 = 0.5*(tmin + tmax); + s.setTemperature(t0); + + // loop up to 5 times + for (int it = 0; it < 5; it++) { + + // set the composition and get p1 + setInitialMoles(s, elMolesGoal); + pval = m_p1->value(s); + + + // If this value of p1 is greater than the specified + // property value, then the current temperature is too + // high. Use it as the new upper bound. Otherwise, it + // is too low, so use it as the new lower bound. + if (pval > xval) { + tmax = t0; + phigh = pval; + } + else { + tmin = t0; + plow = pval; + } + + // Determine the new T estimate by linearly intepolation + // between the upper and lower bounds + slope = (phigh - plow)/(tmax - tmin); + dt = (xval - plow)/slope; + + // If within 100 K, terminate the search + if (fabs(dt) < 100.0) break; + + // update the T estimate + t0 = tmin + dt; + addLogEntry("new T estimate", t0); + + s.setTemperature(t0); + } + endLogGroup("Initial T Estimate"); // initial T estimate + } + + + setInitialMoles(s, elMolesGoal); + + /* + * If requested, get the initial estimate for the + * chemical potentials from the ThermoPhase object + * itself. Or else, create our own estimate. + */ + if (useThermoPhaseElementPotentials) { + bool haveEm = s.getElementPotentials(DATA_PTR(x)); + if (haveEm) { + doublereal rt = GasConstant * s.temperature(); + for (m = 0; m < m_mm; m++) { + x[m] /= rt; + } + } else { + estimateElementPotentials(s, x, elMolesGoal); + } + } else { + /* + * Calculate initial estimates of the element potentials. + * This algorithm uese the MultiPhaseEquil object's + * initialization capabilities to calculate an initial + * estimate of the mole fractions for a set of linearly + * independent component species. Then, the element + * potentials are solved for based on the chemical + * potentials of the component species. + */ + estimateElementPotentials(s, x, elMolesGoal); + } + + /* + * Do a better estimate of the element potentials. + * We have found that the current estimate may not be good + * enough to avoid drastic numerical issues associated with + * the use of a numerically generated jacobian. + * + * The Brinkley algorithm assumes a constant T, P system + * and uses a linearized analytical Jacobian that turns out + * to be very stable. + */ + int info = estimateEP_Brinkley(s, x, elMolesGoal); + if (info != 0) { + if (info == 1) { + addLogEntry("estimateEP_Brinkley didn't converge in given max interations"); + } else if (info == -3) { + addLogEntry("estimateEP_Brinkley had a singular Jacobian. Continuing anyway"); + } + } else { + setToEquilState(s, x, s.temperature()); + // Tempting -> However, nonideal is a problem. Turn on if not worried + // about nonideality and you are having problems with the main + // algorithm. + //if (XY == TP) { + // endLogGroup("ChemEquil::equilibrate"); + // return 0; + //} + } + + /* + * Install the log(temp) into the last solution unknown + * slot. + */ + x[m_mm] = log(s.temperature()); + + /* + * Setting the max and min values for x[]. Also, if element + * abundance vector is zero, setting x[] to -1000. This + * effectively zeroes out all species containing that element. + */ + vector_fp above(nvar); + vector_fp below(nvar); + for (m = 0; m < mm; m++) { + above[m] = 200.0; + below[m] = -2000.0; + if (elMolesGoal[m] < m_elemFracCutoff && m != m_eloc) x[m] = -1000.0; + } + above[mm] = log(s.maxTemp() + 1.0); + below[mm] = log(s.minTemp() - 1.0); + + vector_fp grad(nvar, 0.0); // gradient of f = F*F/2 + vector_fp oldx(nvar, 0.0); // old solution + vector_fp oldresid(nvar, 0.0); + doublereal f, oldf; + + int iter = 0; + doublereal fctr = 1.0, newval; + + goto converge; + next: + + iter++; + if (iter > 1) endLogGroup("Iteration "+int2str(iter-1)); // iteration + beginLogGroup("Iteration "+int2str(iter)); + + // compute the residual and the jacobian using the current + // solution vector + equilResidual(s, x, elMolesGoal, res_trial, xval, yval); + f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin()); + addLogEntry("Residual norm", f); + + // Compute the Jacobian matrix + equilJacobian(s, x, elMolesGoal, jac, xval, yval); + +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + sprintf(sbuf,"Jacobian matrix %d:\n", iter); writelog(sbuf); + for (m = 0; m <= m_mm; m++) { + writelog(" [ "); + for (n = 0; n <= m_mm; n++) { + sprintf(sbuf,"%10.5g ", jac(m,n)); writelog(sbuf); + } + writelog(" ]"); + char xName[32]; + if (m < m_mm) { + string nnn = eNames[m]; + sprintf(xName, "x_%-10s", nnn.c_str()); + } else { + sprintf(xName, "x_XX"); + } + if (m_eloc == m) { + sprintf(xName, "x_ELOC"); + } + if (m == m_skip) { + sprintf(xName, "x_YY"); + } + sprintf(sbuf,"%-12s", xName); writelog(sbuf); + sprintf(sbuf, " = - (%10.5g)\n", res_trial[m]); + writelog(sbuf); + } + } +#endif + + + // compute grad f = F*J + jac.leftMult(DATA_PTR(res_trial), DATA_PTR(grad)); + copy(x.begin(), x.end(), oldx.begin()); + oldf = f; + scale(res_trial.begin(), res_trial.end(), res_trial.begin(), -1.0); + + /* + * Solve the system + */ + try { + info = solve(jac, DATA_PTR(res_trial)); + } + catch (CanteraError) { + addLogEntry("Jacobian is singular."); + endLogGroup(); // iteration + endLogGroup(); // equilibrate + s.restoreState(state); + + throw CanteraError("equilibrate", + "Jacobian is singular. \nTry adding more species, " + "changing the elemental composition slightly, \nor removing " + "unused elements."); + return -3; + } + + // find the factor by which the Newton step can be multiplied + // to keep the solution within bounds. + fctr = 1.0; + for (m = 0; m < nvar; m++) { + newval = x[m] + res_trial[m]; + if (newval > above[m]) { + fctr = fmaxx( 0.0, fminn( fctr, + 0.8*(above[m] - x[m])/(newval - x[m]))); + } + else if (newval < below[m]) { + fctr = fminn(fctr, 0.8*(x[m] - below[m])/(x[m] - newval)); + } + } + if (fctr != 1.0) addLogEntry("factor to keep solution in bounds", + fctr); + + // multiply the step by the scaling factor + scale(res_trial.begin(), res_trial.end(), res_trial.begin(), fctr); + + if (!dampStep(s, oldx, oldf, grad, res_trial, + x, f, elMolesGoal , xval, yval)) + { + fail++; + if (fail > 3) { + addLogEntry("dampStep","Failed 3 times. Giving up."); + endLogGroup(); // iteration + endLogGroup(); // equilibrate + s.restoreState(state); + throw CanteraError("equilibrate", + "Cannot find an acceptable Newton damping coefficient."); + return -4; + } + } + else fail = 0; + + converge: + + // check for convergence. + equilResidual(s, x, elMolesGoal, res_trial, xval, yval); + f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin()); + doublereal xx, yy, deltax, deltay; + xx = m_p1->value(s); + yy = m_p2->value(s); + deltax = (xx - xval)/xval; + deltay = (yy - yval)/yval; + doublereal rmax = 0.0; + bool passThis = true; + for (m = 0; m < nvar; m++) { + double tval = options.relTolerance; + if (m < mm) { + if (m == m_eloc) { + tval = elMolesGoal[m] * options.relTolerance + options.absElemTol + + 1.0E-15; + } else { + tval = elMolesGoal[m] * options.relTolerance + options.absElemTol; + } + } + if (fabs(res_trial[m]) > tval) { + passThis = false; + } + } + if (iter > 0 && passThis + && fabs(deltax) < options.relTolerance + && fabs(deltay) < options.relTolerance) { + options.iterations = iter; + endLogGroup("Iteration "+int2str(iter)); // iteration + beginLogGroup("Converged solution"); + addLogEntry("Iterations",iter); + addLogEntry("Relative error in "+m_p1->symbol(),deltax); + addLogEntry("Relative error in "+m_p2->symbol(),deltay); + addLogEntry("Max residual",rmax); + beginLogGroup("Element potentials"); + doublereal rt = GasConstant* s.temperature(); + for (m = 0; m < m_mm; m++) { + m_lambda[m] = x[m]*rt; + addLogEntry("element "+ s.elementName(m), fp2str(x[m])); + } + + if (m_eloc >= 0) { + adjustEloc(s, elMolesGoal); + } + /* + * Save the calculated and converged element potentials + * to the original ThermoPhase object. + */ + s.setElementPotentials(m_lambda); + addLogEntry("Saving Element Potentials to ThermoPhase Object"); + endLogGroup("Element potentials"); + + if (s.temperature() > s.maxTemp() + 1.0 || + s.temperature() < s.minTemp() - 1.0 ) { + writelog("Warning: Temperature (" + +fp2str(s.temperature())+" K) outside " + "valid range of "+fp2str(s.minTemp())+" K to " + +fp2str(s.maxTemp())+" K\n"); + } + endLogGroup("Converged solution"); + endLogGroup("ChemEquil::equilibrate"); + return 0; + } + + // no convergence + + if (iter > options.maxIterations) { + addLogEntry("equilibrate","no convergence"); + endLogGroup("Iteration "+int2str(iter)); + endLogGroup("ChemEquil::equilibrate"); + s.restoreState(state); + throw CanteraError("equilibrate", + "no convergence in "+int2str(options.maxIterations) + +" iterations."); + return -1; + } + goto next; + } + + + /* + * dampStep: Come up with an acceptable step size. The original implementation + * employed a line search technique that enforced a reduction in the + * norm of the residual at every successful step. Unfortunately, + * this method created false convergence errors near the end of + * a significant number of steps, usually special conditions where + * there were stoichiometric constraints. + * + * This new method just does a delta damping approach, based on limiting + * the jump in the dimensionless element potentials. Mole fractions are + * limited to a factor of 2 jump in the values from this method. + * Near convergence, the delta damping gets out of the way. + */ + int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx, + double oldf, vector_fp& grad, vector_fp& step, vector_fp& x, + double& f, vector_fp& elmols, double xval, double yval ) + { + int nvar = x.size(); + int m; + double damp; + + /* + * Carry out a delta damping approach on the dimensionless element potentials. + */ + damp = 1.0; + for (m = 0; m < m_mm; m++) { + if (m == m_eloc) { + if (step[m] > 1.25) { + damp = MIN(damp, 1.25 /step[m]); + } + if (step[m] < -1.25) { + damp = MIN(damp, -1.25 / step[m]); + } + } else { + if (step[m] > 0.75) { + damp = MIN(damp, 0.75 /step[m]); + } + if (step[m] < -0.75) { + damp = MIN(damp, -0.75 / step[m]); + } + } + } + + /* + * Update the solution unknown + */ + for (m = 0; m < nvar; m++) { + x[m] = oldx[m] + damp * step[m]; + } +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + sprintf(sbuf, "Solution Unknowns: damp = %g\n", damp); writelog(sbuf); + writelog(" X_new X_old Step\n"); + for (m = 0; m < nvar; m++) { + sprintf(sbuf," % -10.5g % -10.5g % -10.5g\n", x[m], oldx[m], step[m]); + writelog(sbuf); + } + } +#endif + return 1; + } + + + /** + * Evaluates the residual vector F, of length mm + */ + void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x, + const vector_fp& elmFracGoal, vector_fp& resid, + doublereal xval, doublereal yval) + { + beginLogGroup("ChemEquil::equilResidual"); + int n, m; + doublereal xx, yy; + doublereal temp = exp(x[m_mm]); + setToEquilState(s, x, temp); + + // residuals are the total element moles + vector_fp& elmFrac = m_elementmolefracs; + for (n = 0; n < m_mm; n++) { + m = m_orderVectorElements[n]; + // drive element potential for absent elements to -1000 + if (elmFracGoal[m] < m_elemFracCutoff && m != m_eloc) { + resid[m] = x[m] + 1000.0; + } else if (n >= m_nComponents) { + resid[m] = x[m]; + } else { + /* + * Change the calculation for small element number, using + * L'Hopital's rule. + * The log formulation is unstable. + */ + if (elmFracGoal[m] < 1.0E-10 || elmFrac[m] < 1.0E-10 || m == m_eloc) { + resid[m] = elmFracGoal[m] - elmFrac[m]; + } else { + resid[m] = log( (1.0 + elmFracGoal[m]) / (1.0 + elmFrac[m]) ); + } + } + addLogEntry(s.elementName(m),fp2str(elmFrac[m])+" (" + +fp2str(elmFracGoal[m])+")"); + } + +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0 && !m_doResPerturb) { + writelog("Residual: ElFracGoal ElFracCurrent Resid\n"); + for (n = 0; n < m_mm; n++) { + sprintf(sbuf," % -14.7E % -14.7E % -10.5E\n", + elmFracGoal[n], elmFrac[n], resid[n]); + writelog(sbuf); + } + } +#endif + + xx = m_p1->value(s); + yy = m_p2->value(s); + resid[m_mm] = xx/xval - 1.0; + resid[m_skip] = yy/yval - 1.0; + string xstr = fp2str(xx)+" ("+fp2str(xval)+")"; + addLogEntry(m_p1->symbol(), xstr); + string ystr = fp2str(yy)+" ("+fp2str(yval)+")"; + addLogEntry(m_p2->symbol(), ystr); + endLogGroup("ChemEquil::equilResidual"); + +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0 && !m_doResPerturb) { + writelog(" Goal Xvalue Resid\n"); + sprintf(sbuf," XX : % -14.7E % -14.7E % -10.5E\n", xval, xx, resid[m_mm]); + writelog(sbuf); + sprintf(sbuf," YY(%1d): % -14.7E % -14.7E % -10.5E\n", m_skip, yval, yy, resid[m_skip]); + writelog(sbuf); + } +#endif + } + + + //-------------------- Jacobian evaluation --------------------------- + + void ChemEquil::equilJacobian(thermo_t& s, vector_fp& x, + const vector_fp& elmols, DenseMatrix& jac, + doublereal xval, doublereal yval) + { + beginLogGroup("equilJacobian"); + int len = x.size(); + vector_fp& r0 = m_jwork1; + vector_fp& r1 = m_jwork2; + r0.resize(len); + r1.resize(len); + int n, m; + doublereal rdx, dx, xsave, dx2; + doublereal atol = 1.e-10; + + equilResidual(s, x, elmols, r0, xval, yval); + + m_doResPerturb = false; + for (n = 0; n < len; n++) { + xsave = x[n]; + dx = atol; + dx2 = fabs(xsave) * 1.0E-7; + if (dx2 > dx) dx = dx2; + x[n] = xsave + dx; + dx = x[n] - xsave; + rdx = 1.0/dx; + + // calculate perturbed residual + + equilResidual(s, x, elmols, r1, xval, yval); + + // compute nth column of Jacobian + + for (m = 0; m < len; m++) { + jac(m, n) = (r1[m] - r0[m])*rdx; + } + x[n] = xsave; + } + m_doResPerturb = false; + endLogGroup("equilJacobian"); + } + + /** + * Given a vector of dimensionless element abundances, + * this routine calculates the moles of the elements and + * the moles of the species. + * Input + * -------- + * x[m] = current dimensionless element potentials.. + */ + double ChemEquil::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 k, m; + double n_t_calc = 0.0; + double tmp; + /* + * Calculate the activity coefficients of the solution, at the + * previous solution state. + */ + vector_fp actCoeff(m_kk, 1.0); + s.setMoleFractions(DATA_PTR(Xmol_i_calc)); + s.setPressure(pressureConst); + s.getActivityCoefficients(DATA_PTR(actCoeff)); + + for (k = 0; k < m_kk; k++) { + tmp = - (m_muSS_RT[k] + log(actCoeff[k])); + for (m = 0; m < m_mm; m++) { + tmp += nAtoms(k,m) * x[m]; + } + if (tmp > 100.) tmp = 100.; + if (tmp < -300.) { + n_i_calc[k] = 0.0; + } else { + n_i_calc[k] = n_t * exp(tmp); + } + n_t_calc += n_i_calc[k]; + } + for (m = 0; m < m_mm; m++) { + eMolesCalc[m] = 0.0; + for (k = 0; k < m_kk; k++) { + eMolesCalc[m] += nAtoms(k,m) * n_i_calc[k]; + } + } + return n_t_calc; + } + + /** + * Do a calculation of the element potentials using + * the Brinkley method, p. 129 Smith and Missen. + * + * We have found that the previous estimate may not be good + * enough to avoid drastic numerical issues associated with + * the use of a numerically generated jacobian used in the + * main algorithm. + * + * The Brinkley algorithm, here, assumes a constant T, P system + * and uses a linearized analytical Jacobian that turns out + * to be very stable even given bad initial guesses. + * + * The pressure and temperature to be used are in the + * ThermoPhase object input into the routine. + * + * The initial guess for the element potentials + * used by this routine is taken from the + * input vector, x. + * + * elMoles is the input element abundance vector to be matched. + * + * Nonideal phases are handled in principle. This is done by + * calculating the activity coefficients and adding them + * into the formula in the correct position. However, + * these are treated as a rhs contribution only. Therefore, + * convergence might be a problem. This has not been tested. + * Also molality based unit systems aren't handled. + * + * On return, int return value contains the success code: + * 0 - successful + * 1 - unsuccessful, max num iterations exceeded + * -3 - unsuccessful, singular jacobian + * + * NOTE: update for activity coefficients. + */ + int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, + vector_fp& elMoles) { + /* + * Before we do anything, we will save the state of the solution. + * Then, if things go drastically wrong, we will restore the + * saved state. + */ + vector_fp state; + s.saveState(state); + double tmp, sum; + bool modifiedMatrix = false; + int neq = m_mm+1; + int retn = 1; + int m, n, k, info, im; + DenseMatrix a1(neq, neq, 0.0); + vector_fp b(neq, 0.0); + vector_fp n_i(m_kk,0.0); + vector_fp n_i_calc(m_kk,0.0); + vector_fp actCoeff(m_kk, 1.0); + + vector_fp Xmol_i_calc(m_kk,0.0); + double beta = 1.0; + + s.getMoleFractions(DATA_PTR(n_i)); + double pressureConst = s.pressure(); + copy(n_i.begin(), n_i.end(), Xmol_i_calc.begin()); + + vector_fp x_old(m_mm+1, 0.0); + vector_fp resid(m_mm+1, 0.0); + vector_int lumpSum(m_mm+1, 0); + + + /* + * Get the nondimensional Gibbs functions for the species + * at their standard states of solution at the current T and P + * of the solution. + */ + s.getGibbs_RT(DATA_PTR(m_muSS_RT)); + + + vector_fp eMolesCalc(m_mm, 0.0); + vector_fp eMolesFix(m_mm, 0.0); + double elMolesTotal = 0.0; + for (m = 0; m < m_mm; m++) { + elMolesTotal += elMoles[m]; + for (k = 0; k < m_kk; k++) { + eMolesFix[m] += nAtoms(k,m) * n_i[k]; + } + } + + for (m = 0; m < m_mm; m++) { + if (x[m] > 50.0) { + x[m] = 50.; + } + if (elMoles[m] > 1.0E-70) { + if (x[m] < -100) { + x[m] = -100.; + } + } else { + if (x[m] < -1000.) { + x[m] = -1000.; + } + } + } + + + double n_t = 0.0; + double sum2 = 0.0; + double nAtomsMax = 1.0; + s.setMoleFractions(DATA_PTR(Xmol_i_calc)); + s.setPressure(pressureConst); + s.getActivityCoefficients(DATA_PTR(actCoeff)); + for (k = 0; k < m_kk; k++) { + tmp = - (m_muSS_RT[k] + log(actCoeff[k])); + sum2 = 0.0; + for (m = 0; m < m_mm; m++) { + sum = nAtoms(k,m); + tmp += sum * x[m]; + sum2 += sum; + if (sum2 > nAtomsMax) { + nAtomsMax = sum2; + } + } + if (tmp > 100.) { + n_t += 2.8E43; + } else { + n_t += exp(tmp); + } + } + + +#ifdef DEBUG_HKM + const vector& eNames = s.elementNames(); + if (ChemEquil_print_lvl > 0) { + writelog("estimateEP_Brinkley::\n\n"); + double temp = s.temperature(); + double pres = s.pressure(); + sprintf(sbuf, "temp = %g\n", temp); writelog(sbuf); + sprintf(sbuf, "pres = %g\n", pres); writelog(sbuf); + writelog("Initial mole numbers and mu_SS:\n"); + writelog(" Name MoleNum mu_SS actCoeff\n"); + for (k = 0; k < m_kk; k++) { + string nnn = s.speciesName(k); + sprintf(sbuf,"%15s %13.5g %13.5g %13.5g\n", + nnn.c_str(), n_i[k], m_muSS_RT[k], actCoeff[k]); + writelog(sbuf); + } + sprintf(sbuf,"Initial n_t = %10.5g\n", n_t); writelog(sbuf); + writelog("Comparison of Goal Element Abundance with Initial Guess:\n"); + writelog(" eName eCurrent eGoal\n"); + for (m = 0; m < m_mm; m++) { + string nnn = s.elementName(m); + sprintf(sbuf,"%5s %13.5g %13.5g\n",nnn.c_str(), eMolesFix[m], elMoles[m]); + writelog(sbuf); + } + } +#endif + for (m = 0; m < m_mm; m++) { + if (m != m_eloc) { + if (elMoles[m] <= options.absElemTol) { + x[m] = -200.; + } + } + } + /* + * ------------------------------------------------------------------- + * Main Loop. + */ + for (int iter = 0; iter < 20* options.maxIterations; iter++) { + /* + * Save the old solution + */ + for (m = 0; m < m_mm; m++) { + x_old[m] = x[m]; + } + x_old[m_mm] = n_t; + /* + * Calculate the mole numbers of species + */ +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + sprintf(sbuf, "START ITERATION %d:\n", iter); writelog(sbuf); + } +#endif + /* + * Calculate the mole numbers of species and elements. + */ + double n_t_calc = calcEmoles(s, x, n_t, Xmol_i_calc, eMolesCalc, n_i_calc, + pressureConst); + for (k = 0; k < m_kk; k++) { + Xmol_i_calc[k] = n_i_calc[k]/n_t_calc; + } + +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + writelog(" Species: Calculated_Moles Calculated_Mole_Fraction\n"); + for (k = 0; k < m_kk; k++) { + string nnn = s.speciesName(k); + sprintf(sbuf,"%15s: %10.5g %10.5g\n", nnn.c_str(), n_i_calc[k], Xmol_i_calc[k]); + writelog(sbuf); + } + sprintf(sbuf,"%15s: %10.5g\n", "Total Molar Sum", n_t_calc); + writelog(sbuf); + sprintf(sbuf,"(iter %d) element moles bal: Goal Calculated\n", iter); + writelog(sbuf); + for (m = 0; m < m_mm; m++) { + string nnn = eNames[m]; + sprintf(sbuf," %8s: %10.5g %10.5g \n", nnn.c_str(), elMoles[m], eMolesCalc[m]); + writelog(sbuf); + } + } +#endif + + double nCutoff; + + bool normalStep = true; + /* + * Decide if we are to do a normal step or a modified step + */ + int iM = -1; + for (m = 0; m < m_mm; m++) { + if (elMoles[m] > 0.001 * elMolesTotal) { + if (eMolesCalc[m] > 1000. * elMoles[m]) { + normalStep = false; + iM = m; + } + if (1000 * eMolesCalc[m] < elMoles[m]) { + normalStep = false; + iM = m; + } + } + } +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + if (!normalStep) { + sprintf(sbuf," NOTE: iter(%d) Doing an abnormal step due to row %d\n", iter, iM); + writelog(sbuf); + } + } +#endif + if (!normalStep) { + beta = 1.0; + resid[m_mm] = 0.0; + for (im = 0; im < m_mm; im++) { + m = m_orderVectorElements[im]; + resid[m] = 0.0; + if (im < m_nComponents) { + if (elMoles[m] > 0.001 * elMolesTotal) { + if (eMolesCalc[m] > 1000. * elMoles[m]) { + resid[m] = -0.5; + resid[m_mm] -= 0.5; + } + if (1000 * eMolesCalc[m] < elMoles[m]) { + resid[m] = 0.5; + resid[m_mm] += 0.5; + } + } + } + } + if (n_t < (elMolesTotal / nAtomsMax)) { + if (resid[m_mm] < 0.0) { + resid[m_mm] = 0.1; + } + } else if (n_t > elMolesTotal) { + if (resid[m_mm] > 0.0) { + resid[m_mm] = 0.0; + } + } + goto updateSolnVector; + } + + + /* + * Determine whether the matrix should be dumbed down because + * the coefficient matrix of species (with significant concentrations) + * is rank deficient. + * + * The basic idea is that at any time during the calculation only a + * small subset of species with sufficient concentration matters. + * If the rank of the element coefficient matrix for that subset of species + * is less than the number of elements, then the matrix created by + * the Brinkley method below may become singular. + * + * The logic below looks for obvious cases where the current element + * coefficient matrix is rank deficient. + * + * The way around rank-deficiency is to lump-sum the corresponding row + * of the matrix. Note, lump-summing seems to work very well in terms of + * its stability properties, i.e., it heads in the right direction, + * albeit with lousy convergence rates. + * + * NOTE: This probably should be extended to a full blown Gauss-Jordon + * factorization scheme in the future. For Example + * the scheme below would fail for the set: HCl NH4Cl, NH3. + * Hopefully, it's caught by the equal rows logic below. + */ + for (m = 0; m < m_mm; m++) { + lumpSum[m] = 1; + } + + nCutoff = 1.0E-9 * n_t_calc; +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + writelog(" Lump Sum Elements Calculation: \n"); + } +#endif + for (m = 0; m < m_mm; m++) { + int kMSp = -1; + int kMSp2 = -1; + int nSpeciesWithElem = 0; + for (k = 0; k < m_kk; k++) { + if (n_i_calc[k] > nCutoff) { + if (fabs(nAtoms(k,m)) > 0.001) { + nSpeciesWithElem++; + if (kMSp != -1) { + kMSp2 = k; + double factor = fabs(nAtoms(kMSp,m) / nAtoms(kMSp2,m)); + for (n = 0; n < m_mm; n++) { + if (fabs(factor * nAtoms(kMSp2,n) - nAtoms(kMSp,n)) > 1.0E-8) { + lumpSum[m] = 0; + break; + } + } + } else { + kMSp = k; + } + } + } + } +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + string nnn = eNames[m]; + sprintf(sbuf," %5s %3d : %5d %5d\n",nnn.c_str(), lumpSum[m], kMSp, kMSp2); + writelog(sbuf); + } +#endif + } + + /* + * Formulate the matrix. + */ + for (im = 0; im < m_mm; im++) { + m = m_orderVectorElements[im]; + if (im < m_nComponents) { + for (n = 0; n < m_mm; n++) { + a1(m,n) = 0.0; + for (k = 0; k < m_kk; k++) { + a1(m,n) += nAtoms(k,m) * nAtoms(k,n) * n_i_calc[k]; + } + } + a1(m,m_mm) = eMolesCalc[m]; + a1(m_mm, m) = eMolesCalc[m]; + } else { + for (n = 0; n <= m_mm; n++) { + a1(m,n) = 0.0; + } + a1(m,m) = 1.0; + } + } + a1(m_mm, m_mm) = 0.0; + + /* + * Formulate the residual, resid, and the estimate for the convergence criteria, sum + */ + sum = 0.0; + for (im = 0; im < m_mm; im++) { + m = m_orderVectorElements[im]; + if (im < m_nComponents) { + resid[m] = elMoles[m] - eMolesCalc[m]; + } else { + resid[m] = 0.0; + } + /* + * For equations with positive and negative coefficients, (electronic charge), + * we must mitigate the convergence criteria by a condition limited by + * finite precision of inverting a matrix. + * Other equations with just positive coefficients aren't limited by this. + */ + if (m == m_eloc) { + tmp = resid[m] / (elMoles[m] + elMolesTotal*1.0E-6 + options.absElemTol); + } else { + tmp = resid[m] / (elMoles[m] + options.absElemTol); + } + sum += tmp * tmp; + } + + for (m = 0; m < m_mm; m++) { + if (a1(m,m) < 1.0E-50) { +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + sprintf(sbuf," NOTE: Diagonalizing the analytical Jac row %d\n", m); + writelog(sbuf); + } +#endif + for (n = 0; n < m_mm; n++) { + a1(m,n) = 0.0; + } + a1(m,m) = 1.0; + if (resid[m] > 0.0) { + resid[m] = 1.0; + } else if (resid[m] < 0.0) { + resid[m] = -1.0; + } else { + resid[m] = 0.0; + } + } + } + + + resid[m_mm] = n_t - n_t_calc; + +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + writelog("Matrix:\n"); + for (m = 0; m <= m_mm; m++) { + writelog(" ["); + for (n = 0; n <= m_mm; n++) { + sprintf(sbuf," %10.5g", a1(m,n)); writelog(sbuf); + } + sprintf(sbuf,"] = %10.5g\n", resid[m]); writelog(sbuf); + } + } +#endif + + tmp = resid[m_mm] /(n_t + 1.0E-15); + sum += tmp * tmp; +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + sprintf(sbuf,"(it %d) Convergence = %g\n", iter, sum); + writelog(sbuf); + } +#endif + /* + * Insist on 20x accuracy compared to the top routine. + * There are instances, for ill-conditioned or + * singular matrices where this is needed to move + * the system to a point where the matrices aren't + * singular. + */ + if (sum < 0.05 * options.relTolerance) { + retn = 0; + goto exit; + } + + /* + * Row Sum scaling + */ + for (m = 0; m <= m_mm; m++) { + tmp = 0.0; + for (n = 0; n <= m_mm; n++) { + tmp += fabs(a1(m,n)); + } + if (m < m_mm && tmp < 1.0E-30) { +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + sprintf(sbuf," NOTE: Diagonalizing row %d\n", m); + writelog(sbuf); + } +#endif + for (n = 0; n <= m_mm; n++) { + if (n != m) { + a1(m,n) = 0.0; + a1(n,m) = 0.0; + } + } + } + tmp = 1.0/tmp; + for (n = 0; n <= m_mm; n++) { + a1(m,n) *= tmp; + } + resid[m] *= tmp; + } + +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + writelog("Row Summed Matrix:\n"); + for (m = 0; m <= m_mm; m++) { + writelog(" ["); + for (n = 0; n <= m_mm; n++) { + sprintf(sbuf," %10.5g", a1(m,n)); writelog(sbuf); + } + sprintf(sbuf,"] = %10.5g\n", resid[m]); writelog(sbuf); + } + } +#endif + + /* + * Next Step: We have row-summed the equations. + * However, there are some degenerate cases where two + * rows will be multiplies of each other in terms of + * 0 < m, 0 < m part of the matrix. This occurs on a case + * by case basis, and depends upon the current state of the + * element potential values, which affect the concentrations + * of species. + * So, the way we have found to eliminate this problem is to + * lump-sum one of the rows of the matrix, except for the + * last column, and stick it all on the diagonal. + * Then, we at least have a non-singular matrix, and the + * modified equation moves the corresponding unknown in the + * correct direction. + * The previous row-sum operation has made the identification + * of identical rows much simpler. + * + * Note at least 6E-4 is necessary for the comparison. + * I'm guessing 1.0E-3. If two rows are anywhere close to being + * equivalent, the algorithm can get stuck in an oscillatory mode. + */ + modifiedMatrix = false; + for (m = 0; m < m_mm; m++) { + int sameAsRow = -1; + for (int im = 0; im < m; im++) { + bool theSame = true; + for (n = 0; n < m_mm; n++) { + if (fabs(a1(m,n) - a1(im,n)) > 1.0E-7) { + theSame = false; + break; + } + } + if (theSame) { + sameAsRow = im; + } + } + if (sameAsRow >= 0 || lumpSum[m]) { +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + if (lumpSum[m]) { + sprintf(sbuf,"Lump summing row %d, due to rank deficiency analysis\n", m); + writelog(sbuf); + } else if (sameAsRow >= 0) { + sprintf(sbuf,"Identified that rows %d and %d are the same\n", m, sameAsRow); + writelog(sbuf); + } + } +#endif + modifiedMatrix = true; + for (n = 0; n < m_mm; n++) { + if (n != m) { + a1(m,m) += fabs(a1(m,n)); + a1(m,n) = 0.0; + } + } + } + } + +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0 && modifiedMatrix) { + writelog("Row Summed, MODIFIED Matrix:\n"); + for (m = 0; m <= m_mm; m++) { + writelog(" ["); + for (n = 0; n <= m_mm; n++) { + sprintf(sbuf," %10.5g", a1(m,n)); writelog(sbuf); + } + sprintf(sbuf,"] = %10.5g\n", resid[m]); writelog(sbuf); + } + } +#endif + + try { + info = solve(a1, DATA_PTR(resid)); + } + catch (CanteraError) { + addLogEntry("estimateEP_Brinkley:Jacobian is singular."); +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + writelog("Matrix is SINGULAR.ERROR\n"); + } +#endif + s.restoreState(state); + throw CanteraError("equilibrate:estimateEP_Brinkley()", + "Jacobian is singular. \nTry adding more species, " + "changing the elemental composition slightly, \nor removing " + "unused elements."); + return -3; + } + + /* + * Figure out the damping coefficient: Use a delta damping + * coefficient formulation: magnitude of change is capped + * to exp(1). + */ + beta = 1.0; + for (m = 0; m < m_mm; m++) { + if (resid[m] > 1.0) { + double betat = 1.0 / resid[m]; + if (betat < beta) { + beta = betat; + } + } + if (resid[m] < -1.0) { + double betat = -1.0 / resid[m]; + if (betat < beta) { + beta = betat; + } + } + } +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + if (beta != 1.0) { + sprintf(sbuf,"(it %d) Beta = %g\n", iter, beta); writelog(sbuf); + } + } +#endif + + /* + * Update the solution vector + */ + updateSolnVector: + for (m = 0; m < m_mm; m++) { + x[m] += beta * resid[m]; + } + n_t *= exp(beta * resid[m_mm]); + + +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + sprintf(sbuf,"(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter); + writelog(sbuf); + for (m = 0; m < m_mm; m++) { + string eee = eNames[m]; + sprintf(sbuf," %5s %10.5g %10.5g %10.5g\n", eee.c_str(), x_old[m], x[m], resid[m]); + writelog(sbuf); + } + sprintf(sbuf," n_t %10.5g %10.5g %10.5g \n", x_old[m_mm], n_t, exp(resid[m_mm])); + writelog(sbuf); + } +#endif + } + exit: +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + double temp = s.temperature(); + double pres = s.pressure(); + + if (retn == 0) { + sprintf(sbuf," ChemEquil::estimateEP_Brinkley() SUCCESS: equilibrium found at T = %g, Pres = %g\n", + temp, pres); + writelog(sbuf); + } else { + sprintf(sbuf," ChemEquil::estimateEP_Brinkley() FAILURE: equilibrium not found at T = %g, Pres = %g\n", + temp, pres); + writelog(sbuf); + } + } +#endif + return retn; + } + + + /* + * + */ + void ChemEquil::adjustEloc(thermo_t &s, vector_fp & elMolesGoal) { + if (m_eloc < 0) return; + if (fabs(elMolesGoal[m_eloc]) > 1.0E-20) return; + s.getMoleFractions(DATA_PTR(m_molefractions)); + int k; + +#ifdef DEBUG_HKM + int maxPosEloc = -1; + int maxNegEloc = -1; + double maxPosVal = -1.0; + double maxNegVal = -1.0; + if (ChemEquil_print_lvl > 0) { + for (k = 0; k < m_kk; k++) { + if (nAtoms(k,m_eloc) > 0.0) { + if (m_molefractions[k] > maxPosVal && m_molefractions[k] > 0.0) { + maxPosVal = m_molefractions[k]; + maxPosEloc = k; + } + } + if (nAtoms(k,m_eloc) < 0.0) { + if (m_molefractions[k] > maxNegVal && m_molefractions[k] > 0.0) { + maxNegVal = m_molefractions[k]; + maxNegEloc = k; + } + } + } + } +#endif + + double sumPos = 0.0; + double sumNeg = 0.0; + for (k = 0; k < m_kk; k++) { + if (nAtoms(k,m_eloc) > 0.0) { + sumPos += nAtoms(k,m_eloc) * m_molefractions[k]; + } + if (nAtoms(k,m_eloc) < 0.0) { + sumNeg += nAtoms(k,m_eloc) * m_molefractions[k]; + } + } + sumNeg = - sumNeg; + + if (sumPos >= sumNeg) { + if ( sumPos <= 0.0) return; + double factor = (elMolesGoal[m_eloc] + sumNeg) / sumPos; +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + if (factor < 0.9999999999) { + string nnn = s.speciesName(maxPosEloc); + sprintf(sbuf,"adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n", + nnn.c_str(), + m_molefractions[maxPosEloc], m_molefractions[maxPosEloc]*factor); + writelog(sbuf); + } + } +#endif + for (k = 0; k < m_kk; k++) { + if (nAtoms(k,m_eloc) > 0.0) { + m_molefractions[k] *= factor; + } + } + } else { + double factor = (-elMolesGoal[m_eloc] + sumPos) / sumNeg; +#ifdef DEBUG_HKM + if (ChemEquil_print_lvl > 0) { + if (factor < 0.9999999999) { + string nnn = s.speciesName(maxNegEloc); + sprintf(sbuf,"adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n", + nnn.c_str(), + m_molefractions[maxNegEloc], m_molefractions[maxNegEloc]*factor); + writelog(sbuf); + } + } +#endif + for (k = 0; k < m_kk; k++) { + if (nAtoms(k,m_eloc) < 0.0) { + m_molefractions[k] *= factor; + } + } + } + + s.setMoleFractions(DATA_PTR(m_molefractions)); + s.getMoleFractions(DATA_PTR(m_molefractions)); + + } + +} // namespace diff --git a/Cantera/src/equil/ChemEquil.h b/Cantera/src/equil/ChemEquil.h new file mode 100755 index 000000000..8ce080eb4 --- /dev/null +++ b/Cantera/src/equil/ChemEquil.h @@ -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 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 *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 diff --git a/Cantera/src/equil/MultiPhase.cpp b/Cantera/src/equil/MultiPhase.cpp new file mode 100644 index 000000000..fbd79c4cb --- /dev/null +++ b/Cantera/src/equil/MultiPhase.cpp @@ -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 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; + } + } + +} + diff --git a/Cantera/src/equil/MultiPhase.h b/Cantera/src/equil/MultiPhase.h new file mode 100644 index 000000000..ca9440eab --- /dev/null +++ b/Cantera/src/equil/MultiPhase.h @@ -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_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 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 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 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 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 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 diff --git a/Cantera/src/equil/MultiPhaseEquil.cpp b/Cantera/src/equil/MultiPhaseEquil.cpp new file mode 100644 index 000000000..5fd7e1e3f --- /dev/null +++ b/Cantera/src/equil/MultiPhaseEquil.cpp @@ -0,0 +1,832 @@ +#include "MultiPhaseEquil.h" +#include "MultiPhase.h" +#include "sort.h" +#include "global.h" + +#include +#include +using namespace std; + + +namespace Cantera { + + const doublereal TINY = 1.0e-20; + + /// Used to print reaction equations. Given a stoichiometric + /// coefficient 'nu' and a chemical symbol 'sym', return a string + /// for this species in the reaction. + /// @param first if this is false, then a " + " string will be + /// added to the beginning of the string. + /// @param nu Stoichiometric coefficient. May be positive or negative. 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; + } +} diff --git a/Cantera/src/equil/MultiPhaseEquil.h b/Cantera/src/equil/MultiPhaseEquil.h new file mode 100644 index 000000000..07a979a0c --- /dev/null +++ b/Cantera/src/equil/MultiPhaseEquil.h @@ -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 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 m_solnrxn; + bool m_force; + }; + +} + + +#endif diff --git a/Cantera/src/equil/PropertyCalculator.h b/Cantera/src/equil/PropertyCalculator.h new file mode 100755 index 000000000..be69dbe29 --- /dev/null +++ b/Cantera/src/equil/PropertyCalculator.h @@ -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 PropertyCalculator { + public: + virtual ~PropertyCalculator(){} + virtual doublereal value(const M& s) =0; + virtual string symbol() =0; + }; + + template + class EnthalpyCalculator : public PropertyCalculator { + public: + virtual doublereal value(const M& s) { + return s.enthalpy_mass(); + } + virtual string symbol() { return "H"; } + }; + + template + class EntropyCalculator : public PropertyCalculator { + public: + virtual doublereal value(const M& s) { + return s.entropy_mass(); + } + virtual string symbol() { return "S"; } + }; + + template + class TemperatureCalculator : public PropertyCalculator { + public: + virtual doublereal value(const M& s) { + return s.temperature(); + } + virtual string symbol() { return "T"; } + }; + + template + class PressureCalculator : public PropertyCalculator { + public: + virtual doublereal value(const M& s) { + return s.pressure(); + } + virtual string symbol() { return "P"; } + }; + + template + class DensityCalculator : public PropertyCalculator { + public: + virtual doublereal value(const M& s) { + return s.density(); + } + virtual string symbol() { return "V"; } + }; + + template + class IntEnergyCalculator : public PropertyCalculator { + public: + virtual doublereal value(const M& s) { + return s.intEnergy_mass(); + } + virtual string symbol() { return "U"; } + }; +} + +#endif + diff --git a/Cantera/src/equil/equil.h b/Cantera/src/equil/equil.h new file mode 100644 index 000000000..cd43ba7c5 --- /dev/null +++ b/Cantera/src/equil/equil.h @@ -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 diff --git a/Cantera/src/equil/equilibrate.cpp b/Cantera/src/equil/equilibrate.cpp new file mode 100644 index 000000000..e0769d809 --- /dev/null +++ b/Cantera/src/equil/equilibrate.cpp @@ -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; + } +}