20 string Phase::name()
const
25 void Phase::setName(
const string& name)
30 size_t Phase::nElements()
const
35 void Phase::checkElementIndex(
size_t m)
const
38 throw IndexError(
"Phase::checkElementIndex",
"elements", m, m_mm-1);
42 void Phase::checkElementArraySize(
size_t mm)
const
49 string Phase::elementName(
size_t m)
const
52 return m_elementNames[m];
55 size_t Phase::elementIndex(
const string& elementName)
const
57 for (
size_t i = 0; i < m_mm; i++) {
58 if (m_elementNames[i] == elementName) {
67 return m_elementNames;
70 double Phase::atomicWeight(
size_t m)
const
72 return m_atomicWeights[m];
75 double Phase::entropyElement298(
size_t m)
const
78 return m_entropy298[m];
81 const vector<double>& Phase::atomicWeights()
const
83 return m_atomicWeights;
86 int Phase::atomicNumber(
size_t m)
const
88 return m_atomicNumbers[m];
91 int Phase::elementType(
size_t m)
const
93 return m_elem_type[m];
96 int Phase::changeElementType(
int m,
int elem_type)
98 int old = m_elem_type[m];
99 m_elem_type[m] = elem_type;
103 double Phase::nAtoms(
size_t k,
size_t m)
const
105 checkElementIndex(m);
106 checkSpeciesIndex(k);
107 return m_speciesComp[m_mm * k + m];
110 size_t Phase::findSpeciesLower(
const string& name)
const
114 auto it = m_speciesLower.find(nLower);
115 if (it == m_speciesLower.end()) {
121 "Lowercase species name '{}' is not unique. "
122 "Set Phase::caseSensitiveSpecies to true to "
123 "enforce case sensitive species names", nLower);
129 size_t Phase::speciesIndex(
const string& nameStr)
const
133 auto it = m_speciesIndices.find(nameStr);
134 if (it != m_speciesIndices.end()) {
136 }
else if (!m_caseSensitiveSpecies) {
137 loc = findSpeciesLower(nameStr);
142 string Phase::speciesName(
size_t k)
const
144 checkSpeciesIndex(k);
145 return m_speciesNames[k];
148 const vector<string>& Phase::speciesNames()
const
150 return m_speciesNames;
153 void Phase::checkSpeciesIndex(
size_t k)
const
156 throw IndexError(
"Phase::checkSpeciesIndex",
"species", k, m_kk-1);
160 void Phase::checkSpeciesArraySize(
size_t kk)
const
167 map<string, size_t> Phase::nativeState()
const
170 if (isCompressible()) {
171 return { {
"T", 0}, {
"D", 1} };
173 return { {
"T", 0}, {
"P", 1} };
176 if (isCompressible()) {
177 return { {
"T", 0}, {
"D", 1}, {
"Y", 2} };
179 return { {
"T", 0}, {
"P", 1}, {
"Y", 2} };
184 string Phase::nativeMode()
const
186 map<size_t, string> states;
187 for (
auto& [name, value] : nativeState()) {
188 states[value] = name;
191 for (
auto& [value, name] : states) {
197 vector<string> Phase::fullStates()
const
200 if (isCompressible()) {
201 return {
"TD",
"TP",
"UV",
"DP",
"HP",
"SP",
"SV"};
203 return {
"TP",
"HP",
"SP"};
206 if (isCompressible()) {
207 return {
"TDX",
"TDY",
"TPX",
"TPY",
"UVX",
"UVY",
"DPX",
"DPY",
208 "HPX",
"HPY",
"SPX",
"SPY",
"SVX",
"SVY"};
210 return {
"TPX",
"TPY",
"HPX",
"HPY",
"SPX",
"SPY"};
215 vector<string> Phase::partialStates()
const
220 if (isCompressible()) {
221 return {
"TD",
"TP",
"UV",
"DP",
"HP",
"SP",
"SV"};
223 return {
"TP",
"HP",
"SP"};
228 size_t Phase::stateSize()
const {
232 return nSpecies() + 2;
236 void Phase::saveState(vector<double>& state)
const
238 state.resize(stateSize());
239 saveState(state.size(), &state[0]);
242 void Phase::saveState(
size_t lenstate,
double* state)
const
244 auto native = nativeState();
247 state[native.at(
"T")] = temperature();
248 if (isCompressible()) {
249 state[native.at(
"D")] = density();
251 state[native.at(
"P")] = pressure();
253 if (native.count(
"X")) {
254 getMoleFractions(state + native[
"X"]);
255 }
else if (native.count(
"Y")) {
256 getMassFractions(state + native[
"Y"]);
260 void Phase::restoreState(
const vector<double>& state)
262 restoreState(state.size(),&state[0]);
265 void Phase::restoreState(
size_t lenstate,
const double* state)
267 size_t ls = stateSize();
273 auto native = nativeState();
274 setTemperature(state[native.at(
"T")]);
275 if (isCompressible()) {
276 setDensity(state[native.at(
"D")]);
278 setPressure(state[native.at(
"P")]);
281 if (native.count(
"X")) {
282 setMoleFractions_NoNorm(state + native[
"X"]);
283 }
else if (native.count(
"Y")) {
284 setMassFractions_NoNorm(state + native[
"Y"]);
286 compositionChanged();
289 void Phase::setMoleFractions(
const double*
const x)
295 for (
size_t k = 0; k < m_kk; k++) {
296 double xk = std::max(x[k], 0.0);
299 sum += m_molwts[k] * xk;
305 const double invSum = 1.0/sum;
306 for (
size_t k=0; k < m_kk; k++) {
307 m_ym[k] = m_y[k]*invSum;
312 for (
size_t k=0; k < m_kk; k++) {
313 m_y[k] = m_ym[k] * m_molwts[k];
318 compositionChanged();
321 void Phase::setMoleFractions_NoNorm(
const double*
const x)
323 m_mmw =
dot(x, x + m_kk, m_molwts.begin());
324 scale(x, x + m_kk, m_ym.begin(), 1.0/m_mmw);
325 transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(),
326 m_y.begin(), multiplies<double>());
327 compositionChanged();
332 vector<double> mf = getCompositionFromMap(xMap);
333 setMoleFractions(mf.data());
336 void Phase::setMoleFractionsByName(
const string& x)
341 void Phase::setMassFractions(
const double*
const y)
343 for (
size_t k = 0; k < m_kk; k++) {
344 m_y[k] = std::max(y[k], 0.0);
346 double norm = accumulate(m_y.begin(), m_y.end(), 0.0);
347 scale(m_y.begin(), m_y.end(), m_y.begin(), 1.0/norm);
349 transform(m_y.begin(), m_y.end(), m_rmolwts.begin(),
350 m_ym.begin(), multiplies<double>());
351 m_mmw = 1.0 / accumulate(m_ym.begin(), m_ym.end(), 0.0);
352 compositionChanged();
355 void Phase::setMassFractions_NoNorm(
const double*
const y)
358 copy(y, y + m_kk, m_y.begin());
359 transform(m_y.begin(), m_y.end(), m_rmolwts.begin(), m_ym.begin(),
360 multiplies<double>());
361 sum = accumulate(m_ym.begin(), m_ym.end(), 0.0);
363 compositionChanged();
368 vector<double> mf = getCompositionFromMap(yMap);
369 setMassFractions(mf.data());
372 void Phase::setMassFractionsByName(
const string& y)
377 void Phase::setState_TD(
double t,
double rho)
383 double Phase::molecularWeight(
size_t k)
const
385 checkSpeciesIndex(k);
389 void Phase::getMolecularWeights(
double* weights)
const
391 const vector<double>& mw = molecularWeights();
392 copy(mw.begin(), mw.end(), weights);
395 const vector<double>& Phase::molecularWeights()
const
400 const vector<double>& Phase::inverseMolecularWeights()
const
405 void Phase::getCharges(
double* charges)
const
407 copy(m_speciesCharge.begin(), m_speciesCharge.end(), charges);
413 for (
size_t k = 0; k < m_kk; k++) {
414 double x = moleFraction(k);
416 comp[speciesName(k)] = x;
425 for (
size_t k = 0; k < m_kk; k++) {
426 double x = massFraction(k);
428 comp[speciesName(k)] = x;
434 void Phase::getMoleFractions(
double*
const x)
const
436 scale(m_ym.begin(), m_ym.end(), x, m_mmw);
439 double Phase::moleFraction(
size_t k)
const
441 checkSpeciesIndex(k);
442 return m_ym[k] * m_mmw;
445 double Phase::moleFraction(
const string& nameSpec)
const
447 size_t iloc = speciesIndex(nameSpec);
449 return moleFraction(iloc);
455 double Phase::massFraction(
size_t k)
const
457 checkSpeciesIndex(k);
461 double Phase::massFraction(
const string& nameSpec)
const
463 size_t iloc = speciesIndex(nameSpec);
465 return massFractions()[iloc];
471 void Phase::getMassFractions(
double*
const y)
const
473 copy(m_y.begin(), m_y.end(), y);
476 double Phase::concentration(
const size_t k)
const
478 checkSpeciesIndex(k);
479 return m_y[k] * m_dens * m_rmolwts[k];
482 void Phase::getConcentrations(
double*
const c)
const
484 scale(m_ym.begin(), m_ym.end(), c, m_dens);
487 void Phase::setConcentrations(
const double*
const conc)
489 assertCompressible(
"setConcentrations");
492 double sum = 0.0, norm = 0.0;
493 for (
size_t k = 0; k != m_kk; ++k) {
494 double ck = std::max(conc[k], 0.0);
496 sum += ck * m_molwts[k];
501 double rsum = 1.0/sum;
502 for (
size_t k = 0; k != m_kk; ++k) {
503 m_ym[k] = m_y[k] * rsum;
504 m_y[k] = m_ym[k] * m_molwts[k];
506 compositionChanged();
509 void Phase::setConcentrationsNoNorm(
const double*
const conc)
511 assertCompressible(
"setConcentrationsNoNorm");
513 double sum = 0.0, norm = 0.0;
514 for (
size_t k = 0; k != m_kk; ++k) {
515 sum += conc[k] * m_molwts[k];
520 double rsum = 1.0/sum;
521 for (
size_t k = 0; k != m_kk; ++k) {
522 m_ym[k] = conc[k] * rsum;
523 m_y[k] = m_ym[k] * m_molwts[k];
525 compositionChanged();
528 void Phase::setMolesNoTruncate(
const double*
const N)
531 copy(N, N + m_kk, m_ym.begin());
532 double totalMoles = accumulate(m_ym.begin(), m_ym.end(), 0.0);
534 copy(N, N + m_kk, m_y.begin());
535 transform(m_y.begin(), m_y.end(), m_molwts.begin(), m_y.begin(), multiplies<double>());
536 double totalMass = accumulate(m_y.begin(), m_y.end(), 0.0);
538 m_mmw = totalMass/totalMoles;
540 scale(m_y.begin(), m_y.end(), m_y.begin(), 1/totalMass);
542 scale(m_ym.begin(), m_ym.end(), m_ym.begin(), 1/(m_mmw * totalMoles));
544 compositionChanged();
547 double Phase::elementalMassFraction(
const size_t m)
const
549 checkElementIndex(m);
551 for (
size_t k = 0; k != m_kk; ++k) {
552 Z_m += nAtoms(k, m) * atomicWeight(m) / molecularWeight(k)
558 double Phase::elementalMoleFraction(
const size_t m)
const
560 checkElementIndex(m);
562 for (
size_t k = 0; k < m_kk; k++) {
564 for (
size_t j = 0; j < nElements(); j++) {
565 atoms += nAtoms(k, j);
567 denom += atoms * moleFraction(k);
569 double numerator = 0.0;
570 for (
size_t k = 0; k != m_kk; ++k) {
571 numerator += nAtoms(k, m) * moleFraction(k);
573 return numerator / denom;
576 double Phase::molarDensity()
const
578 return density()/meanMolecularWeight();
581 double Phase::molarVolume()
const
583 return 1.0/molarDensity();
586 void Phase::setDensity(
const double density_)
588 assertCompressible(
"setDensity");
589 if (density_ > 0.0) {
593 "density must be positive. density = {}", density_);
597 void Phase::assignDensity(
const double density_)
599 if (density_ > 0.0) {
603 "density must be positive. density = {}", density_);
607 double Phase::chargeDensity()
const
610 for (
size_t k = 0; k < m_kk; k++) {
611 cdens += charge(k)*moleFraction(k);
616 double Phase::mean_X(
const double*
const Q)
const
618 return m_mmw*std::inner_product(m_ym.begin(), m_ym.end(), Q, 0.0);
621 double Phase::mean_X(
const vector<double>& Q)
const
623 return m_mmw*std::inner_product(m_ym.begin(), m_ym.end(), Q.begin(), 0.0);
626 double Phase::sum_xlogx()
const
629 for (
size_t k = 0; k < m_kk; k++) {
630 sumxlogx += m_ym[k] * std::log(std::max(m_ym[k],
SmallNumber));
632 return m_mmw * sumxlogx + std::log(m_mmw);
635 size_t Phase::addElement(
const string& symbol,
double weight,
int atomic_number,
636 double entropy298,
int elem_type)
645 }
else if (weight == -12345.0) {
652 const static AnyMap db = AnyMap::fromYamlFile(
653 "element-standard-entropies.yaml");
654 const AnyMap& elem = db[
"elements"].getMapWhere(
"symbol", symbol);
661 auto iter = find(m_elementNames.begin(), m_elementNames.end(), symbol);
662 if (iter != m_elementNames.end()) {
663 size_t m = iter - m_elementNames.begin();
664 if (m_atomicWeights[m] != weight) {
666 "Duplicate elements ({}) have different weights", symbol);
674 m_atomicWeights.push_back(weight);
675 m_elementNames.push_back(symbol);
676 m_atomicNumbers.push_back(atomic_number);
677 m_entropy298.push_back(entropy298);
681 m_elem_type.push_back(elem_type);
687 vector<double> old(m_speciesComp);
688 m_speciesComp.resize(m_kk*m_mm, 0.0);
689 for (
size_t k = 0; k < m_kk; k++) {
690 size_t m_old = m_mm - 1;
691 for (
size_t m = 0; m < m_old; m++) {
692 m_speciesComp[k * m_mm + m] = old[k * (m_old) + m];
694 m_speciesComp[k * (m_mm) + (m_mm-1)] = 0.0;
701 bool Phase::addSpecies(shared_ptr<Species> spec) {
702 if (m_species.find(spec->name) != m_species.end()) {
704 "Phase '{}' already contains a species named '{}'.",
708 vector<double> comp(nElements());
709 for (
const auto& [eName, stoich] : spec->composition) {
710 size_t m = elementIndex(eName);
712 switch (m_undefinedElementBehavior) {
713 case UndefElement::ignore:
716 case UndefElement::add:
718 comp.resize(nElements());
719 m = elementIndex(eName);
722 case UndefElement::error:
725 "Species '{}' contains an undefined element '{}'.",
732 size_t ne = nElements();
733 const vector<double>& aw = atomicWeights();
734 if (spec->charge != 0.0) {
735 size_t eindex = elementIndex(
"E");
736 if (eindex !=
npos) {
737 double ecomp = comp[eindex];
738 if (fabs(spec->charge + ecomp) > 0.001) {
741 "Input charge and element E compositions differ "
742 "for species " + spec->name);
746 comp[eindex] = -spec->charge;
752 eindex = elementIndex(
"E");
754 comp[ne - 1] = - spec->charge;
759 for (
size_t m = 0; m < ne; m++) {
760 wt += comp[m] * aw[m];
766 wt = std::max(wt,
Tiny);
768 spec->setMolecularWeight(wt);
770 m_molwts.push_back(wt);
771 m_rmolwts.push_back(1.0/wt);
773 for (
size_t m = 0; m < ne; m++) {
774 m_speciesComp.push_back(comp[m]);
777 m_speciesNames.push_back(spec->name);
778 m_species[spec->name] = spec;
779 m_speciesIndices[spec->name] = m_kk;
780 m_speciesCharge.push_back(spec->charge);
783 if (m_speciesLower.find(nLower) == m_speciesLower.end()) {
784 m_speciesLower[nLower] = m_kk;
786 m_speciesLower[nLower] =
npos;
795 m_ym.push_back(m_rmolwts[0]);
796 m_mmw = 1.0 / m_ym[0];
805 void Phase::modifySpecies(
size_t k, shared_ptr<Species> spec)
807 if (speciesName(k) != spec->name) {
809 "New species name '{}' does not match existing name '{}'",
810 spec->name, speciesName(k));
812 const shared_ptr<Species>& old = m_species[spec->name];
813 if (spec->composition != old->composition) {
815 "New composition for '{}' does not match existing composition",
818 m_species[spec->name] = spec;
822 void Phase::addSpeciesAlias(
const string& name,
const string& alias)
824 if (speciesIndex(alias) !=
npos) {
826 "Invalid alias '{}': species already exists", alias);
828 size_t k = speciesIndex(name);
830 m_speciesIndices[alias] = k;
833 "Unable to add alias '{}' "
834 "(original species '{}' not found).", alias, name);
838 vector<string> Phase::findIsomers(
const Composition& compMap)
const
840 vector<string> isomerNames;
842 for (
const auto& [name, species] : m_species) {
843 if (species->composition == compMap) {
844 isomerNames.emplace_back(name);
851 vector<string> Phase::findIsomers(
const string& comp)
const
856 shared_ptr<Species> Phase::species(
const string& name)
const
858 size_t k = speciesIndex(name);
860 return m_species.at(speciesName(k));
863 "Unknown species '{}'", name);
867 shared_ptr<Species> Phase::species(
size_t k)
const
869 checkSpeciesIndex(k);
870 return m_species.at(m_speciesNames[k]);
873 void Phase::ignoreUndefinedElements() {
874 m_undefinedElementBehavior = UndefElement::ignore;
877 void Phase::addUndefinedElements() {
878 m_undefinedElementBehavior = UndefElement::add;
881 void Phase::throwUndefinedElements() {
882 m_undefinedElementBehavior = UndefElement::error;
885 bool Phase::ready()
const
890 void Phase::invalidateCache() {
895 void Phase::setMolecularWeight(
const int k,
const double mw)
898 m_rmolwts[k] = 1.0/mw;
900 transform(m_y.begin(), m_y.end(), m_rmolwts.begin(), m_ym.begin(),
901 multiplies<double>());
902 m_mmw = 1.0 / accumulate(m_ym.begin(), m_ym.end(), 0.0);
905 void Phase::compositionChanged() {
909 vector<double> Phase::getCompositionFromMap(
const Composition& comp)
const
911 vector<double> X(m_kk);
912 for (
const auto& [name, value] : comp) {
913 size_t loc = speciesIndex(name);
916 "Unknown species '{}'", name);
923 void Phase::massFractionsToMoleFractions(
const double* Y,
double* X)
const
926 for (
size_t k = 0; k != m_kk; ++k) {
927 rmmw += Y[k]/m_molwts[k];
930 throw CanteraError(
"Phase::massFractionsToMoleFractions",
931 "no input composition given");
933 for (
size_t k = 0; k != m_kk; ++k) {
934 X[k] = Y[k]/(rmmw*m_molwts[k]);
938 void Phase::moleFractionsToMassFractions(
const double* X,
double* Y)
const
940 double mmw =
dot(X, X+m_kk, m_molwts.data());
942 throw CanteraError(
"Phase::moleFractionsToMassFractions",
943 "no input composition given");
945 double rmmw = 1.0/mmw;
946 for (
size_t k = 0; k != m_kk; ++k) {
947 Y[k] = X[k]*m_molwts[k]*rmmw;
#define CT_ELEM_TYPE_ELECTRONCHARGE
This refers to conservation of electrons.
#define ENTROPY298_UNKNOWN
Number indicating we don't know the entropy of the element in its most stable state at 298....
Header file for class Phase.
Declaration for class Cantera::Species.
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
A map of string keys to values whose type can vary at runtime.
double convert(const string &key, const string &units) const
Convert the item stored by the given key to the units specified in units.
Base class for exceptions thrown by Cantera classes.
An array index is out of range.
string toLowerCopy(const string &input)
Convert to lower case.
Composition parseCompString(const string &ss, const vector< string > &names)
Parse a composition string into a map consisting of individual key:composition pairs.
double dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
const double Faraday
Faraday constant [C/kmol].
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
double getElementWeight(const string &ename)
Get the atomic weight of an element.
const double Tiny
Small number to compare differences of mole fractions against.
const vector< string > & elementNames()
Get a vector of the names of the elements defined in Cantera.
const double SmallNumber
smallest number to compare to zero.
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...