25 for (
size_t n = 0; n < mix.
nPhases(); n++) {
30 void MultiPhase::addPhases(vector<ThermoPhase*>& phases,
31 const vector<double>& phaseMoles)
33 for (
size_t n = 0; n < phases.size(); n++) {
34 addPhase(phases[n], phaseMoles[n]);
43 "phases cannot be added after init() has been called.");
47 throw CanteraError(
"MultiPhase::addPhase",
"Phase '{}'' is not "
48 "compatible with MultiPhase equilibrium solver", p->
name());
55 m_moles.push_back(moles);
56 m_temp_OK.push_back(
true);
65 for (
size_t m = 0; m < p->
nElements(); m++) {
71 if (m_enamemap.find(ename) == m_enamemap.end()) {
72 m_enamemap[ename] = m_nel + 1;
73 m_enames.push_back(ename);
77 if (ename ==
"E" || ename ==
"e") {
88 if (m_temp == 298.15 && p->
temperature() > 2.0E-3) {
101 m_Tmin = std::max(p->
minTemp(), m_Tmin);
102 m_Tmax = std::min(p->
maxTemp(), m_Tmax);
106 void MultiPhase::init()
113 m_atoms.resize(m_nel, m_nsp, 0.0);
114 m_moleFractions.resize(m_nsp, 0.0);
115 m_elemAbundances.resize(m_nel, 0.0);
119 for (
size_t m = 0; m < m_nel; m++) {
122 for (
size_t ip = 0; ip < nPhases(); ip++) {
126 for (
size_t kp = 0; kp < nsp; kp++) {
127 if (mlocal !=
npos) {
128 m_atoms(m, k) = p->
nAtoms(kp, mlocal);
133 m_spstart.push_back(m_spphase.size());
135 m_spphase.push_back(ip);
145 uploadMoleFractionsFromPhases();
154 m_phase[n]->setTemperature(m_temp);
155 m_phase[n]->setMoleFractions_NoNorm(&m_moleFractions[m_spstart[n]]);
156 m_phase[n]->setPressure(m_press);
160 void MultiPhase::checkPhaseIndex(
size_t m)
const
162 if (m >= nPhases()) {
163 throw IndexError(
"MultiPhase::checkPhaseIndex",
"phase", m, nPhases()-1);
167 void MultiPhase::checkPhaseArraySize(
size_t mm)
const
169 if (nPhases() > mm) {
170 throw ArraySizeError(
"MultiPhase::checkPhaseIndex", mm, nPhases());
174 double MultiPhase::speciesMoles(
size_t k)
const
176 size_t ip = m_spphase[k];
177 return m_moles[ip]*m_moleFractions[k];
180 double MultiPhase::elementMoles(
size_t m)
const
183 for (
size_t i = 0; i < nPhases(); i++) {
184 double phasesum = 0.0;
185 size_t nsp = m_phase[i]->nSpecies();
186 for (
size_t ik = 0; ik < nsp; ik++) {
187 size_t k = speciesIndex(ik, i);
188 phasesum += m_atoms(m,k)*m_moleFractions[k];
190 sum += phasesum * m_moles[i];
195 double MultiPhase::charge()
const
198 for (
size_t i = 0; i < nPhases(); i++) {
199 sum += phaseCharge(i);
204 size_t MultiPhase::speciesIndex(
const string& speciesName,
const string& phaseName)
209 size_t p = phaseIndex(phaseName);
211 throw CanteraError(
"MultiPhase::speciesIndex",
"phase not found: " + phaseName);
213 size_t k = m_phase[p]->speciesIndex(speciesName);
215 throw CanteraError(
"MultiPhase::speciesIndex",
"species not found: " + speciesName);
217 return m_spstart[p] + k;
220 double MultiPhase::phaseCharge(
size_t p)
const
222 double phasesum = 0.0;
223 size_t nsp = m_phase[p]->nSpecies();
224 for (
size_t ik = 0; ik < nsp; ik++) {
225 size_t k = speciesIndex(ik, p);
226 phasesum += m_phase[p]->charge(ik)*m_moleFractions[k];
228 return Faraday*phasesum*m_moles[p];
231 void MultiPhase::getChemPotentials(
double* mu)
const
235 for (
size_t i = 0; i < nPhases(); i++) {
236 m_phase[i]->getChemPotentials(mu + loc);
237 loc += m_phase[i]->nSpecies();
241 void MultiPhase::getValidChemPotentials(
double not_mu,
double* mu,
bool standard)
const
246 for (
size_t i = 0; i < nPhases(); i++) {
247 if (tempOK(i) || m_phase[i]->nSpecies() > 1) {
249 m_phase[i]->getChemPotentials(mu + loc);
251 m_phase[i]->getStandardChemPotentials(mu + loc);
254 fill(mu + loc, mu + loc + m_phase[i]->nSpecies(), not_mu);
256 loc += m_phase[i]->nSpecies();
260 bool MultiPhase::solutionSpecies(
size_t k)
const
262 if (m_phase[m_spphase[k]]->nSpecies() > 1) {
269 double MultiPhase::gibbs()
const
273 for (
size_t i = 0; i < nPhases(); i++) {
274 if (m_moles[i] > 0.0) {
275 sum += m_phase[i]->gibbs_mole() * m_moles[i];
281 double MultiPhase::enthalpy()
const
285 for (
size_t i = 0; i < nPhases(); i++) {
286 if (m_moles[i] > 0.0) {
287 sum += m_phase[i]->enthalpy_mole() * m_moles[i];
293 double MultiPhase::IntEnergy()
const
297 for (
size_t i = 0; i < nPhases(); i++) {
298 if (m_moles[i] > 0.0) {
299 sum += m_phase[i]->intEnergy_mole() * m_moles[i];
305 double MultiPhase::entropy()
const
309 for (
size_t i = 0; i < nPhases(); i++) {
310 if (m_moles[i] > 0.0) {
311 sum += m_phase[i]->entropy_mole() * m_moles[i];
317 double MultiPhase::cp()
const
321 for (
size_t i = 0; i < nPhases(); i++) {
322 if (m_moles[i] > 0.0) {
323 sum += m_phase[i]->cp_mole() * m_moles[i];
329 void MultiPhase::setPhaseMoleFractions(
const size_t n,
const double*
const x)
336 size_t istart = m_spstart[n];
337 for (
size_t k = 0; k < p->
nSpecies(); k++) {
338 m_moleFractions[istart+k] = x[k];
344 size_t kk = nSpecies();
345 vector<double> moles(kk, 0.0);
346 for (
size_t k = 0; k < kk; k++) {
347 moles[k] = std::max(
getValue(xMap, speciesName(k), 0.0), 0.0);
349 setMoles(moles.data());
352 void MultiPhase::setMolesByName(
const string& x)
359 void MultiPhase::getMoles(
double* molNum)
const
362 copy(m_moleFractions.begin(), m_moleFractions.end(), molNum);
363 double* dtmp = molNum;
364 for (
size_t ip = 0; ip < nPhases(); ip++) {
365 double phasemoles = m_moles[ip];
368 for (
size_t ik = 0; ik < nsp; ik++) {
369 *(dtmp++) *= phasemoles;
374 void MultiPhase::setMoles(
const double* n)
381 for (
size_t ip = 0; ip < nPhases(); ip++) {
384 double phasemoles = 0.0;
385 for (
size_t ik = 0; ik < nsp; ik++) {
389 m_moles[ip] = phasemoles;
391 if (phasemoles > 0.0) {
398 m_moleFractions[loc] = 1.0;
404 void MultiPhase::addSpeciesMoles(
const int indexS,
const double addedMoles)
406 vector<double> tmpMoles(m_nsp, 0.0);
407 getMoles(tmpMoles.data());
408 tmpMoles[indexS] += addedMoles;
409 tmpMoles[indexS] = std::max(tmpMoles[indexS], 0.0);
410 setMoles(tmpMoles.data());
413 void MultiPhase::setState_TP(
const double T,
const double Pres)
423 void MultiPhase::setState_TPMoles(
const double T,
const double Pres,
const double* n)
430 void MultiPhase::getElemAbundances(
double* elemAbundances)
const
432 calcElemAbundances();
433 for (
size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) {
434 elemAbundances[eGlobal] = m_elemAbundances[eGlobal];
438 void MultiPhase::calcElemAbundances()
const
442 for (
size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) {
443 m_elemAbundances[eGlobal] = 0.0;
445 for (
size_t ip = 0; ip < nPhases(); ip++) {
448 double phasemoles = m_moles[ip];
449 for (
size_t ik = 0; ik < nspPhase; ik++) {
450 size_t kGlobal = loc + ik;
451 spMoles = m_moleFractions[kGlobal] * phasemoles;
452 for (
size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) {
453 m_elemAbundances[eGlobal] += m_atoms(eGlobal, kGlobal) * spMoles;
460 double MultiPhase::volume()
const
463 for (
size_t i = 0; i < nPhases(); i++) {
464 double vol = 1.0/m_phase[i]->molarDensity();
465 sum += m_moles[i] * vol;
470 double MultiPhase::equilibrate_MultiPhaseEquil(
int XY,
double err,
int maxsteps,
471 int maxiter,
int loglevel)
482 return e.equilibrate(XY, err, maxsteps, loglevel);
483 }
else if (XY == HP) {
484 double h0 = enthalpy();
485 double Tlow = 0.5*m_Tmin;
486 double Thigh = 2.0*m_Tmax;
488 for (
int n = 0; n < maxiter; n++) {
496 e.equilibrate(TP, err, maxsteps, loglevel);
497 double hnow = enthalpy();
509 if (m_temp < Thigh) {
516 double cpb = (Hhigh - Hlow)/(Thigh - Tlow);
517 dt = (h0 - hnow)/cpb;
519 double dtmax = 0.5*fabs(Thigh - Tlow);
524 double tnew = sqrt(Tlow*Thigh);
528 double herr = fabs((h0 - hnow)/h0);
533 double tnew = m_temp + dt;
537 setTemperature(tnew);
549 double tnew = 0.5*(m_temp + Thigh);
550 if (fabs(tnew - m_temp) < 1.0) {
553 setTemperature(tnew);
557 throw CanteraError(
"MultiPhase::equilibrate_MultiPhaseEquil",
558 "No convergence for T");
559 }
else if (XY == SP) {
560 double s0 = entropy();
562 double Thigh = 1.0e6;
563 for (
int n = 0; n < maxiter; n++) {
567 e.equilibrate(TP, err, maxsteps, loglevel);
568 double snow = entropy();
570 Tlow = std::max(Tlow, m_temp);
572 Thigh = std::min(Thigh, m_temp);
574 double dt = (s0 - snow)*m_temp/cp();
575 double dtmax = 0.5*fabs(Thigh - Tlow);
576 dtmax = (dtmax > 500.0 ? 500.0 : dtmax);
584 double tnew = m_temp + dt;
585 setTemperature(tnew);
596 double tnew = 0.5*(m_temp + Thigh);
597 setTemperature(tnew);
601 throw CanteraError(
"MultiPhase::equilibrate_MultiPhaseEquil",
602 "No convergence for T");
603 }
else if (XY == TV) {
604 double v0 = volume();
606 for (
int n = 0; n < maxiter; n++) {
607 double pnow = pressure();
611 e.equilibrate(TP, err, maxsteps, loglevel);
612 double vnow = volume();
613 double verr = fabs((v0 - vnow)/v0);
619 setPressure(pnow*1.01);
620 double dVdP = (volume() - vnow)/(0.01*pnow);
621 setPressure(pnow + 0.5*(v0 - vnow)/dVdP);
624 throw CanteraError(
"MultiPhase::equilibrate_MultiPhaseEquil",
630 void MultiPhase::equilibrate(
const string& XY,
const string& solver,
631 double rtol,
int max_steps,
int max_iter,
632 int estimate_equil,
int log_level)
636 vector<double> initial_moleFractions = m_moleFractions;
637 vector<double> initial_moles = m_moles;
638 double initial_T = m_temp;
639 double initial_P = m_press;
641 if (solver ==
"auto" || solver ==
"vcs") {
643 debuglog(
"Trying VCS equilibrium solver\n", log_level);
645 int ret = eqsolve.
equilibrate(ixy, estimate_equil, log_level-1,
649 "VCS solver failed. Return code: {}", ret);
651 debuglog(
"VCS solver succeeded\n", log_level);
653 }
catch (std::exception& err) {
654 debuglog(
"VCS solver failed.\n", log_level);
656 m_moleFractions = initial_moleFractions;
657 m_moles = initial_moles;
661 if (solver ==
"auto") {
668 if (solver ==
"auto" || solver ==
"gibbs") {
670 debuglog(
"Trying MultiPhaseEquil (Gibbs) equilibrium solver\n",
672 equilibrate_MultiPhaseEquil(ixy, rtol, max_steps, max_iter,
674 debuglog(
"MultiPhaseEquil solver succeeded\n", log_level);
676 }
catch (std::exception& err) {
677 debuglog(
"MultiPhaseEquil solver failed.\n", log_level);
679 m_moleFractions = initial_moleFractions;
680 m_moles = initial_moles;
688 if (solver !=
"auto") {
690 "Invalid solver specified: '" + solver +
"'");
694 void MultiPhase::setTemperature(
const double T)
703 void MultiPhase::checkElementIndex(
size_t m)
const
706 throw IndexError(
"MultiPhase::checkElementIndex",
"elements", m, m_nel-1);
710 void MultiPhase::checkElementArraySize(
size_t mm)
const
713 throw ArraySizeError(
"MultiPhase::checkElementArraySize", mm, m_nel);
717 string MultiPhase::elementName(
size_t m)
const
722 size_t MultiPhase::elementIndex(
const string& name)
const
724 for (
size_t e = 0; e < m_nel; e++) {
725 if (m_enames[e] == name) {
732 void MultiPhase::checkSpeciesIndex(
size_t k)
const
735 throw IndexError(
"MultiPhase::checkSpeciesIndex",
"species", k, m_nsp-1);
739 void MultiPhase::checkSpeciesArraySize(
size_t kk)
const
742 throw ArraySizeError(
"MultiPhase::checkSpeciesArraySize", kk, m_nsp);
746 string MultiPhase::speciesName(
const size_t k)
const
751 double MultiPhase::nAtoms(
const size_t kGlob,
const size_t mGlob)
const
753 return m_atoms(mGlob, kGlob);
756 void MultiPhase::getMoleFractions(
double*
const x)
const
758 std::copy(m_moleFractions.begin(), m_moleFractions.end(), x);
761 string MultiPhase::phaseName(
const size_t iph)
const
763 return m_phase[iph]->name();
766 int MultiPhase::phaseIndex(
const string& pName)
const
768 for (
int iph = 0; iph < (int) nPhases(); iph++) {
769 if (m_phase[iph]->name() == pName) {
776 double MultiPhase::phaseMoles(
const size_t n)
const
781 void MultiPhase::setPhaseMoles(
const size_t n,
const double moles)
786 size_t MultiPhase::speciesPhaseIndex(
const size_t kGlob)
const
788 return m_spphase[kGlob];
791 double MultiPhase::moleFraction(
const size_t kGlob)
const
793 return m_moleFractions[kGlob];
796 bool MultiPhase::tempOK(
const size_t p)
const
801 void MultiPhase::uploadMoleFractionsFromPhases()
804 for (
size_t ip = 0; ip < nPhases(); ip++) {
809 calcElemAbundances();
812 void MultiPhase::updatePhases()
const
815 for (
size_t p = 0; p < nPhases(); p++) {
816 m_phase[p]->setState_TPX(m_temp, m_press, m_moleFractions.data() + loc);
817 loc += m_phase[p]->nSpecies();
819 if (m_temp < m_phase[p]->minTemp() || m_temp > m_phase[p]->maxTemp()) {
820 m_temp_OK[p] =
false;
828 for (
size_t ip = 0; ip < x.
nPhases(); ip++) {
830 s <<
"*************** " << x.
phase(ip).
name() <<
" *****************" << std::endl;
832 s <<
"*************** Phase " << ip <<
" *****************" << std::endl;
834 s <<
"Moles: " << x.
phaseMoles(ip) << std::endl;
Headers for the MultiPhase object that is used to set up multiphase equilibrium problems (see Chemica...
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
Base class for exceptions thrown by Cantera classes.
An array index is out of range.
Multiphase chemical equilibrium solver.
A class for multiphase mixtures.
vector< ThermoPhase * > m_phase
Vector of the ThermoPhase pointers.
size_t nPhases() const
Number of phases.
double phaseMoles(const size_t n) const
Return the number of moles in phase n.
void updatePhases() const
Set the states of the phase objects to the locally-stored state within this MultiPhase object.
ThermoPhase & phase(size_t n)
Return a reference to phase n.
vector< double > m_moles
Vector of the number of moles in each phase.
size_t nSpecies() const
Returns the number of species in the phase.
int atomicNumber(size_t m) const
Atomic number of element m.
double temperature() const
Temperature (K).
size_t elementIndex(const string &name) const
Return the index of element named 'name'.
string speciesName(size_t k) const
Name of the species with index k.
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
double nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
size_t nElements() const
Number of elements.
virtual double pressure() const
Return the thermodynamic pressure (Pa).
string elementName(size_t m) const
Name of the element with index m.
string name() const
Return the name of the phase.
Base class for a phase with thermodynamic properties.
virtual double minTemp(size_t k=npos) const
Minimum temperature for which the thermodynamic data for the species or phase are valid.
virtual void setState_TPX(double t, double p, const double *x)
Set the temperature (K), pressure (Pa), and mole fractions.
virtual string report(bool show_thermo=true, double threshold=-1e-14) const
returns a summary of the state of the phase as a string
virtual double maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
Cantera's Interface to the Multiphase chemical equilibrium solver.
int equilibrate(int XY, int estimateEquil=0, int printLvl=0, double err=1.0e-6, int maxsteps=VCS_MAXSTEPS, int loglevel=-99)
Equilibrate the solution using the current element abundances stored in the MultiPhase object.
Composition parseCompString(const string &ss, const vector< string > &names)
Parse a composition string into a map consisting of individual key:composition pairs.
virtual bool compatibleWithMultiPhase() const
Indicates whether this phase type can be used with class MultiPhase for equilibrium calculations.
void debuglog(const string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
const double Faraday
Faraday constant [C/kmol].
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
const double Undef
Fairly random number to be used to initialize variables against to see if they are subsequently defin...
const U & getValue(const map< T, U > &m, const T &key, const U &default_val)
Const accessor for a value in a map.
std::ostream & operator<<(std::ostream &s, MultiPhase &x)
Function to output a MultiPhase description to a stream.
int _equilflag(const char *xy)
map property strings to integers
map< string, double > Composition
Map from string names to doubles.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...
Interface class for the vcsnonlinear solver.