diff --git a/Cantera/src/ChemEquil.h b/Cantera/src/ChemEquil.h index 4caa68a93..bcff25c07 100755 --- a/Cantera/src/ChemEquil.h +++ b/Cantera/src/ChemEquil.h @@ -159,6 +159,7 @@ namespace Cantera { if (solver > 0) { MultiPhase mix; mix.addPhase(&s, 1.0); + mix.init(); mix.setTemperature(s.temperature()); mix.setPressure(s.pressure()); equilibrate(mix, XY, rtol, maxsteps, loglevel); diff --git a/Cantera/src/MultiPhase.cpp b/Cantera/src/MultiPhase.cpp index eaab1e02d..bd747ffb9 100644 --- a/Cantera/src/MultiPhase.cpp +++ b/Cantera/src/MultiPhase.cpp @@ -1,11 +1,3 @@ -/** - * @file MultiPhase.cpp - * - * $Author$ - * $Date$ - * $Revision$ - */ - #include "MultiPhase.h" #include "MultiPhaseEquil.h" @@ -15,16 +7,24 @@ 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) { + m_Tmin(1.0), m_Tmax(100000.0) { } void MultiPhase:: - addPhase(phase_t* p, doublereal moles) { + 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."); @@ -56,33 +56,47 @@ namespace Cantera { // 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[ename] == 0) { + 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(); } - //cout << "min, max = " << m_Tmin << " " << m_Tmax << endl; + + // 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; - //cout << p->name() << " " << t << " " << m_Tmax << endl; } - } - /// Process phases and build atomic composition array. After - /// init() has been called, no more phases may be added. + /// 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; @@ -149,7 +163,6 @@ namespace Cantera { /// Moles of species \c k. doublereal MultiPhase::speciesMoles(index_t k) { - if (!m_init) init(); index_t ip = m_spphase[k]; return m_moles[ip]*m_moleFractions[k]; } @@ -159,7 +172,6 @@ namespace Cantera { doublereal MultiPhase::elementMoles(index_t m) { doublereal sum = 0.0, phasesum; index_t i, k = 0, ik, nsp; - if (!m_init) init(); for (i = 0; i < m_np; i++) { phasesum = 0.0; nsp = m_phase[i]->nSpecies(); @@ -182,7 +194,10 @@ namespace Cantera { return sum; } - /// Charge of one phase + /// 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) { doublereal phasesum = 0.0; int ik, k, nsp = m_phase[p]->nSpecies(); @@ -194,6 +209,7 @@ namespace Cantera { } + /// Get the chemical potentials of all species in all phases. void MultiPhase::getChemPotentials(doublereal* mu) { index_t i, loc = 0; updatePhases(); @@ -203,10 +219,39 @@ namespace Cantera { } } + /// 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) { index_t i, loc = 0; - updatePhases(); + + updatePhases(); + // iterate over the phases for (i = 0; i < m_np; i++) { if (tempOK(i) || m_phase[i]->nSpecies() > 1) { if (!standard) @@ -220,7 +265,7 @@ namespace Cantera { } } - + /// True if species \a k belongs to a solution phase. bool MultiPhase::solutionSpecies(index_t k) { if (m_phase[m_spphase[k]]->nSpecies() > 1) return true; @@ -228,6 +273,7 @@ namespace Cantera { return false; } + /// The Gibbs free energy of the mixture (J). doublereal MultiPhase::gibbs() const { index_t i; doublereal sum = 0.0; @@ -237,6 +283,7 @@ namespace Cantera { return sum; } + /// The enthalpy of the mixture (J). doublereal MultiPhase::enthalpy() const { index_t i; doublereal sum = 0.0; @@ -246,6 +293,7 @@ namespace Cantera { return sum; } + /// The entropy of the mixture (J/K). doublereal MultiPhase::entropy() const { index_t i; doublereal sum = 0.0; @@ -255,6 +303,9 @@ namespace Cantera { 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; @@ -264,45 +315,52 @@ namespace Cantera { return sum; } - void MultiPhase::updateMoleFractions() { - if (!m_init) init(); - // save the current mole fractions for each phase - index_t ip, loc = 0; - for (ip = 0; ip < m_np; ip++) { - phase_t* p = m_phase[ip]; - p->getMoleFractions(m_moleFractions.begin() + loc); - loc += p->nSpecies(); - } - } + + /// 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) { - if (!m_init) init(); int kk = nSpecies(); doublereal x; - vector_fp mf(kk, 0.0); + vector_fp moles(kk, 0.0); for (int k = 0; k < kk; k++) { x = xMap[speciesName(k)]; - if (x > 0.0) mf[k] = x; + if (x > 0.0) moles[k] = x; } - setMoles(mf.begin()); + setMoles(moles.begin()); } + + /// Set the species moles using a string. Unspecified species are + /// set to zero. void MultiPhase::setMolesByName(const string& x) { compositionMap xx; - if (!m_init) init(); + + // 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); } + /// 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; @@ -328,6 +386,7 @@ namespace Cantera { } } + /// The total mixture volume [m^3]. doublereal MultiPhase::volume() { int i; doublereal sum = 0; @@ -337,30 +396,6 @@ namespace Cantera { return sum; } - /** - * updatePhases() const: - * - * This routine take the state information storred in - * this current object, and sets the state in the underlying - * ThermoPhase objects. - * Currently, this state information is the Temperature - * Pressure, and the mole fraction vector. Temperature and - * pressure are assumed uniform for all phases. - */ - 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(); - const doublereal* x = m_moleFractions.begin() + 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; - } - } - doublereal MultiPhase::equilibrate(int XY, doublereal err, int maxsteps, int maxiter, int loglevel) { doublereal error; @@ -370,9 +405,11 @@ namespace Cantera { int n; bool start, once; doublereal ferr, hnow, herr = 1.0; + doublereal snow, serr = 1.0, s0; doublereal Tlow = -1.0, Thigh = -1.0; - doublereal hlow = 0.0, hhigh = 0.0, tnew; - doublereal dta, dtmax; + doublereal Hlow = Undef, Hhigh = Undef, tnew; + doublereal dta, dtmax, cpb; + MultiPhaseEquil* e = 0; if (!m_init) init(); if (loglevel > 0) { @@ -383,59 +420,177 @@ namespace Cantera { addLogEntry("problem type","fixed T,P"); } // create an equilibrium manager - MultiPhaseEquil e(this); - error = e.equilibrate(XY, err, maxsteps, loglevel-1); - // if (loglevel > 0) e.printInfo(); + e = new MultiPhaseEquil(this); + try { + error = e->equilibrate(XY, err, maxsteps, loglevel-1); + } + catch (CanteraError err) { + if (loglevel > 0) { + endLogGroup(); + //write_logfile("equil_err.html"); + } + delete e; + e = 0; + throw err; + } goto done; } + else if (XY == HP) { h0 = enthalpy(); - start = true; - Tlow = m_Tmin; // lower bound on T - Thigh = m_Tmax; // upper bound on T - hlow = 0.0; - hhigh = 0.0; - once = true; + Tlow = 0.5*m_Tmin; // lower bound on T + Thigh = 2.0*m_Tmax; // upper bound on T if (loglevel > 0) { addLogEntry("problem type","fixed H,P"); addLogEntry("H target",fp2str(h0)); - addLogEntry("min T",fp2str(Tlow)); - addLogEntry("max T",fp2str(Thigh)); + //addLogEntry("min T",fp2str(Tlow)); + //addLogEntry("max T",fp2str(Thigh)); } for (n = 0; n < maxiter; n++) { - MultiPhaseEquil e(this, strt); - ferr = 0.1; - if (fabs(dt) < 1.0) ferr = err; - start = false; + + // if 'strt' is false, the current composition will be used as + // the starting estimate; otherwise it will be estimated + e = new MultiPhaseEquil(this, strt); + // start with a loose error tolerance, but tighten it as we get + // close to the final temperature if (loglevel > 0) { beginLogGroup("iteration "+int2str(n)); } try { - error = e.equilibrate(TP, err, maxsteps, loglevel-1); - + error = e->equilibrate(TP, err, maxsteps, loglevel-1); 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; + 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 { - if (m_temp < Thigh) Thigh = m_temp; + tnew = sqrt(Tlow*Thigh); + dt = tnew - m_temp; + //cpb = cp(); } + herr = fabs((h0 - hnow)/h0); if (loglevel > 0) { 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(); } - dt = (h0 - hnow)/cp(); - dtmax = 0.5*(Thigh - Tlow); + + + if (herr < err) { // || dta < 1.0e-4) { + if (loglevel > 0) { + 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) { + if (loglevel > 0) + addLogEntry("no convergence", + "try estimating starting composition"); + strt = true; + } + else { + tnew = 0.5*(Tlow + Thigh); + setTemperature(tnew); + if (loglevel > 0) + addLogEntry("no convergence", + "trying T = "+fp2str(m_temp)); + } + endLogGroup(); + } + } + delete e; + e = 0; + if (loglevel > 0) { + addLogEntry("reached max number of T iterations",int2str(maxiter)); + endLogGroup(); + } + throw CanteraError("MultiPhase::equilibrate", + "No convergence for T"); + } + else if (XY == SP) { + writelog("SP\n"); + s0 = entropy(); + start = true; + Tlow = 1.0; // m_Tmin; // lower bound on T + Thigh = 1.0e6; // m_Tmax; // upper bound on T + if (loglevel > 0) { + 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++) { + e = new MultiPhaseEquil(this, strt); + ferr = 0.1; + if (fabs(dt) < 1.0) ferr = err; + //start = false; + if (loglevel > 0) { + beginLogGroup("iteration "+int2str(n)); + } + try { + error = e->equilibrate(TP, err, maxsteps, loglevel-1); + 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); + if (loglevel > 0) { + 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) { if (loglevel > 0) { addLogEntry("T iterations",int2str(n)); addLogEntry("Final T",fp2str(temperature())); - addLogEntry("H rel error",fp2str(herr)); + addLogEntry("S rel error",fp2str(serr)); } goto done; } @@ -446,7 +601,8 @@ namespace Cantera { // the current composition as the starting estimate if (dta < 100.0) strt = false; } - catch (CanteraError e) { + + catch (CanteraError err) { if (!strt) { if (loglevel > 0) addLogEntry("no convergence", @@ -461,58 +617,89 @@ namespace Cantera { "trying T = "+fp2str(m_temp)); } + endLogGroup(); } } + delete e; + e = 0; + if (loglevel > 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++) { - MultiPhaseEquil e(this, start); - ferr = 0.1; - start = false; - if (fabs(dt) < 1.0) ferr = err; - if (loglevel > 1) { - beginLogGroup("iteration "+int2str(n)); - } - 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); - } - } + + +// 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) { if (loglevel > 0) { addLogEntry("problem type","fixed T, V"); } - //doublereal dt = 1.0e3; + doublereal dt = 1.0e3; doublereal v0 = volume(); doublereal dVdP; int n; bool start = true; - doublereal error, vnow, pnow, verr; + doublereal error, ferr, vnow, pnow, verr, tnew; for (n = 0; n < maxiter; n++) { pnow = pressure(); MultiPhaseEquil e(this, start); @@ -549,10 +736,50 @@ namespace Cantera { } return -1.0; done: + delete e; + e = 0; if (loglevel > 0) { endLogGroup(); } return err; } + + + + //------------------------------------------------------------- + // + // 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(m_moleFractions.begin() + loc); + loc += p->nSpecies(); + } + } + + + /// 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 = m_moleFractions.begin() + 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/MultiPhase.h b/Cantera/src/MultiPhase.h index dc82ee712..e2caf96b9 100644 --- a/Cantera/src/MultiPhase.h +++ b/Cantera/src/MultiPhase.h @@ -7,8 +7,6 @@ namespace Cantera { - class MultiPhaseEquil; - /// 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. @@ -22,9 +20,11 @@ namespace Cantera { public: + // some typedefs for convenience typedef size_t index_t; typedef ThermoPhase phase_t; typedef DenseMatrix array_t; + typedef vector phase_list; /// Constructor. The constructor takes no arguments, since /// phases are added using method addPhase. @@ -35,6 +35,8 @@ namespace Cantera { /// phase objects. virtual ~MultiPhase() {} + void addPhases(phase_list& phases, const vector_fp& phaseMoles); + /// 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 @@ -215,6 +217,8 @@ namespace Cantera { protected: + // These methods are meant for internal use. + /// update the locally-stored composition to match the current /// compositions of the phase objects. void updateMoleFractions(); @@ -242,7 +246,6 @@ namespace Cantera { bool m_init; int m_eloc; mutable vector m_temp_OK; - MultiPhaseEquil* m_equil; doublereal m_Tmin, m_Tmax; };