From 88d45066c7904027f7bdb266ad4b406b118dd87c Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Mon, 27 Jun 2005 14:51:05 +0000 Subject: [PATCH] cleanup --- Cantera/src/MultiPhase.cpp | 50 +++---- Cantera/src/MultiPhase.h | 98 ++++++++++---- Cantera/src/MultiPhaseEquil.cpp | 228 +++++++++++--------------------- Cantera/src/MultiPhaseEquil.h | 10 +- 4 files changed, 174 insertions(+), 212 deletions(-) diff --git a/Cantera/src/MultiPhase.cpp b/Cantera/src/MultiPhase.cpp index 5ce60ca68..0e093ce8f 100644 --- a/Cantera/src/MultiPhase.cpp +++ b/Cantera/src/MultiPhase.cpp @@ -8,13 +8,12 @@ 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_equil(0), m_Tmin(1.0), m_Tmax(100000.0) { } - - void MultiPhase:: addPhase(phase_t* p, doublereal moles) { @@ -187,8 +186,6 @@ namespace Cantera { } - /// Chemical potentials. Write into array \c mu the chemical - /// potentials of all species [J/kmol]. void MultiPhase::getChemPotentials(doublereal* mu) { index_t i, loc = 0; updatePhases(); @@ -198,8 +195,6 @@ namespace Cantera { } } - /// Chemical potentials. Write into array \c mu the chemical - /// potentials of all species [J/kmol]. void MultiPhase::getValidChemPotentials(doublereal not_mu, doublereal* mu, bool standard) { index_t i, loc = 0; @@ -218,17 +213,6 @@ namespace Cantera { } - /// Chemical potentials. Write into array \c mu the chemical - /// potentials of all species [J/kmol]. - void MultiPhase::getStandardChemPotentials(doublereal* mu) { - index_t i, loc = 0; - updatePhases(); - for (i = 0; i < m_np; i++) { - m_phase[i]->getStandardChemPotentials(mu + loc); - loc += m_phase[i]->nSpecies(); - } - } - bool MultiPhase::solutionSpecies(index_t k) { if (m_phase[m_spphase[k]]->nSpecies() > 1) return true; @@ -236,7 +220,7 @@ namespace Cantera { return false; } - doublereal MultiPhase::gibbs() { + doublereal MultiPhase::gibbs() const { index_t i; doublereal sum = 0.0; updatePhases(); @@ -245,7 +229,7 @@ namespace Cantera { return sum; } - doublereal MultiPhase::enthalpy() { + doublereal MultiPhase::enthalpy() const { index_t i; doublereal sum = 0.0; updatePhases(); @@ -254,7 +238,7 @@ namespace Cantera { return sum; } - doublereal MultiPhase::entropy() { + doublereal MultiPhase::entropy() const { index_t i; doublereal sum = 0.0; updatePhases(); @@ -263,7 +247,7 @@ namespace Cantera { return sum; } - doublereal MultiPhase::cp() { + doublereal MultiPhase::cp() const { index_t i; doublereal sum = 0.0; updatePhases(); @@ -345,12 +329,12 @@ namespace Cantera { return sum; } - void MultiPhase::updatePhases() { + void MultiPhase::updatePhases() const { if (!m_init) init(); index_t p, nsp, loc = 0; for (p = 0; p < m_np; p++) { nsp = m_phase[p]->nSpecies(); - doublereal* x = m_moleFractions.begin() + loc; + const doublereal* x = m_moleFractions.begin() + loc; loc += nsp; m_phase[p]->setState_TPX(m_temp, m_press, x); m_temp_OK[p] = true; @@ -383,11 +367,10 @@ namespace Cantera { // create an equilibrium manager MultiPhaseEquil e(this); error = e.equilibrate(XY, err, maxsteps, loglevel-1); - if (loglevel > 0) e.printInfo(); + // if (loglevel > 0) e.printInfo(); goto done; } else if (XY == HP) { - dt = 1.0e2; h0 = enthalpy(); start = true; Tlow = m_Tmin; // lower bound on T @@ -401,9 +384,10 @@ namespace Cantera { addLogEntry("min T",fp2str(Tlow)); addLogEntry("max T",fp2str(Thigh)); } - ferr = 0.1; for (n = 0; n < maxiter; n++) { MultiPhaseEquil e(this, strt); + ferr = 0.1; + if (fabs(dt) < 1.0) ferr = err; start = false; if (loglevel > 0) { beginLogGroup("iteration "+int2str(n)); @@ -413,14 +397,10 @@ namespace Cantera { hnow = enthalpy(); if (hnow < h0) { - if (m_temp > Tlow) { - Tlow = m_temp; - } + if (m_temp > Tlow) Tlow = m_temp; } else { - if (m_temp < Thigh) { - Thigh = m_temp; - } + if (m_temp < Thigh) Thigh = m_temp; } herr = fabs((h0 - hnow)/h0); if (loglevel > 0) { @@ -443,12 +423,16 @@ namespace Cantera { } 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 e) { if (!strt) { if (loglevel > 0) - addLogEntry("no convergence","setting strt to True"); + addLogEntry("no convergence", + "setting strt to True"); strt = true; } else { diff --git a/Cantera/src/MultiPhase.h b/Cantera/src/MultiPhase.h index a44de170d..dc82ee712 100644 --- a/Cantera/src/MultiPhase.h +++ b/Cantera/src/MultiPhase.h @@ -22,6 +22,10 @@ namespace Cantera { public: + typedef size_t index_t; + typedef ThermoPhase phase_t; + typedef DenseMatrix array_t; + /// Constructor. The constructor takes no arguments, since /// phases are added using method addPhase. MultiPhase(); @@ -31,22 +35,27 @@ namespace Cantera { /// phase objects. virtual ~MultiPhase() {} - typedef size_t index_t; - typedef ThermoPhase phase_t; - typedef DenseMatrix array_t; - - /// Add a phase to the mixture. /// @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() { return int(m_nel); } + + /// Name of element \a m. string elementName(int m) { return m_enames[m]; } + + /// Index of element with name \a name. int elementIndex(string name) { return m_enamemap[name] - 1;} + /// Number of species, summed over all phases. int nSpecies() { return int(m_nsp); } + + /// Name of species with index \a k. string speciesName(int k) { return m_snames[k]; } + + /// Number of atoms of element \a m in species \a k. doublereal nAtoms(int k, int m) { if (!m_init) init(); return m_atoms(m,k); @@ -87,36 +96,56 @@ namespace Cantera { 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(); + + /// 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(); + /// Total charge (Coulombs). doublereal charge(); + + /// Charge (Coulombs) of phase with index \a p. doublereal phaseCharge(index_t p); - /// Total moles of element m, summed over all - /// phases + /// Total moles of element \a m, summed over all phases. doublereal elementMoles(index_t m); - /// Chemical potentials. Write into array \c mu the chemical - /// potentials of all species [J/kmol]. + /// 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]. void getChemPotentials(doublereal* mu); - /// Valid chemical potentials. Write into array \c mu the + /// 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 \c not_mu. + /// 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. void getValidChemPotentials(doublereal not_mu, doublereal* mu, bool standard = false); - /// Chemical potentials. Write into array \c mu the chemical - /// potentials of all species [J/kmol]. - void getStandardChemPotentials(doublereal* mu); - /// Temperature [K]. - doublereal temperature() { - return m_temp; - } + doublereal temperature() { 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 = 0); @@ -127,26 +156,39 @@ namespace Cantera { updatePhases(); } + /// Pressure [Pa]. doublereal pressure() { return m_press; } + /// Volume [m^3]. doublereal volume(); + /// Set the pressure [Pa]. void setPressure(doublereal P) { m_press = P; updatePhases(); } - doublereal enthalpy(); - doublereal entropy(); - doublereal gibbs(); - doublereal cp(); + /// 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() { return m_np; } + /// Return true is species \a k is a species in a + /// multicomponent solution phase. bool solutionSpecies(index_t k); index_t speciesPhaseIndex(index_t k) { @@ -157,8 +199,6 @@ namespace Cantera { return m_moleFractions[k]; } - void updateMoleFractions(); - void setPhaseMoleFractions(index_t n, doublereal* x); void setMolesByName(compositionMap& xMap); @@ -167,16 +207,22 @@ namespace Cantera { void setMoles(doublereal* n); + /// Return true if the phase \a p has valid thermo data for + /// the current temperature. bool tempOK(index_t p) { return m_temp_OK[p]; } protected: + /// update the locally-stored composition to match the current + /// compositions of the phase objects. + void updateMoleFractions(); + /// 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(); + void updatePhases() const; vector_fp m_moles; vector m_phase; @@ -195,7 +241,7 @@ namespace Cantera { index_t m_nsp; bool m_init; int m_eloc; - vector m_temp_OK; + mutable vector m_temp_OK; MultiPhaseEquil* m_equil; doublereal m_Tmin, m_Tmax; }; diff --git a/Cantera/src/MultiPhaseEquil.cpp b/Cantera/src/MultiPhaseEquil.cpp index ba347b2be..1d4878b4e 100644 --- a/Cantera/src/MultiPhaseEquil.cpp +++ b/Cantera/src/MultiPhaseEquil.cpp @@ -1,7 +1,6 @@ #include "MultiPhaseEquil.h" #include "MultiPhase.h" #include "sort.h" -#include "recipes.h" #include "global.h" #include @@ -19,8 +18,6 @@ using namespace std; #endif #endif -#undef DEBUG_MULTIPHASE_EQUIL - namespace Cantera { const doublereal TINY = 1.0e-20; @@ -140,8 +137,9 @@ namespace Cantera { m_N.resize(m_nsp, m_nsp - m_nel); m_order.resize(m_nsp, 0); - if (start) + if (start) { setInitialMoles(); + } computeN(); vector_fp dxi(m_nsp - m_nel, 1.0e-20); @@ -157,9 +155,10 @@ namespace Cantera { m_dsoln.push_back(0); } m_force = false; - setMoles(); + updateMixMoles(); } + doublereal MultiPhaseEquil::equilibrate(int XY, doublereal err, int maxsteps, int loglevel) { int i; @@ -175,10 +174,6 @@ namespace Cantera { endLogGroup(); } if (loglevel > 2) printInfo(); - //if (error() == 0.0) { - // write_logfile("equil_err.html"); - // Cantera::error("stopping"); - //} if (error() < err) break; } if (i >= maxsteps) { @@ -202,8 +197,7 @@ namespace Cantera { return error(); } - void MultiPhaseEquil::setMoles() { - //vector_fp n(m_nsp_mix, 0.0); + void MultiPhaseEquil::updateMixMoles() { fill(m_work3.begin(), m_work3.end(), 0.0); index_t k; for (k = 0; k < m_nsp; k++) { @@ -212,10 +206,10 @@ namespace Cantera { m_mix->setMoles(m_work3.begin()); } - /// Clean up the composition by setting species with negative mole - /// numbers to zero. The solution algorithm can leave some species - /// in stoichiometric condensed phases with very small negative - /// mole numbers. This method simply sets these to zero. + /// 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; @@ -226,73 +220,69 @@ namespace Cantera { } - /** - * Estimate the initial mole fractions. Uses the Simplex method - * to estimate the initial number of moles of each species. The - * linear Gibbs minimization problem is solved, neglecting the - * free energy of mixing terms. This procedure produces a good - * estimate of the low-temperature equilibrium composition. - * - * @param s phase object - * @param elementMoles vector of elemental moles - */ + /// 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 m, n; - doublereal lp = log(m_press/OneAtm); - - DenseMatrix aa(m_nel+2, m_nsp+1, 0.0); - - // first column contains fixed element moles - for (m = 0; m < m_nel; m++) { - aa(m+1,0) = m_mix->elementMoles(m_element[m]); - } + index_t m, n, ik, j; - // get the array of non-dimensional Gibbs functions for the pure - // species - //m_mix->getStandardChemPotentials(m_mu.begin()); double not_mu = 1.0e12; m_mix->getValidChemPotentials(not_mu, m_mu.begin(), true); - - int kpp = 0; - index_t k, q; - doublereal rt = GasConstant * m_temp; - for (k = 0; k < m_nsp; k++) { - kpp++; - aa(0, kpp) = -m_mu[m_species[k]]/rt; - aa(0, kpp) -= m_dsoln[k]*lp; // ideal gas - for (q = 0; q < m_nel; q++) - aa(q+1, kpp) = -m_mix->nAtoms(m_species[k], m_element[q]); - } + doublereal dg_rt; - integer mp = m_nel+2; // parameters for SIMPLX - integer np = m_nsp+1; - integer m1 = 0; - integer m2 = 0; - integer m3 = m_nel; - integer icase=0; - integer nel = m_nel; - integer nsp = m_nsp; - vector_int iposv(m_nel); - vector_int izrov(m_nsp); - - // solve the linear programming problem + int idir; + double nu; + double delta_xi, dxi_min = 1.0e10; + bool redo = true; + int iter = 0; + while (redo) { - simplx_(&aa(0,0), &nel, &nsp, &mp, &np, &m1, &m2, &m3, - &icase, izrov.begin(), iposv.begin()); - - fill(m_moles.begin(), m_moles.end(), 0.0); - for (n = 0; n < m_nel; n++) { - int ksp = 0; - int ip = iposv[n] - 1; - for (int k = 0; k < int(m_nsp); k++) { - if (ip == ksp) { - m_moles[k] = aa(n+1, 0); + // choose a set of components based on the current + // composition + computeN(); + + 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 (delta_xi < SmallNumber && ik < m_nel) 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; } - ksp++; } + // set the moles of the phase objects to match + updateMixMoles(); } - setMoles(); - return icase; + return 0; } @@ -307,30 +297,20 @@ namespace Cantera { /// The constituent species are taken to be the first M species /// in array 'species' that have linearly-independent compositions. /// - /// Arguments: + /// @param order On entry, vector \a order should contain species + /// index numbers in the order of decreasing desirability as a + /// constituent. For example, if it is desired to choose the + /// constituents from among the major species, this array might + /// list species index numbers in decreasing order of mole + /// fraction. If array 'species' does not have length = + /// nSpecies(), then the species will be considered as candidates + /// to be constituents in declaration order, beginning with the + /// first phase added. /// - /// On entry, vector species shold contain species index numbers - /// in the order of decreasing desirability as a constituent. For - /// example, if it is desired to choose the constituents from - /// among the major species, this array might list species index - /// numbers in decreasing order of mole fraction. If array - /// 'species' does not have length = nSpecies(), then the species - /// will be considered as candidates to be constituents in - /// declaration order, beginning with the first phase added. - /// - /// On return, the first M entries of array 'species' contain the index - /// numbers of the constituent species. - /// - /// Matrix nu is an output array that contains the stoichiometric - /// coefficents for a set of K - M formation reactions for the - /// non-constituent species, such that nu(k,i) is the net - /// stoichiometric coefficent of species k in reaction i. Matrix - /// nu will be resized to (K, K-M) and its initial values, if - /// any, will be erased. - void MultiPhaseEquil::getComponents(const vector_int& order) { index_t m, k, j; int n; + // if the input species array has the wrong size, ignore it // and consider the species for constituents in declarationi order. if (order.size() != m_nsp) { @@ -346,13 +326,6 @@ namespace Cantera { index_t nColumns = m_nsp; doublereal fctr; -#ifdef DEBUG_MULTIPHASE_EQUIL - cout << "most abundant:" << endl; - for (m = 0; m < nRows; m++) { - cout << m_mix->speciesName(m_species[m_order[m]]) << " " << m_moles[m_order[m]] << endl; - } -#endif - // set up the atomic composition matrix for (m = 0; m < nRows; m++) { for (k = 0; k < nColumns; k++) { @@ -380,17 +353,10 @@ namespace Cantera { m_A(n, kmax) = tmp; } // exchange the species labels on the columns -#ifdef DEBUG_MULTIPHASE_EQUIL - cout << "in row " << m << ", pivot is zero" << endl; - cout << "exchanging " << m_mix->speciesName(m_species[m_order[m]]) << " for " << m_mix->speciesName(m_species[m_order[kmax]]) << endl; -#endif itmp = m_order[m]; m_order[m] = m_order[kmax]; m_order[kmax] = itmp; - // throw an exception if the entire row is zero - // if (k >= m_nsp) - // throw CanteraError("getComponents","all zeros!"); } // scale row m so that the diagonal element is unity @@ -423,23 +389,6 @@ namespace Cantera { } } -#ifdef DEBUG_MULTIPHASE_EQUIL - // check - bool ok = true; - for (m = 0; m < nRows; m++) { - cout << m_mix->speciesName(m_species[m_order[m]]) << " " << m_moles[m_order[m]] << endl; - if (m_A(m,m) != 1.0) ok = false; - for (n = 0; n < nRows; n++) { - if (n != m && fabs(m_A(m,n)) > TINY) - ok = false; - } - } - if (!ok) { - cout << m_A << endl; - throw CanteraError("getComponents","error in A matrix"); - } -#endif - // create stoichometric coefficient matrix. for (n = 0; n < int(m_nsp); n++) { if (n < int(m_nel)) @@ -462,17 +411,6 @@ namespace Cantera { } } - - /// Re-arrange a vector of species properties in sequential form - /// into sorted (components first) form. - void MultiPhaseEquil::sort(vector_fp& x) { - copy(x.begin(), x.end(), m_work2.begin()); - index_t k; - for (k = 0; k < m_nsp; k++) { - x[k] = m_work2[m_order[k]]; - } - } - /// Re-arrange a vector of species properties in sorted form /// (components first) into unsorted, sequential form. void MultiPhaseEquil::unsort(vector_fp& x) { @@ -549,12 +487,14 @@ namespace Cantera { m_moles[k] += omega * deltaN[k]; } else { - m_moles[k] = fabs(m_moles[k])*fminn(10.0, exp(-m_deltaG_RT[ik - m_nel])); + m_moles[k] = fabs(m_moles[k])*fminn(10.0, + exp(-m_deltaG_RT[ik - m_nel])); } } - setMoles(); + updateMixMoles(); } + /// Take one step in composition, given the gradient of G at the /// starting point, and a vector of reaction steps dxi. doublereal MultiPhaseEquil:: @@ -567,14 +507,7 @@ namespace Cantera { doublereal grad0 = computeReactionSteps(m_dxi); // compute the mole fraction changes. - //multiply(m_N, dxi.begin(), m_work.begin()); - for (ik = 0; ik < m_nsp; ik++) { - m_work[ik] = 0.0; - k = m_order[ik]; - for (j = 0; j < m_nsp - m_nel; j++) { - m_work[ik] += m_N(ik, j) * m_dxi[j]; - } - } + multiply(m_N, m_dxi.begin(), m_work.begin()); // change to sequential form unsort(m_work); @@ -643,7 +576,6 @@ namespace Cantera { for (k = 0; k < m_nsp; k++) { grad1 += m_work[k] * m_mu[m_species[k]]; } - // doublereal grad1 = dot(m_work.begin(), m_work.end(), m_work2.begin()); omega = omegamax; if (grad1 > 0.0) { @@ -739,10 +671,6 @@ namespace Cantera { fctr = 1.0; else fctr = 1.0/(term1 + csum + sum); - //if (fctr < -999.0 || fctr > 999.0) { - // cout << "fctr, term1, csum, sum = " << fctr << " " << term1 << " " << csum << " " << sum << endl; - // cout << reactionString(j) << endl; - //} } dxi[j] = -fctr*dg_rt; index_t m; diff --git a/Cantera/src/MultiPhaseEquil.h b/Cantera/src/MultiPhaseEquil.h index e67348154..9fa86726d 100644 --- a/Cantera/src/MultiPhaseEquil.h +++ b/Cantera/src/MultiPhaseEquil.h @@ -47,13 +47,14 @@ namespace Cantera { void getComponents(const vector_int& order); int setInitialMoles(); + int setInitialMoles2(); void computeN(); doublereal stepComposition(int loglevel); - void sort(vector_fp& x); + //void sort(vector_fp& x); void unsort(vector_fp& x); void step(doublereal omega, vector_fp& deltaN); doublereal computeReactionSteps(vector_fp& dxi); - void setMoles(); + void updateMixMoles(); void finish(); // moles of the species with sorted index ns @@ -61,7 +62,10 @@ namespace Cantera { 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]]]; } + 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;