25 for (
size_t n = 0; n < mix.
nPhases(); n++) {
31 const vector<double>& phaseMoles)
33 for (
size_t n = 0; n < phases.size(); 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());
65 for (
size_t m = 0; m < p->
nElements(); m++) {
77 if (ename ==
"E" || ename ==
"e") {
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) {
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++) {
198 for (
size_t i = 0; i <
nPhases(); i++) {
222 double phasesum = 0.0;
223 size_t nsp =
m_phase[p]->nSpecies();
224 for (
size_t ik = 0; ik < nsp; ik++) {
235 for (
size_t i = 0; i <
nPhases(); i++) {
236 m_phase[i]->getChemPotentials(mu + loc);
246 for (
size_t i = 0; i <
nPhases(); i++) {
249 m_phase[i]->getChemPotentials(mu + loc);
251 m_phase[i]->getStandardChemPotentials(mu + loc);
273 for (
size_t i = 0; i <
nPhases(); i++) {
285 for (
size_t i = 0; i <
nPhases(); i++) {
297 for (
size_t i = 0; i <
nPhases(); i++) {
309 for (
size_t i = 0; i <
nPhases(); i++) {
321 for (
size_t i = 0; i <
nPhases(); i++) {
337 for (
size_t k = 0; k < p->
nSpecies(); k++) {
345 vector<double> moles(kk, 0.0);
346 for (
size_t k = 0; k < kk; k++) {
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;
381 for (
size_t ip = 0; ip <
nPhases(); ip++) {
384 double phasemoles = 0.0;
385 for (
size_t ik = 0; ik < nsp; ik++) {
391 if (phasemoles > 0.0) {
406 vector<double> tmpMoles(
m_nsp, 0.0);
408 tmpMoles[indexS] += addedMoles;
409 tmpMoles[indexS] = std::max(tmpMoles[indexS], 0.0);
433 for (
size_t eGlobal = 0; eGlobal <
m_nel; eGlobal++) {
442 for (
size_t eGlobal = 0; eGlobal <
m_nel; eGlobal++) {
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;
452 for (
size_t eGlobal = 0; eGlobal <
m_nel; eGlobal++) {
463 for (
size_t i = 0; i <
nPhases(); i++) {
464 double vol = 1.0/
m_phase[i]->molarDensity();
471 int maxiter,
int loglevel)
482 return e.equilibrate(XY, err, maxsteps, loglevel);
483 }
else if (XY == HP) {
486 double Thigh = 2.0*
m_Tmax;
488 for (
int n = 0; n < maxiter; n++) {
496 e.equilibrate(TP, err, maxsteps, loglevel);
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;
549 double tnew = 0.5*(
m_temp + Thigh);
550 if (fabs(tnew -
m_temp) < 1.0) {
557 throw CanteraError(
"MultiPhase::equilibrate_MultiPhaseEquil",
558 "No convergence for T");
559 }
else if (XY == SP) {
562 double Thigh = 1.0e6;
563 for (
int n = 0; n < maxiter; n++) {
567 e.equilibrate(TP, err, maxsteps, loglevel);
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;
596 double tnew = 0.5*(
m_temp + Thigh);
601 throw CanteraError(
"MultiPhase::equilibrate_MultiPhaseEquil",
602 "No convergence for T");
603 }
else if (XY == TV) {
606 for (
int n = 0; n < maxiter; n++) {
611 e.equilibrate(TP, err, maxsteps, loglevel);
613 double verr = fabs((v0 - vnow)/v0);
620 double dVdP = (
volume() - vnow)/(0.01*pnow);
624 throw CanteraError(
"MultiPhase::equilibrate_MultiPhaseEquil",
631 double rtol,
int max_steps,
int max_iter,
632 int estimate_equil,
int log_level)
637 vector<double> initial_moles =
m_moles;
638 double initial_T =
m_temp;
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);
661 if (solver ==
"auto") {
668 if (solver ==
"auto" || solver ==
"gibbs") {
670 debuglog(
"Trying MultiPhaseEquil (Gibbs) equilibrium solver\n",
674 debuglog(
"MultiPhaseEquil solver succeeded\n", log_level);
676 }
catch (std::exception& err) {
677 debuglog(
"MultiPhaseEquil solver failed.\n", log_level);
688 if (solver !=
"auto") {
690 "Invalid solver specified: '" + solver +
"'");
706 throw IndexError(
"MultiPhase::checkElementIndex",
"elements", m,
m_nel-1);
724 for (
size_t e = 0; e <
m_nel; e++) {
735 throw IndexError(
"MultiPhase::checkSpeciesIndex",
"species", k,
m_nsp-1);
768 for (
int iph = 0; iph < (int)
nPhases(); iph++) {
769 if (
m_phase[iph]->name() == pName) {
804 for (
size_t ip = 0; ip <
nPhases(); ip++) {
815 for (
size_t p = 0; p <
nPhases(); p++) {
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.
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
An array index is out of range.
Multiphase chemical equilibrium solver.
A class for multiphase mixtures.
void init()
Process phases and build atomic composition array.
void checkPhaseIndex(size_t m) const
Check that the specified phase index is in range Throws an exception if m is greater than nPhases()
size_t speciesIndex(size_t k, size_t p) const
Return the global index of the species belonging to phase number p with local index k within the phas...
bool solutionSpecies(size_t kGlob) const
Return true is species kGlob is a species in a multicomponent solution phase.
vector< ThermoPhase * > m_phase
Vector of the ThermoPhase pointers.
double nAtoms(const size_t kGlob, const size_t mGlob) const
Returns the Number of atoms of global element mGlob in global species kGlob.
void setMolesByName(const Composition &xMap)
Set the number of moles of species in the mixture.
void setMoles(const double *n)
Sets all of the global species mole numbers.
DenseMatrix m_atoms
Global Stoichiometric Coefficient array.
double gibbs() const
The Gibbs function of the mixture [J].
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range.
size_t m_nel
Number of distinct elements in all of the phases.
double speciesMoles(size_t kGlob) const
Returns the moles of global species k. units = kmol.
void getValidChemPotentials(double not_mu, double *mu, bool standard=false) const
Returns a vector of Valid chemical potentials.
double m_temp
Current value of the temperature (kelvin)
void calcElemAbundances() const
Calculate the element abundance vector.
size_t nSpecies() const
Number of species, summed over all phases.
double enthalpy() const
The enthalpy of the mixture [J].
double pressure() const
Pressure [Pa].
vector< size_t > m_spstart
Vector of ints containing of first species index in the global list of species for each phase.
vector< size_t > m_spphase
Mapping between the global species number and the phase ID.
void getMoles(double *molNum) const
Get the mole numbers of all species in the multiphase object.
double minTemp() const
Minimum temperature for which all solution phases have valid thermo data.
vector< double > m_moleFractions
Locally stored vector of mole fractions of all species comprising the MultiPhase object.
vector< double > m_elemAbundances
Vector of element abundances.
double equilibrate_MultiPhaseEquil(int XY, double err, int maxsteps, int maxiter, int loglevel)
Set the mixture to a state of chemical equilibrium using the MultiPhaseEquil solver.
void checkPhaseArraySize(size_t mm) const
Check that an array size is at least nPhases() Throws an exception if mm is less than nPhases().
vector< bool > m_temp_OK
Vector of bools indicating whether temperatures are ok for phases.
double phaseCharge(size_t p) const
Charge (Coulombs) of phase with index p.
int phaseIndex(const string &pName) const
Returns the index, given the phase name.
size_t nPhases() const
Number of phases.
size_t m_eloc
Global ID of the element corresponding to the electronic charge.
double entropy() const
The entropy of the mixture [J/K].
void getChemPotentials(double *mu) const
Returns a vector of Chemical potentials.
double moleFraction(const size_t kGlob) const
Returns the mole fraction of global species k.
double m_press
Current value of the pressure (Pa)
bool tempOK(size_t p) const
Return true if the phase p has valid thermo data for the current temperature.
void addPhases(vector< ThermoPhase * > &phases, const vector< double > &phaseMoles)
Add a vector of phases to the mixture.
void addPhase(ThermoPhase *p, double moles)
Add a phase to the mixture.
map< string, size_t > m_enamemap
Returns the global element index, given the element string name.
void addSpeciesMoles(const int indexS, const double addedMoles)
Adds moles of a certain species to the mixture.
size_t elementIndex(const string &name) const
Returns the index of the element with name name.
void checkElementArraySize(size_t mm) const
Check that an array size is at least nElements().
size_t speciesPhaseIndex(const size_t kGlob) const
Returns the phase index of the Kth "global" species.
void setPressure(double P)
Set the pressure [Pa].
void setState_TPMoles(const double T, const double Pres, const double *Moles)
Set the state of the underlying ThermoPhase objects in one call.
vector< string > m_enames
String names of the global elements.
string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
double phaseMoles(const size_t n) const
Return the number of moles in phase n.
void getMoleFractions(double *const x) const
Returns the global Species mole fractions.
size_t m_nsp
Number of distinct species in all of the phases.
void setPhaseMoleFractions(const size_t n, const double *const x)
Set the Mole fractions of the nth phase.
vector< int > m_atomicNumber
Atomic number of each global element.
double volume() const
The total mixture volume [m^3].
void uploadMoleFractionsFromPhases()
Update the locally-stored composition within this object to match the current compositions of the pha...
bool m_init
True if the init() routine has been called, and the MultiPhase frozen.
vector< string > m_snames
Vector of species names in the problem.
double m_Tmin
Minimum temperature for which thermo parameterizations are valid.
void updatePhases() const
Set the states of the phase objects to the locally-stored state within this MultiPhase object.
void getElemAbundances(double *elemAbundances) const
Retrieves a vector of element abundances.
void checkSpeciesArraySize(size_t kk) const
Check that an array size is at least nSpecies().
double charge() const
Total charge summed over all phases (Coulombs).
string phaseName(const size_t iph) const
Returns the name of the n'th phase.
double cp() const
Heat capacity at constant pressure [J/K].
void setTemperature(const double T)
Set the temperature [K].
void setState_TP(const double T, const double Pres)
Set the state of the underlying ThermoPhase objects in one call.
void checkElementIndex(size_t m) const
Check that the specified element index is in range.
ThermoPhase & phase(size_t n)
Return a reference to phase n.
double IntEnergy() const
The internal energy of the mixture [J].
string elementName(size_t m) const
Returns the name of the global element m.
vector< double > m_moles
Vector of the number of moles in each phase.
void setPhaseMoles(const size_t n, const double moles)
Set the number of moles of phase with index n.
double elementMoles(size_t m) const
Total moles of global element m, summed over all phases.
double m_Tmax
Minimum temperature for which thermo parameterizations are valid.
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.
void equilibrate(const string &XY, const string &solver="auto", double rtol=1e-9, int max_steps=50000, int max_iter=100, int estimate_equil=0, int log_level=0)
Equilibrate a MultiPhase object.
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...
int _equilflag(const char *xy)
map property strings to integers
std::ostream & operator<<(std::ostream &s, const Array2D &m)
Output the current contents of the Array2D object.
const U & getValue(const map< T, U > &m, const T &key, const U &default_val)
Const accessor for a value in a map.
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.