diff --git a/Cantera/src/BasisOptimize.cpp b/Cantera/src/BasisOptimize.cpp new file mode 100644 index 000000000..972aa04a0 --- /dev/null +++ b/Cantera/src/BasisOptimize.cpp @@ -0,0 +1,531 @@ +/** + * @file basopt.cpp + * + * $Author$ + * $Date$ + * $Revision$ + */ + +#include "ct_defs.h" +#include "ThermoPhase.h" +#include "MultiPhase.h" + +using namespace Cantera; +using namespace std; +#ifdef DEBUG_HKM +extern int debug_print_lvl; +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 + * orderVectorElement + * orderVectorSpecies + * + * 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 + if (debug_print_lvl >= 1) { + printf(" "); for(i=0; i<77; i++) printf("-"); printf("\n"); + printf(" --- Subroutine BASOPT called to "); + printf("calculate the number of components and "); + printf("evaluate the formation matrix\n"); + if (debug_print_lvl > 0) { + printf(" ---\n"); + + printf(" --- Formula Matrix used in BASOPT calculation\n"); + printf(" --- Species | Order | "); + for (j = 0; j < ne; j++) { + jj = orderVectorElements[j]; + printf(" "); + ename = mphase->elementName(jj); + print_stringTrunc(ename.c_str(), 4, 1); + printf("(%1d)", j); + } + printf("\n"); + for (k = 0; k < nspecies; k++) { + kk = orderVectorSpecies[k]; + printf(" --- "); + sname = mphase->speciesName(kk); + print_stringTrunc(sname.c_str(), 11, 1); + printf(" | %4d |", k); + for (j = 0; j < ne; j++) { + jj = orderVectorElements[j]; + double num = mphase->nAtoms(kk,jj); + printf("%6.1g ", num); + } + printf("\n"); + } + printf(" --- \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. + */ + k = amax(DATA_PTR(molNum), jr, nspecies); + if (molNum[k] == 0.0) *usedZeroedSpecies = true; + /* + * If the largest molNum is negative, then we are done. + */ + if (molNum[k] == 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. + */ + molNum[k] = 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(k,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 (debug_print_lvl >= 1) { + kk = orderVectorSpecies[k]; + sname = mphase->speciesName(kk); + printf(" --- %-12.12s", sname.c_str()); + jj = orderVectorSpecies[jr]; + ename = mphase->speciesName(jj); + printf("(%9.2g) replaces %-12.12s", molNum[kk], ename.c_str()); + printf("(%9.2g) as component %3d\n", molNum[jj], jr); + } +#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) { + printf("ERROR: mlequ returned an error condition\n"); + throw CanteraError("basopt", "mlequ returned an error condition"); + } + +#ifdef DEBUG_HKM + if (debug_print_lvl >= 1) { + printf(" ---\n"); + printf(" --- Number of Components = %d\n", nComponents); + printf(" --- Formula Matrix:\n"); + printf(" --- Components: "); + for (k = 0; k < nComponents; k++) { + kk = orderVectorSpecies[k]; + printf(" %3d (%3d) ", k, kk); + } + printf("\n --- Components Moles: "); + for (k = 0; k < nComponents; k++) { + kk = orderVectorSpecies[k]; + printf("%-11.3g", molNumBase[kk]); + } + printf("\n --- NonComponent | Moles | "); + for (i = 0; i < nComponents; i++) { + kk = orderVectorSpecies[i]; + sname = mphase->speciesName(kk); + printf("%-11.10s", sname.c_str()); + } + printf("\n"); + + for (i = 0; i < nNonComponents; i++) { + k = i + nComponents; + kk = orderVectorSpecies[k]; + printf(" --- %3d (%3d) ", k, kk); + sname = mphase->speciesName(kk); + printf("%-10.10s", sname.c_str()); + printf("|%10.3g|", molNumBase[kk]); + /* + * Print the negative of formRxnMatrix[]; it's easier to interpret. + */ + for (j = 0; j < nComponents; j++) { + printf(" %6.2f", - formRxnMatrix[j + i * ne]); + } + printf("\n"); + } + printf(" "); for (i=0; i<77; i++) printf("-"); printf("\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++) { + printf("%c", str[i]); + } + } 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++) printf(" "); + } + printf("%s", str); + if (rs != 0) { + for (i = 0; i < rs; i++) printf(" "); + } + } +} +#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; + } + printf("vcs_mlequ ERROR: Encountered a zero column: %d\n", i); + 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() *************************************************************/ diff --git a/Cantera/src/ChemEquil.cpp b/Cantera/src/ChemEquil.cpp index 77a719b51..f86dda9e5 100755 --- a/Cantera/src/ChemEquil.cpp +++ b/Cantera/src/ChemEquil.cpp @@ -103,6 +103,7 @@ namespace Cantera { 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); // set up elemental composition matrix @@ -1016,6 +1017,50 @@ namespace Cantera { 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) { + 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.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. @@ -1068,28 +1113,29 @@ namespace Cantera { int m, n, k, info; DenseMatrix a1(neq, neq, 0.0); vector_fp b(neq, 0.0); - vector_fp muSS_RT(m_kk, 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 muSS_RT_mod(m_kk, 0.0); + + vector_fp Xmol_i_calc(m_kk,0.0); double beta = 1.0; s.getMoleFractions(DATA_PTR(n_i)); + 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(muSS_RT)); - copy(muSS_RT.begin(), muSS_RT.end(), muSS_RT_mod.begin()); - + 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; @@ -1119,10 +1165,10 @@ namespace Cantera { double n_t = 0.0; double sum2 = 0.0; double nAtomsMax = 1.0; - s.setMoleFractions(DATA_PTR(x)); + s.setMoleFractions(DATA_PTR(Xmol_i_calc)); s.getActivityCoefficients(DATA_PTR(actCoeff)); for (k = 0; k < m_kk; k++) { - tmp = - (muSS_RT[k] + log(actCoeff[k])); + tmp = - (m_muSS_RT[k] + log(actCoeff[k])); sum2 = 0.0; for (m = 0; m < m_mm; m++) { sum = nAtoms(k,m); @@ -1191,44 +1237,21 @@ namespace Cantera { printf("START ITERATION %d:\n", iter); } #endif - - double n_t_calc = 0.0; - s.setMoleFractions(DATA_PTR(x)); - s.getActivityCoefficients(DATA_PTR(actCoeff)); - sum2 = 0.0; + /* + * Calculate the mole numbers of species and elements. + */ + double n_t_calc = calcEmoles(s, x, n_t, Xmol_i_calc, eMolesCalc, n_i_calc); for (k = 0; k < m_kk; k++) { - tmp = - (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]; -#ifdef DEBUG_HKM_EPEQUIL - if (debug_prnt_lvl > 0) { - string nnn = s.speciesName(k); - printf("%15s: %10.5g (%10.5g)\n", nnn.c_str(), - n_i_calc[k], tmp); - } -#endif + Xmol_i_calc[k] = n_i_calc[k]/n_t_calc; } + #ifdef DEBUG_HKM_EPEQUIL if (debug_prnt_lvl > 0) { - printf("%15s: %10.5g\n", "Total Molar Sum", n_t_calc); - } -#endif - 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]; + string nnn = s.speciesName(k); + printf("%15s: %10.5g %10.5g\n", nnn.c_str(), n_i_calc[k], Xmol_i_calc[k]); } - } -#ifdef DEBUG_HKM_EPEQUIL - if (debug_prnt_lvl > 0) { + printf("%15s: %10.5g\n", "Total Molar Sum", n_t_calc); printf("(iter %d) element moles bal: Goal Calculated\n", iter); for (m = 0; m < m_mm; m++) { string nnn = eNames[m]; diff --git a/Cantera/src/ChemEquil.h b/Cantera/src/ChemEquil.h index 06a1a90f4..a5443078e 100755 --- a/Cantera/src/ChemEquil.h +++ b/Cantera/src/ChemEquil.h @@ -147,6 +147,10 @@ namespace Cantera { 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); + int m_mm; int m_kk; int m_skip; @@ -192,6 +196,12 @@ namespace Cantera { 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; /* diff --git a/Cantera/src/Makefile.in b/Cantera/src/Makefile.in index 72dc85bd5..678aa0501 100755 --- a/Cantera/src/Makefile.in +++ b/Cantera/src/Makefile.in @@ -78,7 +78,7 @@ HETEROKIN = $(THERMO) $(HETEROKIN_OBJ) CK = $(KINETICS) # chemical equilibrium -EQUIL_OBJ = ChemEquil.o MultiPhaseEquil.o sort.o MultiPhase.o equilibrate.o +EQUIL_OBJ = ChemEquil.o MultiPhaseEquil.o sort.o MultiPhase.o equilibrate.o BasisOptimize.o EQUIL_H = ChemEquil.h MultiPhaseEquil.h MultiPhase.h Crystal.h sort.h equil.h EQUIL = $(THERMO) $(EQUIL_OBJ) diff --git a/Cantera/src/MultiPhase.cpp b/Cantera/src/MultiPhase.cpp index 9790db104..36ed50690 100644 --- a/Cantera/src/MultiPhase.cpp +++ b/Cantera/src/MultiPhase.cpp @@ -367,6 +367,25 @@ namespace Cantera { 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 @@ -392,7 +411,7 @@ namespace Cantera { else { m_moleFractions[loc] = 1.0; } - loc += p->nSpecies(); + loc += nsp; } } diff --git a/Cantera/src/MultiPhase.h b/Cantera/src/MultiPhase.h index bcf1f40ef..0bbd2ed01 100644 --- a/Cantera/src/MultiPhase.h +++ b/Cantera/src/MultiPhase.h @@ -218,6 +218,7 @@ namespace Cantera { void setMolesByName(const std::string& x); + void getMoles(doublereal * molNum) const; void setMoles(doublereal* n); /// Return true if the phase \a p has valid thermo data for @@ -296,6 +297,12 @@ namespace Cantera { } return s; } + + +int BasisOptimize( int *usedZeroedSpecies, bool doFormRxn, + MultiPhase *mphase, vector_int & orderVectorSpecies, + vector_int & orderVectorElements, + vector_fp & formRxnMatrix); } #endif