MultiPhase.cpp

Go to the documentation of this file.
00001 /**
00002  * @file MultiPhase.cpp
00003  * Definitions for the \link Cantera::MultiPhase MultiPhase\endlink 
00004  * object that is used to set up multiphase equilibrium problems (see \ref equilfunctions).
00005  */
00006 /*
00007  *
00008  *  $Date: 2010-01-03 19:46:26 -0500 (Sun, 03 Jan 2010) $
00009  *  $Revision: 368 $
00010  */
00011 
00012 #include "MultiPhase.h"
00013 #include "MultiPhaseEquil.h"
00014 
00015 #include "ThermoPhase.h"
00016 #include "DenseMatrix.h"
00017 #include "stringUtils.h"
00018 #include "global.h"
00019 
00020 using namespace std;
00021 
00022 namespace Cantera {
00023 
00024   /// Constructor.
00025   MultiPhase::MultiPhase() :
00026     m_np(0),
00027     m_temp(0.0),
00028     m_press(0.0), 
00029     m_nel(0),
00030     m_nsp(0),
00031     m_init(false),
00032     m_eloc(-1), 
00033     m_Tmin(1.0),
00034     m_Tmax(100000.0)
00035   {
00036   }
00037 
00038   void MultiPhase::
00039   addPhases(MultiPhase& mix) {
00040     index_t n;
00041     for (n = 0; n < mix.m_np; n++) {
00042       addPhase(mix.m_phase[n], mix.m_moles[n]);
00043     }
00044   }
00045 
00046   void MultiPhase::
00047   addPhases(phase_list& phases, const vector_fp& phaseMoles) {
00048     index_t np = phases.size();
00049     index_t n;
00050     for (n = 0; n < np; n++) {
00051       addPhase(phases[n], phaseMoles[n]);
00052     }
00053     init();
00054   }
00055 
00056   void MultiPhase::
00057   addPhase(phase_t* p, doublereal moles) {
00058     if (m_init) {
00059       throw CanteraError("addPhase",
00060                          "phases cannot be added after init() has been called.");
00061     }
00062 
00063     // save the pointer to the phase object
00064     m_phase.push_back(p);
00065 
00066     // store its number of moles
00067     m_moles.push_back(moles);
00068     m_temp_OK.push_back(true);
00069 
00070     // update the number of phases and the total number of
00071     // species
00072     m_np = m_phase.size();
00073     m_nsp += p->nSpecies();
00074 
00075     // determine if this phase has new elements
00076     // for each new element, add an entry in the map
00077     // from names to index number + 1:
00078 
00079     string ename;
00080     // iterate over the elements in this phase
00081     index_t m, nel = p->nElements();
00082     for (m = 0; m < nel; m++) {
00083       ename = p->elementName(m);
00084 
00085       // if no entry is found for this element name, then
00086       // it is a new element. In this case, add the name
00087       // to the list of names, increment the element count, 
00088       // and add an entry to the name->(index+1) map.
00089       if (m_enamemap.find(ename) == m_enamemap.end()) {
00090         m_enamemap[ename] = m_nel + 1;
00091         m_enames.push_back(ename);
00092         m_atomicNumber.push_back(p->atomicNumber(m));
00093 
00094         // Element 'E' (or 'e') is special. Note its location.
00095         if (ename == "E" || ename == "e") m_eloc = m_nel;
00096 
00097         m_nel++;
00098       }
00099     }
00100 
00101     // If the mixture temperature hasn't been set, then set the
00102     // temperature and pressure to the values for the phase being
00103     // added.
00104     if (m_temp == 0.0 && p->temperature() > 0.0) {
00105       m_temp = p->temperature();
00106       m_press = p->pressure();
00107     }
00108 
00109     // If this is a solution phase, update the minimum and maximum
00110     // mixture temperatures. Stoichiometric phases are excluded,
00111     // since a mixture may define multiple stoichiometric phases,
00112     // each of which has thermo data valid only over a limited
00113     // range. For example, a mixture might be defined to contain a
00114     // phase representing water ice and one representing liquid
00115     // water, only one of which should be present if the mixture
00116     // represents an equilibrium state. 
00117     if (p->nSpecies() > 1) {
00118       double t = p->minTemp();
00119       if (t > m_Tmin) m_Tmin = t;
00120       t = p->maxTemp();
00121       if (t < m_Tmax) m_Tmax = t;
00122     }
00123   }
00124 
00125 
00126   // Process phases and build atomic composition array. This method
00127   // must be called after all phases are added, before doing
00128   // anything else with the mixture. After init() has been called,
00129   // no more phases may be added.
00130   void MultiPhase::init() {
00131     if (m_init) return;
00132     index_t ip, kp, k = 0, nsp, m;
00133     int mlocal;
00134     string sym;
00135 
00136     // allocate space for the atomic composition matrix
00137     m_atoms.resize(m_nel, m_nsp, 0.0);
00138     m_moleFractions.resize(m_nsp, 0.0);
00139     m_elemAbundances.resize(m_nel, 0.0);
00140 
00141     // iterate over the elements
00142     //   -> fill in m_atoms(m,k), m_snames(k), m_spphase(k),
00143     //              m_sptart(ip)
00144     for (m = 0; m < m_nel; m++) {
00145       sym = m_enames[m];
00146       k = 0;
00147       // iterate over the phases
00148       for (ip = 0; ip < m_np; ip++) {
00149         phase_t* p = m_phase[ip];
00150         nsp = p->nSpecies();
00151         mlocal = p->elementIndex(sym);    
00152         for (kp = 0; kp < nsp; kp++) {
00153           if (mlocal >= 0) {
00154             m_atoms(m, k) = p->nAtoms(kp, mlocal);
00155           }
00156           if (m == 0) {
00157             m_snames.push_back(p->speciesName(kp));
00158             if (kp == 0) {
00159               m_spstart.push_back(m_spphase.size());
00160             }
00161             m_spphase.push_back(ip);
00162           }
00163           k++;
00164         }
00165       }
00166     }
00167 
00168     if (m_eloc >= 0) {
00169       doublereal esum;
00170       for (k = 0; k < m_nsp; k++) {
00171         esum = 0.0;
00172         for (m = 0; m < m_nel; m++) {
00173           if (int(m) != m_eloc)
00174             esum += m_atoms(m,k) * m_atomicNumber[m];
00175         }
00176         //m_atoms(m_eloc, k) += esum;
00177       }
00178     }
00179 
00180     /// set the initial composition within each phase to the
00181     /// mole fractions stored in the phase objects
00182     m_init = true;
00183     
00184     updateMoleFractions();
00185 
00186   }
00187 
00188 
00189   // Return a reference to phase n. The state of phase n is
00190   // also updated to match the state stored locally in the 
00191   // mixture object.
00192   MultiPhase::phase_t& MultiPhase::phase(index_t n) {
00193     if (!m_init) init();
00194     m_phase[n]->setTemperature(m_temp);
00195     m_phase[n]->setMoleFractions_NoNorm(DATA_PTR(m_moleFractions) + m_spstart[n]);
00196     m_phase[n]->setPressure(m_press);
00197     return *m_phase[n];
00198   }
00199 
00200   /// Moles of species \c k.
00201   doublereal MultiPhase::speciesMoles(index_t k) const {
00202     index_t ip = m_spphase[k];
00203     return m_moles[ip]*m_moleFractions[k];
00204   }
00205 
00206   /// Total moles of element m, summed over all
00207   /// phases
00208   doublereal MultiPhase::elementMoles(index_t m) const {
00209     doublereal sum = 0.0, phasesum;
00210     index_t i, k = 0, ik, nsp;
00211     for (i = 0; i < m_np; i++) {
00212       phasesum = 0.0;
00213       nsp = m_phase[i]->nSpecies();
00214       for (ik = 0; ik < nsp; ik++) {
00215         k = speciesIndex(ik, i);
00216         phasesum += m_atoms(m,k)*m_moleFractions[k];
00217       }
00218       sum += phasesum * m_moles[i];
00219     }
00220     return sum;
00221   }
00222 
00223   /// Total charge, summed over all phases
00224   doublereal MultiPhase::charge() const {
00225     doublereal sum = 0.0;
00226     index_t i;
00227     for (i = 0; i < m_np; i++) {
00228       sum += phaseCharge(i);
00229     }
00230     return sum;
00231   }
00232 
00233   int MultiPhase::speciesIndex(std::string speciesName, std::string phaseName) {
00234     int p = phaseIndex(phaseName);
00235     if (p < 0) {
00236       throw CanteraError("MultiPhase::speciesIndex", "phase not found: " + phaseName);
00237     }
00238     int k = m_phase[p]->speciesIndex(speciesName);
00239     if (k < 0) {
00240       throw CanteraError("MultiPhase::speciesIndex", "species not found: " + speciesName);
00241     }
00242     return m_spstart[p] + k;
00243   }
00244 
00245   /// Net charge of one phase (Coulombs). The net charge is computed as
00246   /// \f[ Q_p = N_p \sum_k F z_k X_k \f]
00247   /// where the sum runs only over species in phase \a p.
00248   /// @param p index of the phase for which the charge is desired.   
00249   doublereal MultiPhase::phaseCharge(index_t p) const {
00250     doublereal phasesum = 0.0;
00251     int ik, k, nsp = m_phase[p]->nSpecies();
00252     for (ik = 0; ik < nsp; ik++) {
00253       k = speciesIndex(ik, p);
00254       phasesum += m_phase[p]->charge(ik)*m_moleFractions[k];
00255     }
00256     return Faraday*phasesum*m_moles[p];
00257   }
00258 
00259 
00260   /// Get the chemical potentials of all species in all phases. 
00261   void MultiPhase::getChemPotentials(doublereal* mu) const {
00262     index_t i, loc = 0;
00263     updatePhases();            
00264     for (i = 0; i < m_np; i++) {
00265       m_phase[i]->getChemPotentials(mu + loc);
00266       loc += m_phase[i]->nSpecies();
00267     }
00268   }
00269 
00270   // Get chemical potentials of species with valid thermo
00271   // data. This method is designed for use in computing chemical
00272   // equilibrium by Gibbs minimization. For solution phases (more
00273   // than one species), this does the same thing as
00274   // getChemPotentials. But for stoichiometric phases, this writes
00275   // into array \a mu the user-specified value \a not_mu instead of
00276   // the chemical potential if the temperature is outside the range
00277   // for which the thermo data for the one species in the phase are
00278   // valid. The need for this arises since many condensed phases
00279   // have thermo data fit only for the temperature range for which
00280   // they are stable. For example, in the NASA database, the fits
00281   // for H2O(s) are only done up to 0 C, the fits for H2O(L) are
00282   // only done from 0 C to 100 C, etc. Using the polynomial fits outside
00283   // the range for which the fits were done can result in spurious
00284   // chemical potentials, and can lead to condensed phases
00285   // appearing when in fact they should be absent.
00286   //
00287   // By setting \a not_mu to a large positive value, it is possible
00288   // to force routines which seek to minimize the Gibbs free energy
00289   // of the mixture to zero out any phases outside the temperature
00290   // range for which their thermo data are valid.
00291   //
00292   // If this method is called with \a standard set to true, then
00293   // the composition-independent standard chemical potentials are
00294   // returned instead of the composition-dependent chemical
00295   // potentials.
00296   //
00297   void MultiPhase::getValidChemPotentials(doublereal not_mu,
00298                                           doublereal* mu, bool standard) const {
00299     index_t i, loc = 0;
00300 
00301     updatePhases();           
00302     // iterate over the phases 
00303     for (i = 0; i < m_np; i++) {
00304       if (tempOK(i) || m_phase[i]->nSpecies() > 1) {
00305         if (!standard)
00306           m_phase[i]->getChemPotentials(mu + loc);
00307         else
00308           m_phase[i]->getStandardChemPotentials(mu + loc);
00309       }
00310       else
00311         fill(mu + loc, mu + loc + m_phase[i]->nSpecies(), not_mu);
00312       loc += m_phase[i]->nSpecies();
00313     }
00314   }
00315 
00316   /// True if species \a k belongs to a solution phase.
00317   bool MultiPhase::solutionSpecies(index_t k) const {
00318     if (m_phase[m_spphase[k]]->nSpecies() > 1)
00319       return true;
00320     else
00321       return false;
00322   }
00323 
00324   /// The Gibbs free energy of the mixture (J).
00325   doublereal MultiPhase::gibbs() const {
00326     index_t i;
00327     doublereal sum = 0.0;
00328     updatePhases();
00329     for (i = 0; i < m_np; i++) 
00330       sum += m_phase[i]->gibbs_mole() * m_moles[i];
00331     return sum;
00332   }
00333 
00334   /// The enthalpy of the mixture (J).
00335   doublereal MultiPhase::enthalpy() const {
00336     index_t i;
00337     doublereal sum = 0.0;
00338     updatePhases();
00339     for (i = 0; i < m_np; i++) 
00340       sum += m_phase[i]->enthalpy_mole() * m_moles[i];
00341     return sum;
00342   }
00343 
00344   /// The internal energy of the mixture (J).
00345   doublereal MultiPhase::IntEnergy() const {
00346     index_t i;
00347     doublereal sum = 0.0;
00348     updatePhases();
00349     for (i = 0; i < m_np; i++) 
00350       sum += m_phase[i]->intEnergy_mole() * m_moles[i];
00351     return sum;
00352   }
00353 
00354   /// The entropy of the mixture (J/K).
00355   doublereal MultiPhase::entropy() const {
00356     index_t i;
00357     doublereal sum = 0.0;
00358     updatePhases();
00359     for (i = 0; i < m_np; i++) 
00360       sum += m_phase[i]->entropy_mole() * m_moles[i];
00361     return sum;
00362   }
00363 
00364   /// The specific heat at constant pressure and composition (J/K).
00365   /// Note that this does not account for changes in composition of
00366   /// the mixture with temperature.
00367   doublereal MultiPhase::cp() const {
00368     index_t i;
00369     doublereal sum = 0.0;
00370     updatePhases();
00371     for (i = 0; i < m_np; i++) 
00372       sum += m_phase[i]->cp_mole() * m_moles[i];
00373     return sum;
00374   }
00375 
00376 
00377 
00378   /// Set the mole fractions of phase \a n to the values in 
00379   /// array \a x.  
00380   void MultiPhase::setPhaseMoleFractions(const index_t n, const doublereal* const x) {
00381     phase_t* p = m_phase[n];
00382     p->setState_TPX(m_temp, m_press, x);
00383     int nsp = p->nSpecies();
00384     int istart = m_spstart[n];
00385     for (int k = 0; k < nsp; k++) {
00386       m_moleFractions[istart+k] = x[k];
00387     }
00388   }
00389 
00390   // Set the species moles using a map. The map \a xMap maps
00391   // species name strings to mole numbers. Mole numbers that are
00392   // less than or equal to zero will be set to zero.
00393   void MultiPhase::setMolesByName(compositionMap& xMap) {
00394     int kk = nSpecies();
00395     doublereal x;
00396     vector_fp moles(kk, 0.0);
00397     for (int k = 0; k < kk; k++) {
00398       x = xMap[speciesName(k)];
00399       if (x > 0.0) moles[k] = x;
00400     }
00401     setMoles(DATA_PTR(moles));
00402   }
00403 
00404   // Set the species moles using a string. Unspecified species are
00405   // set to zero.
00406   void MultiPhase::setMolesByName(const std::string& x) {
00407     compositionMap xx;
00408 
00409     // add an entry in the map for every species, with value -1.0.
00410     // Function parseCompString (stringUtils.cpp) uses the names
00411     // in the map to specify the allowed species.
00412     int kk = nSpecies();
00413     for (int k = 0; k < kk; k++) { 
00414       xx[speciesName(k)] = -1.0;
00415     }
00416 
00417     // build the composition map from the string, and then set the
00418     // moles.
00419     parseCompString(x, xx);
00420     setMolesByName(xx); 
00421   }
00422  
00423   // Get the mole numbers of all species in the multiphase
00424   // object
00425   void MultiPhase::getMoles(doublereal * molNum) const {
00426     /*
00427      * First copy in the mole fractions
00428      */
00429     copy(m_moleFractions.begin(), m_moleFractions.end(), molNum);
00430     index_t ik;
00431     doublereal *dtmp = molNum;
00432     for (index_t ip = 0; ip < m_np; ip++) {
00433       doublereal phasemoles = m_moles[ip];
00434       phase_t* p = m_phase[ip];
00435       index_t nsp = p->nSpecies();
00436       for (ik = 0; ik < nsp; ik++) {
00437         *(dtmp++) *= phasemoles;
00438       }
00439     }
00440   }
00441 
00442   /// Set the species moles to the values in array \a n. The state
00443   /// of each phase object is also updated to have the specified
00444   /// composition and the mixture temperature and pressure.
00445   void MultiPhase::setMoles(const doublereal* n) {
00446     if (!m_init) init();
00447     index_t ip, loc = 0;
00448     index_t ik, k = 0, nsp;
00449     doublereal phasemoles;
00450     for (ip = 0; ip < m_np; ip++) {
00451       phase_t* p = m_phase[ip];
00452       nsp = p->nSpecies();
00453       phasemoles = 0.0;
00454       for (ik = 0; ik < nsp; ik++) {
00455         phasemoles += n[k];
00456         k++;
00457       }
00458       m_moles[ip] = phasemoles;
00459       if (nsp > 1) {
00460         if (phasemoles > 0.0) {
00461           p->setState_TPX(m_temp, m_press, n + loc);
00462           p->getMoleFractions(DATA_PTR(m_moleFractions) + loc);
00463         } else {
00464           p->getMoleFractions(DATA_PTR(m_moleFractions) + loc);
00465         }
00466       }
00467       else {
00468         m_moleFractions[loc] = 1.0;
00469       }
00470       loc += nsp;
00471     }
00472   }
00473 
00474   void MultiPhase::addSpeciesMoles(const int indexS, const doublereal addedMoles) {
00475     vector_fp tmpMoles(m_nsp, 0.0);
00476     getMoles(DATA_PTR(tmpMoles));
00477     tmpMoles[indexS] += addedMoles;
00478     if (tmpMoles[indexS] < 0.0) {
00479       tmpMoles[indexS] = 0.0;
00480     }
00481     setMoles(DATA_PTR(tmpMoles));
00482   }
00483   
00484   void MultiPhase::setState_TP(const doublereal T, const doublereal Pres) {
00485     if (!m_init) init();
00486     m_temp  = T;
00487     m_press = Pres;
00488     updatePhases();
00489   }
00490 
00491   void MultiPhase::setState_TPMoles(const doublereal T, const doublereal Pres,
00492                                     const doublereal *n) {
00493     m_temp  = T;
00494     m_press = Pres;
00495     setMoles(n);       
00496   }
00497 
00498   void MultiPhase::getElemAbundances(doublereal *elemAbundances) const {
00499     index_t eGlobal;
00500     calcElemAbundances();
00501     for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
00502       elemAbundances[eGlobal] = m_elemAbundances[eGlobal];
00503     }
00504   }
00505 
00506   // Internal routine to calculate the element abundance vector
00507   void MultiPhase::calcElemAbundances() const {
00508     index_t loc = 0;
00509     index_t eGlobal;
00510     int ik, kGlobal;
00511     doublereal spMoles;
00512     for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
00513       m_elemAbundances[eGlobal] = 0.0;
00514     }
00515     for (index_t ip = 0; ip < m_np; ip++) {
00516       phase_t* p = m_phase[ip];
00517       int nspPhase = p->nSpecies();
00518       doublereal phasemoles = m_moles[ip];
00519       for (ik = 0; ik < nspPhase; ik++) {
00520         kGlobal = loc + ik;
00521         spMoles = m_moleFractions[kGlobal] * phasemoles;
00522         for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
00523           m_elemAbundances[eGlobal] += m_atoms(eGlobal, kGlobal) * spMoles;
00524         }
00525       }
00526       loc += nspPhase;
00527     }
00528   }
00529 
00530   /// The total mixture volume [m^3].
00531   doublereal MultiPhase::volume() const {
00532     int i;
00533     doublereal sum = 0;
00534     for (i = 0; i < int(m_np); i++) {
00535       sum += m_moles[i]/m_phase[i]->molarDensity();
00536     }
00537     return sum;
00538   }
00539 
00540   doublereal MultiPhase::equilibrate(int XY, doublereal err, 
00541                                      int maxsteps, int maxiter, int loglevel) {
00542     doublereal error;
00543     bool strt = false;
00544     doublereal dt;
00545     doublereal h0;
00546     int n;
00547     bool start;
00548     doublereal ferr, hnow, herr = 1.0;
00549     doublereal snow, serr = 1.0, s0;
00550     doublereal Tlow = -1.0, Thigh = -1.0;
00551     doublereal Hlow = Undef, Hhigh = Undef, tnew;
00552     doublereal dta=0.0, dtmax, cpb;
00553     MultiPhaseEquil* e = 0;
00554 
00555     if (!m_init) init();
00556     if (loglevel > 0)
00557       beginLogGroup("MultiPhase::equilibrate", loglevel);
00558 
00559     if (XY == TP) {
00560       if (loglevel > 0) {
00561         addLogEntry("problem type","fixed T,P");
00562         addLogEntry("Temperature",temperature());
00563         addLogEntry("Pressure", pressure());
00564       }
00565 
00566       // create an equilibrium manager 
00567       e = new MultiPhaseEquil(this);
00568       try {
00569         error = e->equilibrate(XY, err, maxsteps, loglevel);
00570       }
00571       catch (CanteraError &err) {
00572         if (loglevel > 0)
00573           endLogGroup();
00574         delete e;
00575         e = 0;
00576         throw err;
00577       }
00578       goto done;
00579     }
00580 
00581     else if (XY == HP) {
00582       h0 = enthalpy();
00583       Tlow = 0.5*m_Tmin;      // lower bound on T
00584       Thigh = 2.0*m_Tmax;     // upper bound on T
00585       if (loglevel > 0) {
00586         addLogEntry("problem type","fixed H,P");
00587         addLogEntry("H target",fp2str(h0));
00588       }
00589       for (n = 0; n < maxiter; n++) {
00590 
00591         // if 'strt' is false, the current composition will be used as
00592         // the starting estimate; otherwise it will be estimated
00593         //                if (e) {
00594         //    cout << "e should be zero, but it is not!" << endl;
00595         //    delete e;
00596         // }
00597         e = new MultiPhaseEquil(this, strt);
00598         // start with a loose error tolerance, but tighten it as we get 
00599         // close to the final temperature
00600         if (loglevel > 0)   
00601           beginLogGroup("iteration "+int2str(n));
00602 
00603         try {
00604           error = e->equilibrate(TP, err, maxsteps, loglevel);
00605           hnow = enthalpy();
00606           // the equilibrium enthalpy monotonically increases with T; 
00607           // if the current value is below the target, the we know the
00608           // current temperature is too low. Set 
00609           if (hnow < h0) {
00610             if (m_temp > Tlow) {
00611               Tlow = m_temp;
00612               Hlow = hnow;
00613             }
00614           }
00615           // the current enthalpy is greater than the target; therefore the
00616           // current temperature is too high.
00617           else {
00618             if (m_temp < Thigh) {
00619               Thigh = m_temp;
00620               Hhigh = hnow;
00621             }
00622           }
00623           if (Hlow != Undef && Hhigh != Undef) {
00624             cpb = (Hhigh - Hlow)/(Thigh - Tlow);
00625             dt = (h0 - hnow)/cpb;    
00626             dta = fabs(dt);
00627             dtmax = 0.5*fabs(Thigh - Tlow);
00628             if (dta > dtmax) dt *= dtmax/dta;  
00629           }
00630           else {
00631             tnew = sqrt(Tlow*Thigh);
00632             dt = tnew - m_temp;
00633             //cpb = cp();
00634           }
00635 
00636           herr = fabs((h0 - hnow)/h0);
00637           if (loglevel > 0) {
00638             addLogEntry("T",fp2str(temperature()));
00639             addLogEntry("H",fp2str(hnow));
00640             addLogEntry("H rel error",fp2str(herr));
00641             addLogEntry("lower T bound",fp2str(Tlow));
00642             addLogEntry("upper T bound",fp2str(Thigh));
00643             endLogGroup(); // iteration
00644           }
00645 
00646 
00647           if (herr < err) { // || dta < 1.0e-4) {
00648             if (loglevel > 0) {
00649               addLogEntry("T iterations",int2str(n));
00650               addLogEntry("Final T",fp2str(temperature()));
00651               addLogEntry("H rel error",fp2str(herr));
00652             }
00653             goto done;
00654           }
00655           tnew = m_temp + dt;
00656           if (tnew < 0.0) tnew = 0.5*m_temp;
00657           //dta = fabs(tnew - m_temp);
00658           setTemperature(tnew);
00659 
00660           // if the size of Delta T is not too large, use
00661           // the current composition as the starting estimate
00662           if (dta < 100.0) strt = false;
00663 
00664         }
00665 
00666         catch (CanteraError err) {
00667           if (!strt) {        
00668             if (loglevel > 0)
00669               addLogEntry("no convergence",
00670                           "try estimating starting composition");
00671             strt = true;
00672           }
00673           else {
00674             tnew = 0.5*(m_temp + Thigh);
00675             if (fabs(tnew - m_temp) < 1.0) tnew = m_temp + 1.0;
00676             setTemperature(tnew);
00677             if (loglevel > 0)
00678               addLogEntry("no convergence",
00679                           "trying T = "+fp2str(m_temp));
00680           }
00681           if (loglevel > 0)
00682             endLogGroup();
00683         }
00684         delete e;
00685         e = 0;
00686       }
00687       if (loglevel > 0) {
00688         addLogEntry("reached max number of T iterations",int2str(maxiter));
00689         endLogGroup();
00690       }
00691       throw CanteraError("MultiPhase::equilibrate",
00692                          "No convergence for T");
00693     }
00694     else if (XY == SP) {
00695       s0 = entropy();
00696       start = true;
00697       Tlow = 1.0; // m_Tmin;      // lower bound on T
00698       Thigh = 1.0e6; // m_Tmax;   // upper bound on T
00699       if (loglevel > 0) {
00700         addLogEntry("problem type","fixed S,P");
00701         addLogEntry("S target",fp2str(s0));
00702         addLogEntry("min T",fp2str(Tlow));
00703         addLogEntry("max T",fp2str(Thigh));
00704       }
00705       for (n = 0; n < maxiter; n++) {
00706         if (e) delete e;
00707         e = new MultiPhaseEquil(this, strt);
00708         ferr = 0.1;
00709         if (fabs(dt) < 1.0) ferr = err;
00710         //start = false;
00711         if (loglevel > 0)
00712           beginLogGroup("iteration "+int2str(n));
00713                 
00714         try {
00715           error = e->equilibrate(TP, err, maxsteps, loglevel);
00716           snow = entropy();
00717           if (snow < s0) {
00718             if (m_temp > Tlow) Tlow = m_temp;
00719           }
00720           else {
00721             if (m_temp < Thigh) Thigh = m_temp;
00722           }
00723           serr = fabs((s0 - snow)/s0);
00724           if (loglevel > 0) {
00725             addLogEntry("T",fp2str(temperature()));
00726             addLogEntry("S",fp2str(snow));
00727             addLogEntry("S rel error",fp2str(serr));
00728             endLogGroup();
00729           }
00730           dt = (s0 - snow)*m_temp/cp();
00731           dtmax = 0.5*fabs(Thigh - Tlow);
00732           dtmax = (dtmax > 500.0 ? 500.0 : dtmax);
00733           dta = fabs(dt);
00734           if (dta > dtmax) dt *= dtmax/dta;
00735           if (herr < err || dta < 1.0e-4) {
00736             if (loglevel > 0) {
00737               addLogEntry("T iterations",int2str(n));
00738               addLogEntry("Final T",fp2str(temperature()));
00739               addLogEntry("S rel error",fp2str(serr));
00740             }
00741             goto done;
00742           }
00743           tnew = m_temp + dt;
00744           setTemperature(tnew);
00745 
00746           // if the size of Delta T is not too large, use
00747           // the current composition as the starting estimate
00748           if (dta < 100.0) strt = false;
00749         }
00750 
00751         catch (CanteraError err) {
00752           if (!strt) {
00753             if (loglevel > 0) {
00754               addLogEntry("no convergence",
00755                           "setting strt to True");
00756             }
00757             strt = true;
00758           }
00759           else {
00760             tnew = 0.5*(m_temp + Thigh);
00761             setTemperature(tnew);
00762             if (loglevel > 0) {
00763               addLogEntry("no convergence",
00764                           "trying T = "+fp2str(m_temp));
00765             }
00766           }
00767           if (loglevel > 0)
00768             endLogGroup();
00769         }
00770         delete e;
00771         e = 0;
00772       }
00773       if (loglevel > 0) {
00774         addLogEntry("reached max number of T iterations",int2str(maxiter));
00775         endLogGroup();
00776       }
00777       throw CanteraError("MultiPhase::equilibrate",
00778                          "No convergence for T");
00779     }
00780     else if (XY == TV) {
00781       addLogEntry("problem type","fixed T, V");
00782       //            doublereal dt = 1.0e3;
00783       doublereal v0 = volume();
00784       doublereal dVdP;
00785       int n;
00786       bool start = true;
00787       doublereal error, vnow, pnow, verr;
00788       for (n = 0; n < maxiter; n++) {
00789         pnow = pressure();
00790         MultiPhaseEquil e(this, start);
00791         start = false;
00792         beginLogGroup("iteration "+int2str(n));
00793                 
00794         error = e.equilibrate(TP, err, maxsteps, loglevel);
00795         vnow = volume();
00796         verr = fabs((v0 - vnow)/v0);
00797         addLogEntry("P",fp2str(pressure()));
00798         addLogEntry("V rel error",fp2str(verr));
00799         endLogGroup();
00800                 
00801         if (verr < err) {
00802           addLogEntry("P iterations",int2str(n));
00803           addLogEntry("Final P",fp2str(pressure()));
00804           addLogEntry("V rel error",fp2str(verr));
00805           goto done;
00806         }
00807         // find dV/dP
00808         setPressure(pnow*1.01);
00809         dVdP = (volume() - vnow)/(0.01*pnow);
00810         setPressure(pnow + 0.5*(v0 - vnow)/dVdP);
00811       }
00812     }
00813 
00814     else {
00815       if (loglevel > 0)
00816         endLogGroup();
00817       throw CanteraError("MultiPhase::equilibrate","unknown option");
00818     }
00819     return -1.0;
00820   done:
00821     delete e;
00822     e = 0;
00823     if (loglevel > 0)
00824       endLogGroup();
00825     return err;
00826   }
00827 
00828 #ifdef MULTIPHASE_DEVEL
00829   void importFromXML(string infile, string id) {
00830     XML_Node* root = get_XML_File(infile);
00831     if (id == "-") id = "";
00832     XML_Node* x = get_XML_Node(string("#")+id, root);
00833     if (x.name() != "multiphase")
00834       throw CanteraError("MultiPhase::importFromXML",
00835                          "Current XML_Node is not a multiphase element.");
00836     vector<XML_Node*> phases;
00837     x.getChildren("phase",phases);
00838     int np = phases.size();
00839     int n;
00840     ThermoPhase* p;
00841     for (n = 0; n < np; n++) {
00842       XML_Node& ph = *phases[n];
00843       srcfile = infile;
00844       if (ph.hasAttrib("src")) srcfile = ph["src"];
00845       idstr =  ph["id"];
00846       p = newPhase(srcfile, idstr);
00847       if (p) {
00848         addPhase(p, ph.value());
00849       }
00850     }
00851   }
00852 #endif
00853 
00854   void MultiPhase::setTemperature(const doublereal T) {
00855     if (!m_init) init();
00856     m_temp = T;
00857     updatePhases();
00858   }
00859 
00860   // Name of element \a m.
00861   std::string MultiPhase::elementName(int m) const {
00862     return m_enames[m];
00863   }
00864 
00865   // Index of element with name \a name.
00866   int MultiPhase::elementIndex(std::string name) const {
00867     for (size_t e = 0; e < m_nel; e++) {
00868       if (m_enames[e] == name) {
00869         return (int) e;
00870       }
00871     }
00872     return -1;
00873   }
00874   
00875   // Name of species with global index \a k.
00876   std::string MultiPhase::speciesName(const int k) const {
00877     return m_snames[k]; 
00878   }
00879 
00880   doublereal MultiPhase::nAtoms(const int kGlob, const int mGlob) const {
00881     return m_atoms(mGlob, kGlob);
00882   }
00883 
00884   void MultiPhase::getMoleFractions(doublereal* const x) const {
00885     std::copy(m_moleFractions.begin(), m_moleFractions.end(), x);
00886   }
00887 
00888   std::string MultiPhase::phaseName(const index_t iph) const {
00889     const phase_t *tptr = m_phase[iph];
00890     return tptr->id();
00891   }
00892 
00893   int MultiPhase::phaseIndex(const std::string &pName) const {
00894     std::string tmp;
00895     for (int iph = 0; iph < (int) m_np; iph++) {
00896       const phase_t *tptr = m_phase[iph];
00897       tmp = tptr->id();
00898       if (tmp == pName) {
00899         return iph;
00900       }
00901     }
00902     return -1;
00903   }
00904 
00905   doublereal MultiPhase::phaseMoles(const index_t n) const {
00906     return m_moles[n];
00907   }
00908 
00909   void MultiPhase::setPhaseMoles(const index_t n, const doublereal moles) {
00910     m_moles[n] = moles;
00911   }
00912 
00913   int MultiPhase::speciesPhaseIndex(const index_t kGlob) const {
00914     return m_spphase[kGlob];
00915   }
00916 
00917   doublereal MultiPhase::moleFraction(const index_t kGlob) const{
00918     return m_moleFractions[kGlob];
00919   }
00920 
00921 
00922   bool MultiPhase::tempOK(const index_t p) const {
00923     return m_temp_OK[p];
00924   }
00925 
00926   /// Update the locally-stored species mole fractions. 
00927   void MultiPhase::updateMoleFractions() {
00928     uploadMoleFractionsFromPhases();
00929   }
00930   /// Update the locally-stored species mole fractions. 
00931   void MultiPhase::uploadMoleFractionsFromPhases() {
00932     index_t ip, loc = 0;
00933     for (ip = 0; ip < m_np; ip++) {
00934       phase_t* p = m_phase[ip];
00935       p->getMoleFractions(DATA_PTR(m_moleFractions) + loc);
00936       loc += p->nSpecies();
00937     }
00938     calcElemAbundances();
00939   }
00940 
00941   //-------------------------------------------------------------
00942   //
00943   // protected methods
00944   //
00945   //-------------------------------------------------------------
00946 
00947 
00948 
00949   /// synchronize the phase objects with the mixture state. This
00950   /// method sets each phase to the mixture temperature and
00951   /// pressure, and sets the phase mole fractions based on the
00952   /// mixture mole numbers. 
00953   void MultiPhase::updatePhases() const {
00954     index_t p, nsp, loc = 0;
00955     for (p = 0; p < m_np; p++) {
00956       nsp = m_phase[p]->nSpecies();
00957       const doublereal* x = DATA_PTR(m_moleFractions) + loc;
00958       loc += nsp;
00959       m_phase[p]->setState_TPX(m_temp, m_press, x);
00960       m_temp_OK[p] = true;
00961       if (m_temp < m_phase[p]->minTemp() 
00962           || m_temp > m_phase[p]->maxTemp()) {
00963         m_temp_OK[p] = false;
00964       }
00965     }
00966   }            
00967 
00968 }
00969 
Generated by  doxygen 1.6.3