Cantera  3.2.0a4
Loading...
Searching...
No Matches
Phase.cpp
Go to the documentation of this file.
1/**
2 * @file Phase.cpp
3 * Definition file for class Phase.
4 */
5
6// This file is part of Cantera. See License.txt in the top-level directory or
7// at https://cantera.org/license.txt for license and copyright information.
8
14
15using namespace std;
16
17namespace Cantera
18{
19
20string Phase::name() const
21{
22 return m_name;
23}
24
25void Phase::setName(const string& name)
26{
27 m_name = name;
28}
29
30size_t Phase::nElements() const
31{
32 return m_mm;
33}
34
35void Phase::checkElementIndex(size_t m) const
36{
37 if (m >= m_mm) {
38 throw IndexError("Phase::checkElementIndex", "elements", m, m_mm);
39 }
40}
41
42void Phase::checkElementArraySize(size_t mm) const
43{
44 if (m_mm > mm) {
45 throw ArraySizeError("Phase::checkElementArraySize", mm, m_mm);
46 }
47}
48
49string Phase::elementName(size_t m) const
50{
52 return m_elementNames[m];
53}
54
55size_t Phase::elementIndex(const string& elementName) const
56{
57 for (size_t i = 0; i < m_mm; i++) {
58 if (m_elementNames[i] == elementName) {
59 return i;
60 }
61 }
62 return npos;
63}
64
65const vector<string>& Phase::elementNames() const
66{
67 return m_elementNames;
68}
69
70double Phase::atomicWeight(size_t m) const
71{
72 return m_atomicWeights[m];
73}
74
75double Phase::entropyElement298(size_t m) const
76{
78 return m_entropy298[m];
79}
80
81const vector<double>& Phase::atomicWeights() const
82{
83 return m_atomicWeights;
84}
85
86int Phase::atomicNumber(size_t m) const
87{
88 return m_atomicNumbers[m];
89}
90
91int Phase::elementType(size_t m) const
92{
93 return m_elem_type[m];
94}
95
96int Phase::changeElementType(int m, int elem_type)
97{
98 int old = m_elem_type[m];
99 m_elem_type[m] = elem_type;
100 return old;
101}
102
103double Phase::nAtoms(size_t k, size_t m) const
104{
107 return m_speciesComp[m_mm * k + m];
108}
109
110size_t Phase::findSpeciesLower(const string& name) const
111{
112 string nLower = toLowerCopy(name);
113 size_t loc = npos;
114 auto it = m_speciesLower.find(nLower);
115 if (it == m_speciesLower.end()) {
116 return npos;
117 } else {
118 loc = it->second;
119 if (loc == npos) {
120 throw CanteraError("Phase::findSpeciesLower",
121 "Lowercase species name '{}' is not unique. "
122 "Set Phase::caseSensitiveSpecies to true to "
123 "enforce case sensitive species names", nLower);
124 }
125 }
126 return loc;
127}
128
129size_t Phase::speciesIndex(const string& nameStr) const
130{
131 size_t loc = npos;
132
133 auto it = m_speciesIndices.find(nameStr);
134 if (it != m_speciesIndices.end()) {
135 return it->second;
136 } else if (!m_caseSensitiveSpecies) {
137 loc = findSpeciesLower(nameStr);
138 }
139 return loc;
140}
141
142string Phase::speciesName(size_t k) const
143{
145 return m_speciesNames[k];
146}
147
148const vector<string>& Phase::speciesNames() const
149{
150 return m_speciesNames;
151}
152
153void Phase::checkSpeciesIndex(size_t k) const
154{
155 if (k >= m_kk) {
156 throw IndexError("Phase::checkSpeciesIndex", "species", k, m_kk);
157 }
158}
159
160void Phase::checkSpeciesArraySize(size_t kk) const
161{
162 if (m_kk > kk) {
163 throw ArraySizeError("Phase::checkSpeciesArraySize", kk, m_kk);
164 }
165}
166
167map<string, size_t> Phase::nativeState() const
168{
169 if (isPure()) {
170 if (isCompressible()) {
171 return { {"T", 0}, {"D", 1} };
172 } else {
173 return { {"T", 0}, {"P", 1} };
174 }
175 } else {
176 if (isCompressible()) {
177 return { {"T", 0}, {"D", 1}, {"Y", 2} };
178 } else {
179 return { {"T", 0}, {"P", 1}, {"Y", 2} };
180 }
181 }
182}
183
184string Phase::nativeMode() const
185{
186 map<size_t, string> states; // reverse lookup
187 for (auto& [name, value] : nativeState()) {
188 states[value] = name;
189 }
190 string out;
191 for (auto& [value, name] : states) {
192 out += name;
193 }
194 return out;
195}
196
197vector<string> Phase::fullStates() const
198{
199 if (isPure()) {
200 if (isCompressible()) {
201 return {"TD", "TP", "UV", "DP", "HP", "SP", "SV"};
202 } else {
203 return {"TP", "HP", "SP"};
204 }
205 } else {
206 if (isCompressible()) {
207 return {"TDX", "TDY", "TPX", "TPY", "UVX", "UVY", "DPX", "DPY",
208 "HPX", "HPY", "SPX", "SPY", "SVX", "SVY"};
209 } else {
210 return {"TPX", "TPY", "HPX", "HPY", "SPX", "SPY"};
211 }
212 }
213}
214
215vector<string> Phase::partialStates() const
216{
217 if (isPure()) {
218 return {};
219 } else {
220 if (isCompressible()) {
221 return {"TD", "TP", "UV", "DP", "HP", "SP", "SV"};
222 } else {
223 return {"TP", "HP", "SP"};
224 }
225 }
226}
227
228void Phase::savePartialState(size_t lenstate, double* state) const
229{
230 if (lenstate < partialStateSize()) {
231 throw ArraySizeError("Phase::savePartialState", lenstate, partialStateSize());
232 }
233
234 auto native = nativeState();
235 state[native.at("T")] = temperature();
236 if (isCompressible()) {
237 state[native.at("D")] = density();
238 } else {
239 state[native.at("P")] = pressure();
240 }
241}
242
243void Phase::restorePartialState(size_t lenstate, const double* state)
244{
245 if (lenstate < partialStateSize()) {
246 throw ArraySizeError("Phase::restorePartialState", lenstate, partialStateSize());
247 }
248
249 auto native = nativeState();
250 setTemperature(state[native.at("T")]);
251 if (isCompressible()) {
252 setDensity(state[native.at("D")]);
253 } else {
254 setPressure(state[native.at("P")]);
255 }
256}
257
258size_t Phase::stateSize() const {
259 if (isPure()) {
260 return partialStateSize();
261 } else {
262 return partialStateSize() + nSpecies();
263 }
264}
265
266void Phase::saveState(vector<double>& state) const
267{
268 state.resize(stateSize());
269 saveState(state.size(), &state[0]);
270}
271
272void Phase::saveState(size_t lenstate, double* state) const
273{
274 if (lenstate < stateSize()) {
275 throw ArraySizeError("Phase::saveState", lenstate, stateSize());
276 }
277 savePartialState(lenstate, state);
278 auto native = nativeState();
279 if (native.count("X")) {
280 getMoleFractions(state + native["X"]);
281 } else if (native.count("Y")) {
282 getMassFractions(state + native["Y"]);
283 }
284}
285
286void Phase::restoreState(const vector<double>& state)
287{
288 restoreState(state.size(),&state[0]);
289}
290
291void Phase::restoreState(size_t lenstate, const double* state)
292{
293 size_t ls = stateSize();
294 if (lenstate < ls) {
295 throw ArraySizeError("Phase::restoreState",
296 lenstate, ls);
297 }
298 restorePartialState(lenstate, state);
299
300 auto native = nativeState();
301
302 if (native.count("X")) {
303 setMoleFractions_NoNorm(state + native["X"]);
304 } else if (native.count("Y")) {
305 setMassFractions_NoNorm(state + native["Y"]);
306 }
308}
309
310void Phase::setMoleFractions(const double* const x)
311{
312 // Use m_y as a temporary work vector for the non-negative mole fractions
313 double norm = 0.0;
314 // sum is calculated below as the unnormalized molecular weight
315 double sum = 0;
316 for (size_t k = 0; k < m_kk; k++) {
317 double xk = std::max(x[k], 0.0); // Ignore negative mole fractions
318 m_y[k] = xk;
319 norm += xk;
320 sum += m_molwts[k] * xk;
321 }
322
323 // Set m_ym to the normalized mole fractions divided by the normalized mean
324 // molecular weight:
325 // m_ym_k = X_k / (sum_k X_k M_k)
326 const double invSum = 1.0/sum;
327 for (size_t k=0; k < m_kk; k++) {
328 m_ym[k] = m_y[k]*invSum;
329 }
330
331 // Now set m_y to the normalized mass fractions:
332 // m_y = X_k M_k / (sum_k X_k M_k)
333 for (size_t k=0; k < m_kk; k++) {
334 m_y[k] = m_ym[k] * m_molwts[k];
335 }
336
337 // Calculate the normalized molecular weight
338 m_mmw = sum/norm;
340}
341
342void Phase::setMoleFractions_NoNorm(const double* const x)
343{
344 m_mmw = dot(x, x + m_kk, m_molwts.begin());
345 scale(x, x + m_kk, m_ym.begin(), 1.0/m_mmw);
346 transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(),
347 m_y.begin(), multiplies<double>());
349}
350
352{
353 vector<double> mf = getCompositionFromMap(xMap);
354 setMoleFractions(mf.data());
355}
356
357void Phase::setMoleFractionsByName(const string& x)
358{
360}
361
362void Phase::setMassFractions(const double* const y)
363{
364 for (size_t k = 0; k < m_kk; k++) {
365 m_y[k] = std::max(y[k], 0.0); // Ignore negative mass fractions
366 }
367 double norm = accumulate(m_y.begin(), m_y.end(), 0.0);
368 scale(m_y.begin(), m_y.end(), m_y.begin(), 1.0/norm);
369
370 transform(m_y.begin(), m_y.end(), m_rmolwts.begin(),
371 m_ym.begin(), multiplies<double>());
372 m_mmw = 1.0 / accumulate(m_ym.begin(), m_ym.end(), 0.0);
374}
375
376void Phase::setMassFractions_NoNorm(const double* const y)
377{
378 double sum = 0.0;
379 copy(y, y + m_kk, m_y.begin());
380 transform(m_y.begin(), m_y.end(), m_rmolwts.begin(), m_ym.begin(),
381 multiplies<double>());
382 sum = accumulate(m_ym.begin(), m_ym.end(), 0.0);
383 m_mmw = 1.0/sum;
385}
386
388{
389 vector<double> mf = getCompositionFromMap(yMap);
390 setMassFractions(mf.data());
391}
392
393void Phase::setMassFractionsByName(const string& y)
394{
396}
397
398void Phase::setState_TD(double t, double rho)
399{
400 vector<double> state(partialStateSize());
401 savePartialState(state.size(), state.data());
402 try {
404 setDensity(rho);
405 } catch (std::exception&) {
406 restorePartialState(state.size(), state.data());
407 throw;
408 }
409}
410
411double Phase::molecularWeight(size_t k) const
412{
414 return m_molwts[k];
415}
416
417void Phase::getMolecularWeights(double* weights) const
418{
419 const vector<double>& mw = molecularWeights();
420 copy(mw.begin(), mw.end(), weights);
421}
422
423const vector<double>& Phase::molecularWeights() const
424{
425 return m_molwts;
426}
427
428const vector<double>& Phase::inverseMolecularWeights() const
429{
430 return m_rmolwts;
431}
432
433void Phase::getCharges(double* charges) const
434{
435 copy(m_speciesCharge.begin(), m_speciesCharge.end(), charges);
436}
437
439{
440 Composition comp;
441 for (size_t k = 0; k < m_kk; k++) {
442 double x = moleFraction(k);
443 if (x > threshold) {
444 comp[speciesName(k)] = x;
445 }
446 }
447 return comp;
448}
449
451{
452 Composition comp;
453 for (size_t k = 0; k < m_kk; k++) {
454 double x = massFraction(k);
455 if (x > threshold) {
456 comp[speciesName(k)] = x;
457 }
458 }
459 return comp;
460}
461
462void Phase::getMoleFractions(double* const x) const
463{
464 scale(m_ym.begin(), m_ym.end(), x, m_mmw);
465}
466
467double Phase::moleFraction(size_t k) const
468{
470 return m_ym[k] * m_mmw;
471}
472
473double Phase::moleFraction(const string& nameSpec) const
474{
475 size_t iloc = speciesIndex(nameSpec);
476 if (iloc != npos) {
477 return moleFraction(iloc);
478 } else {
479 return 0.0;
480 }
481}
482
483double Phase::massFraction(size_t k) const
484{
486 return m_y[k];
487}
488
489double Phase::massFraction(const string& nameSpec) const
490{
491 size_t iloc = speciesIndex(nameSpec);
492 if (iloc != npos) {
493 return massFractions()[iloc];
494 } else {
495 return 0.0;
496 }
497}
498
499void Phase::getMassFractions(double* const y) const
500{
501 copy(m_y.begin(), m_y.end(), y);
502}
503
504double Phase::concentration(const size_t k) const
505{
507 return m_y[k] * m_dens * m_rmolwts[k];
508}
509
510void Phase::getConcentrations(double* const c) const
511{
512 scale(m_ym.begin(), m_ym.end(), c, m_dens);
513}
514
515void Phase::setConcentrations(const double* const conc)
516{
517 assertCompressible("setConcentrations");
518
519 // Use m_y as temporary storage for non-negative concentrations
520 double sum = 0.0, norm = 0.0;
521 for (size_t k = 0; k != m_kk; ++k) {
522 double ck = std::max(conc[k], 0.0); // Ignore negative concentrations
523 m_y[k] = ck;
524 sum += ck * m_molwts[k];
525 norm += ck;
526 }
527 m_mmw = sum/norm;
528 setDensity(sum);
529 double rsum = 1.0/sum;
530 for (size_t k = 0; k != m_kk; ++k) {
531 m_ym[k] = m_y[k] * rsum;
532 m_y[k] = m_ym[k] * m_molwts[k]; // m_y is now the mass fraction
533 }
535}
536
537void Phase::setConcentrationsNoNorm(const double* const conc)
538{
539 assertCompressible("setConcentrationsNoNorm");
540
541 double sum = 0.0, norm = 0.0;
542 for (size_t k = 0; k != m_kk; ++k) {
543 sum += conc[k] * m_molwts[k];
544 norm += conc[k];
545 }
546 m_mmw = sum/norm;
547 setDensity(sum);
548 double rsum = 1.0/sum;
549 for (size_t k = 0; k != m_kk; ++k) {
550 m_ym[k] = conc[k] * rsum;
551 m_y[k] = m_ym[k] * m_molwts[k];
552 }
554}
555
556void Phase::setMolesNoTruncate(const double* const N)
557{
558 // get total moles
559 copy(N, N + m_kk, m_ym.begin());
560 double totalMoles = accumulate(m_ym.begin(), m_ym.end(), 0.0);
561 // get total mass
562 copy(N, N + m_kk, m_y.begin());
563 transform(m_y.begin(), m_y.end(), m_molwts.begin(), m_y.begin(), multiplies<double>());
564 double totalMass = accumulate(m_y.begin(), m_y.end(), 0.0);
565 // mean molecular weight
566 m_mmw = totalMass/totalMoles;
567 // mass fractions
568 scale(m_y.begin(), m_y.end(), m_y.begin(), 1/totalMass);
569 // moles fractions/m_mmw
570 scale(m_ym.begin(), m_ym.end(), m_ym.begin(), 1/(m_mmw * totalMoles));
571 // composition has changed
573}
574
575double Phase::elementalMassFraction(const size_t m) const
576{
578 double Z_m = 0.0;
579 for (size_t k = 0; k != m_kk; ++k) {
580 Z_m += nAtoms(k, m) * atomicWeight(m) / molecularWeight(k)
581 * massFraction(k);
582 }
583 return Z_m;
584}
585
586double Phase::elementalMoleFraction(const size_t m) const
587{
589 double denom = 0;
590 for (size_t k = 0; k < m_kk; k++) {
591 double atoms = 0;
592 for (size_t j = 0; j < nElements(); j++) {
593 atoms += nAtoms(k, j);
594 }
595 denom += atoms * moleFraction(k);
596 }
597 double numerator = 0.0;
598 for (size_t k = 0; k != m_kk; ++k) {
599 numerator += nAtoms(k, m) * moleFraction(k);
600 }
601 return numerator / denom;
602}
603
605{
606 return density()/meanMolecularWeight();
607}
608
609double Phase::molarVolume() const
610{
611 return 1.0/molarDensity();
612}
613
614void Phase::setDensity(const double density_)
615{
616 assertCompressible("setDensity");
617 if (density_ > 0.0) {
618 m_dens = density_;
619 } else {
620 throw CanteraError("Phase::setDensity",
621 "density must be positive. density = {}", density_);
622 }
623}
624
625void Phase::assignDensity(const double density_)
626{
627 if (density_ > 0.0) {
628 m_dens = density_;
629 } else {
630 throw CanteraError("Phase::assignDensity",
631 "density must be positive. density = {}", density_);
632 }
633}
634
636{
637 double cdens = 0.0;
638 for (size_t k = 0; k < m_kk; k++) {
639 cdens += charge(k)*moleFraction(k);
640 }
641 return cdens * Faraday;
642}
643
644double Phase::mean_X(const double* const Q) const
645{
646 return m_mmw*std::inner_product(m_ym.begin(), m_ym.end(), Q, 0.0);
647}
648
649double Phase::mean_X(const vector<double>& Q) const
650{
651 return m_mmw*std::inner_product(m_ym.begin(), m_ym.end(), Q.begin(), 0.0);
652}
653
654double Phase::sum_xlogx() const
655{
656 double sumxlogx = 0;
657 for (size_t k = 0; k < m_kk; k++) {
658 sumxlogx += m_ym[k] * std::log(std::max(m_ym[k], SmallNumber));
659 }
660 return m_mmw * sumxlogx + std::log(m_mmw);
661}
662
663size_t Phase::addElement(const string& symbol, double weight, int atomic_number,
664 double entropy298, int elem_type)
665{
666 // Look up the atomic weight if not given
667 if (weight == 0.0) {
668 try {
669 weight = getElementWeight(symbol);
670 } catch (CanteraError&) {
671 // assume this is just a custom element with zero atomic weight
672 }
673 } else if (weight == -12345.0) {
674 weight = getElementWeight(symbol);
675 }
676
677 // Try to look up the standard entropy if not given. Fail silently.
678 if (entropy298 == ENTROPY298_UNKNOWN) {
679 try {
680 const static AnyMap db = AnyMap::fromYamlFile(
681 "element-standard-entropies.yaml");
682 const AnyMap& elem = db["elements"].getMapWhere("symbol", symbol);
683 entropy298 = elem.convert("entropy298", "J/kmol/K", ENTROPY298_UNKNOWN);
684 } catch (CanteraError&) {
685 }
686 }
687
688 // Check for duplicates
689 auto iter = find(m_elementNames.begin(), m_elementNames.end(), symbol);
690 if (iter != m_elementNames.end()) {
691 size_t m = iter - m_elementNames.begin();
692 if (m_atomicWeights[m] != weight) {
693 throw CanteraError("Phase::addElement",
694 "Duplicate elements ({}) have different weights", symbol);
695 } else {
696 // Ignore attempt to add duplicate element with the same weight
697 return m;
698 }
699 }
700
701 // Add the new element
702 m_atomicWeights.push_back(weight);
703 m_elementNames.push_back(symbol);
704 m_atomicNumbers.push_back(atomic_number);
705 m_entropy298.push_back(entropy298);
706 if (symbol == "E") {
708 } else {
709 m_elem_type.push_back(elem_type);
710 }
711 m_mm++;
712
713 // Update species compositions
714 if (m_kk) {
715 vector<double> old(m_speciesComp);
716 m_speciesComp.resize(m_kk*m_mm, 0.0);
717 for (size_t k = 0; k < m_kk; k++) {
718 size_t m_old = m_mm - 1;
719 for (size_t m = 0; m < m_old; m++) {
720 m_speciesComp[k * m_mm + m] = old[k * (m_old) + m];
721 }
722 m_speciesComp[k * (m_mm) + (m_mm-1)] = 0.0;
723 }
725
726 return m_mm-1;
727}
728
729bool Phase::addSpecies(shared_ptr<Species> spec)
730{
731 if (m_nSpeciesLocks) {
732 throw CanteraError("Phase::addSpecies",
733 "Cannot add species to ThermoPhase '{}' because it is being "
734 "used by another object,\nsuch as a Reactor, Domain1D (flame), "
735 "SolutionArray, or MultiPhase (Mixture) object.", m_name);
736 }
737
738 if (m_species.find(spec->name) != m_species.end()) {
739 throw CanteraError("Phase::addSpecies",
740 "Phase '{}' already contains a species named '{}'.",
741 m_name, spec->name);
742 }
743
744 vector<double> comp(nElements());
745 for (const auto& [eName, stoich] : spec->composition) {
746 size_t m = elementIndex(eName);
747 if (m == npos) { // Element doesn't exist in this phase
749 case UndefElement::ignore:
750 return false;
751
752 case UndefElement::add:
753 addElement(eName);
754 comp.resize(nElements());
755 m = elementIndex(eName);
756 break;
757
758 case UndefElement::error:
759 default:
760 throw CanteraError("Phase::addSpecies",
761 "Species '{}' contains an undefined element '{}'.",
762 spec->name, eName);
763 }
764 }
765 comp[m] = stoich;
766 }
767
768 size_t ne = nElements();
769 const vector<double>& aw = atomicWeights();
770 if (spec->charge != 0.0) {
771 size_t eindex = elementIndex("E");
772 if (eindex != npos) {
773 double ecomp = comp[eindex];
774 if (fabs(spec->charge + ecomp) > 0.001) {
775 if (ecomp != 0.0) {
776 throw CanteraError("Phase::addSpecies",
777 "Input charge and element E compositions differ "
778 "for species " + spec->name);
779 } else {
780 // Just fix up the element E composition based on the input
781 // species charge
782 comp[eindex] = -spec->charge;
783 }
784 }
785 } else {
786 addElement("E", 0.000545, 0, 0.0, CT_ELEM_TYPE_ELECTRONCHARGE);
787 ne = nElements();
788 eindex = elementIndex("E");
789 comp.resize(ne);
790 comp[ne - 1] = - spec->charge;
791 }
792 }
793
794 double wt = 0.0;
795 for (size_t m = 0; m < ne; m++) {
796 wt += comp[m] * aw[m];
797 }
798
799 // Some surface phases may define species representing empty sites
800 // that have zero molecular weight. Give them a very small molecular
801 // weight to avoid dividing by zero.
802 wt = std::max(wt, Tiny);
803
804 spec->setMolecularWeight(wt);
805
806 m_molwts.push_back(wt);
807 m_rmolwts.push_back(1.0/wt);
808
809 for (size_t m = 0; m < ne; m++) {
810 m_speciesComp.push_back(comp[m]);
811 }
812
813 m_speciesNames.push_back(spec->name);
814 m_species[spec->name] = spec;
815 m_speciesIndices[spec->name] = m_kk;
816 m_speciesCharge.push_back(spec->charge);
817
818 string nLower = toLowerCopy(spec->name);
819 if (m_speciesLower.find(nLower) == m_speciesLower.end()) {
820 m_speciesLower[nLower] = m_kk;
821 } else {
822 m_speciesLower[nLower] = npos;
823 }
824 m_kk++;
825
826 // Ensure that the Phase has a valid mass fraction vector that sums to
827 // one. We will assume that species 0 has a mass fraction of 1.0 and mass
828 // fraction of all other species is 0.0.
829 if (m_kk == 1) {
830 m_y.push_back(1.0);
831 m_ym.push_back(m_rmolwts[0]);
832 m_mmw = 1.0 / m_ym[0];
833 } else {
834 m_y.push_back(0.0);
835 m_ym.push_back(0.0);
836 }
837 m_workS.push_back(0.0);
839 return true;
840}
841
842void Phase::modifySpecies(size_t k, shared_ptr<Species> spec)
843{
844 if (speciesName(k) != spec->name) {
845 throw CanteraError("Phase::modifySpecies",
846 "New species name '{}' does not match existing name '{}'",
847 spec->name, speciesName(k));
848 }
849 const shared_ptr<Species>& old = m_species[spec->name];
850 if (spec->composition != old->composition) {
851 throw CanteraError("Phase::modifySpecies",
852 "New composition for '{}' does not match existing composition",
853 spec->name);
854 }
855 m_species[spec->name] = spec;
857}
858
859void Phase::addSpeciesAlias(const string& name, const string& alias)
860{
861 if (speciesIndex(alias) != npos) {
862 throw CanteraError("Phase::addSpeciesAlias",
863 "Invalid alias '{}': species already exists", alias);
864 }
865 size_t k = speciesIndex(name);
866 if (k != npos) {
867 m_speciesIndices[alias] = k;
868 } else {
869 throw CanteraError("Phase::addSpeciesAlias",
870 "Unable to add alias '{}' "
871 "(original species '{}' not found).", alias, name);
872 }
873}
874
876{
877 if (!m_nSpeciesLocks) {
878 throw CanteraError("Phase::removeSpeciesLock",
879 "ThermoPhase '{}' has no current species locks.", m_name);
880 }
882}
883
884vector<string> Phase::findIsomers(const Composition& compMap) const
885{
886 vector<string> isomerNames;
887
888 for (const auto& [name, species] : m_species) {
889 if (species->composition == compMap) {
890 isomerNames.emplace_back(name);
891 }
892 }
893
894 return isomerNames;
895}
896
897vector<string> Phase::findIsomers(const string& comp) const
898{
899 return findIsomers(parseCompString(comp));
900}
901
902shared_ptr<Species> Phase::species(const string& name) const
903{
904 size_t k = speciesIndex(name);
905 if (k != npos) {
906 return m_species.at(speciesName(k));
907 } else {
908 throw CanteraError("Phase::species",
909 "Unknown species '{}'", name);
910 }
911}
912
913shared_ptr<Species> Phase::species(size_t k) const
914{
916 return m_species.at(m_speciesNames[k]);
917}
918
920 m_undefinedElementBehavior = UndefElement::ignore;
921}
922
924 m_undefinedElementBehavior = UndefElement::add;
925}
926
928 m_undefinedElementBehavior = UndefElement::error;
929}
930
931bool Phase::ready() const
932{
933 return (m_kk > 0);
934}
935
937 m_stateNum++;
938 m_cache.clear();
939}
940
941void Phase::setMolecularWeight(const int k, const double mw)
942{
943 m_molwts[k] = mw;
944 m_rmolwts[k] = 1.0/mw;
945
946 transform(m_y.begin(), m_y.end(), m_rmolwts.begin(), m_ym.begin(),
947 multiplies<double>());
948 m_mmw = 1.0 / accumulate(m_ym.begin(), m_ym.end(), 0.0);
949}
950
952 m_stateNum++;
953}
954
955vector<double> Phase::getCompositionFromMap(const Composition& comp) const
956{
957 vector<double> X(m_kk);
958 for (const auto& [name, value] : comp) {
959 size_t loc = speciesIndex(name);
960 if (loc == npos) {
961 throw CanteraError("Phase::getCompositionFromMap",
962 "Unknown species '{}'", name);
963 }
964 X[loc] = value;
965 }
966 return X;
967}
968
969void Phase::massFractionsToMoleFractions(const double* Y, double* X) const
970{
971 double rmmw = 0.0;
972 for (size_t k = 0; k != m_kk; ++k) {
973 rmmw += Y[k]/m_molwts[k];
974 }
975 if (rmmw == 0.0) {
976 throw CanteraError("Phase::massFractionsToMoleFractions",
977 "no input composition given");
978 }
979 for (size_t k = 0; k != m_kk; ++k) {
980 X[k] = Y[k]/(rmmw*m_molwts[k]);
981 }
982}
983
984void Phase::moleFractionsToMassFractions(const double* X, double* Y) const
985{
986 double mmw = dot(X, X+m_kk, m_molwts.data());
987 if (mmw == 0.0) {
988 throw CanteraError("Phase::moleFractionsToMassFractions",
989 "no input composition given");
990 }
991 double rmmw = 1.0/mmw;
992 for (size_t k = 0; k != m_kk; ++k) {
993 Y[k] = X[k]*m_molwts[k]*rmmw;
994 }
995}
996
997} // namespace Cantera
#define CT_ELEM_TYPE_ELECTRONCHARGE
This refers to conservation of electrons.
Definition Elements.h:41
#define ENTROPY298_UNKNOWN
Number indicating we don't know the entropy of the element in its most stable state at 298....
Definition Elements.h:85
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.
Definition AnyMap.h:431
double convert(const string &key, const string &units) const
Convert the item stored by the given key to the units specified in units.
Definition AnyMap.cpp:1595
static AnyMap fromYamlFile(const string &name, const string &parent_name="")
Create an AnyMap from a YAML file.
Definition AnyMap.cpp:1841
Array size error.
Base class for exceptions thrown by Cantera classes.
An array index is out of range.
virtual vector< string > partialStates() const
Return a vector of settable partial property sets within a phase.
Definition Phase.cpp:215
void getCharges(double *charges) const
Copy the vector of species charges into array charges.
Definition Phase.cpp:433
virtual void getConcentrations(double *const c) const
Get the species concentrations (kmol/m^3).
Definition Phase.cpp:510
map< string, size_t > m_speciesLower
Map of lower-case species names to indices.
Definition Phase.h:952
double massFraction(size_t k) const
Return the mass fraction of a single species.
Definition Phase.cpp:483
virtual double molarDensity() const
Molar density (kmol/m^3).
Definition Phase.cpp:604
void assignDensity(const double density_)
Set the internally stored constant density (kg/m^3) of the phase.
Definition Phase.cpp:625
virtual bool addSpecies(shared_ptr< Species > spec)
Add a Species to this Phase.
Definition Phase.cpp:729
virtual void setMoleFractions(const double *const x)
Set the mole fractions to the specified values.
Definition Phase.cpp:310
int changeElementType(int m, int elem_type)
Change the element type of the mth constraint Reassigns an element type.
Definition Phase.cpp:96
const vector< double > & atomicWeights() const
Return a read-only reference to the vector of atomic weights.
Definition Phase.cpp:81
virtual vector< string > fullStates() const
Return a vector containing full states defining a phase.
Definition Phase.cpp:197
void assertCompressible(const string &setter) const
Ensure that phase is compressible.
Definition Phase.h:839
vector< double > m_workS
Vector of size m_kk, used as a temporary holding area.
Definition Phase.h:903
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range.
Definition Phase.cpp:153
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
Definition Phase.cpp:286
vector< double > m_speciesComp
Atomic composition of the species.
Definition Phase.h:888
ValueCache m_cache
Cached for saved calculations within each ThermoPhase.
Definition Phase.h:859
virtual void restorePartialState(size_t lenstate, const double *state)
Set the internal thermodynamic state of the phase, excluding composition.
Definition Phase.cpp:243
size_t m_nSpeciesLocks
Reference counter preventing species addition.
Definition Phase.h:894
vector< string > m_speciesNames
Vector of the species names.
Definition Phase.h:946
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:232
bool m_caseSensitiveSpecies
Flag determining whether case sensitive species names are enforced.
Definition Phase.h:900
vector< string > m_elementNames
element names
Definition Phase.h:957
void ignoreUndefinedElements()
Set behavior when adding a species containing undefined elements to just skip the species.
Definition Phase.cpp:919
UndefElement::behavior m_undefinedElementBehavior
Flag determining behavior when adding species with an undefined element.
Definition Phase.h:897
virtual void setMassFractions_NoNorm(const double *const y)
Set the mass fractions to the specified values without normalizing.
Definition Phase.cpp:376
virtual map< string, size_t > nativeState() const
Return a map of properties defining the native state of a substance.
Definition Phase.cpp:167
double chargeDensity() const
Charge density [C/m^3].
Definition Phase.cpp:635
void addUndefinedElements()
Set behavior when adding a species containing undefined elements to add those elements to the phase.
Definition Phase.cpp:923
virtual void setConcentrationsNoNorm(const double *const conc)
Set the concentrations without ignoring negative concentrations.
Definition Phase.cpp:537
vector< int > m_atomicNumbers
element atomic numbers
Definition Phase.h:956
size_t m_kk
Number of species in the phase.
Definition Phase.h:879
int atomicNumber(size_t m) const
Atomic number of element m.
Definition Phase.cpp:86
virtual void modifySpecies(size_t k, shared_ptr< Species > spec)
Modify the thermodynamic data associated with a species.
Definition Phase.cpp:842
double m_mmw
mean molecular weight of the mixture (kg kmol-1)
Definition Phase.h:924
double elementalMoleFraction(const size_t m) const
Elemental mole fraction of element m.
Definition Phase.cpp:586
void setState_TD(double t, double rho)
Set the internally stored temperature (K) and density (kg/m^3)
Definition Phase.cpp:398
vector< double > m_rmolwts
inverse of species molecular weights (kmol kg-1)
Definition Phase.h:939
double temperature() const
Temperature (K).
Definition Phase.h:587
virtual void setPressure(double p)
Set the internally stored pressure (Pa) at constant temperature and composition.
Definition Phase.h:641
virtual bool isCompressible() const
Return whether phase represents a compressible substance.
Definition Phase.h:258
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:680
void moleFractionsToMassFractions(const double *X, double *Y) const
Converts a mixture composition from mass fractions to mole fractions.
Definition Phase.cpp:984
virtual void setConcentrations(const double *const conc)
Set the concentrations to the specified values within the phase.
Definition Phase.cpp:515
void removeSpeciesLock()
Decrement species lock counter.
Definition Phase.cpp:875
void saveState(vector< double > &state) const
Save the current internal state of the phase.
Definition Phase.cpp:266
Composition getMoleFractionsByName(double threshold=0.0) const
Get the mole fractions by name.
Definition Phase.cpp:438
size_t elementIndex(const string &name) const
Return the index of element named 'name'.
Definition Phase.cpp:55
virtual double concentration(const size_t k) const
Concentration of species k.
Definition Phase.cpp:504
double atomicWeight(size_t m) const
Atomic weight of element m.
Definition Phase.cpp:70
void checkElementArraySize(size_t mm) const
Check that an array size is at least nElements().
Definition Phase.cpp:42
void setMassFractionsByName(const Composition &yMap)
Set the species mass fractions by name.
Definition Phase.cpp:387
int elementType(size_t m) const
Return the element constraint type Possible types include:
Definition Phase.cpp:91
string speciesName(size_t k) const
Name of the species with index k.
Definition Phase.cpp:142
map< string, size_t > m_speciesIndices
Map of species names to indices.
Definition Phase.h:949
virtual void setDensity(const double density_)
Set the internally stored density (kg/m^3) of the phase.
Definition Phase.cpp:614
Composition getMassFractionsByName(double threshold=0.0) const
Get the mass fractions by name.
Definition Phase.cpp:450
virtual size_t stateSize() const
Return size of vector defining internal state of the phase.
Definition Phase.cpp:258
string nativeMode() const
Return string acronym representing the native state of a Phase.
Definition Phase.cpp:184
vector< double > getCompositionFromMap(const Composition &comp) const
Converts a Composition to a vector with entries for each species Species that are not specified are s...
Definition Phase.cpp:955
map< string, shared_ptr< Species > > m_species
Map of Species objects.
Definition Phase.h:892
size_t findSpeciesLower(const string &nameStr) const
Find lowercase species name in m_speciesIndices when case sensitive species names are not enforced an...
Definition Phase.cpp:110
vector< double > m_molwts
species molecular weights (kg kmol-1)
Definition Phase.h:937
virtual vector< string > findIsomers(const Composition &compMap) const
Return a vector with isomers names matching a given composition map.
Definition Phase.cpp:884
const vector< double > & inverseMolecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition Phase.cpp:428
virtual bool isPure() const
Return whether phase represents a pure (single species) substance.
Definition Phase.h:248
vector< double > m_y
Mass fractions of the species.
Definition Phase.h:935
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition Phase.cpp:462
void setMoleFractionsByName(const Composition &xMap)
Set the species mole fractions by name.
Definition Phase.cpp:351
const double * massFractions() const
Return a const pointer to the mass fraction array.
Definition Phase.h:467
vector< int > m_elem_type
Vector of element types.
Definition Phase.h:958
double sum_xlogx() const
Evaluate .
Definition Phase.cpp:654
string m_name
Name of the phase.
Definition Phase.h:914
const vector< double > & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition Phase.cpp:423
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition Phase.cpp:129
double moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition Phase.cpp:467
double m_dens
Density (kg m-3).
Definition Phase.h:922
const vector< string > & elementNames() const
Return a read-only reference to the vector of element names.
Definition Phase.cpp:65
virtual double density() const
Density (kg/m^3).
Definition Phase.h:612
virtual void compositionChanged()
Apply changes to the state which are needed after the composition changes.
Definition Phase.cpp:951
vector< double > m_atomicWeights
element atomic weights (kg kmol-1)
Definition Phase.h:955
void checkSpeciesArraySize(size_t kk) const
Check that an array size is at least nSpecies().
Definition Phase.cpp:160
void getMolecularWeights(double *weights) const
Copy the vector of molecular weights into array weights.
Definition Phase.cpp:417
double nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
Definition Phase.cpp:103
void addSpeciesAlias(const string &name, const string &alias)
Add a species alias (that is, a user-defined alternative species name).
Definition Phase.cpp:859
virtual void setMolesNoTruncate(const double *const N)
Set the state of the object with moles in [kmol].
Definition Phase.cpp:556
virtual void setTemperature(double temp)
Set the internally stored temperature of the phase (K).
Definition Phase.h:648
size_t nElements() const
Number of elements.
Definition Phase.cpp:30
void setMolecularWeight(const int k, const double mw)
Set the molecular weight of a single species to a given value.
Definition Phase.cpp:941
double mean_X(const double *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition Phase.cpp:644
vector< double > m_entropy298
Entropy at 298.15 K and 1 bar of stable state pure elements (J kmol-1)
Definition Phase.h:961
virtual void savePartialState(size_t lenstate, double *state) const
Save the current thermodynamic state of the phase, excluding composition.
Definition Phase.cpp:228
vector< double > m_ym
m_ym[k] = mole fraction of species k divided by the mean molecular weight of mixture.
Definition Phase.h:928
virtual void setMassFractions(const double *const y)
Set the mass fractions to the specified values and normalize them.
Definition Phase.cpp:362
const vector< string > & speciesNames() const
Return a const reference to the vector of species names.
Definition Phase.cpp:148
virtual bool ready() const
Returns a bool indicating whether the object is ready for use.
Definition Phase.cpp:931
double molecularWeight(size_t k) const
Molecular weight of species k.
Definition Phase.cpp:411
double elementalMassFraction(const size_t m) const
Elemental mass fraction of element m.
Definition Phase.cpp:575
shared_ptr< Species > species(const string &name) const
Return the Species object for the named species.
Definition Phase.cpp:902
virtual double molarVolume() const
Molar volume (m^3/kmol).
Definition Phase.cpp:609
virtual void invalidateCache()
Invalidate any cached values which are normally updated only when a change in state is detected.
Definition Phase.cpp:936
void checkElementIndex(size_t m) const
Check that the specified element index is in range.
Definition Phase.cpp:35
void getMassFractions(double *const y) const
Get the species mass fractions.
Definition Phase.cpp:499
virtual size_t partialStateSize() const
Get the size of the partial state vector of the phase.
Definition Phase.h:305
int m_stateNum
State Change variable.
Definition Phase.h:943
void throwUndefinedElements()
Set the behavior when adding a species containing undefined elements to throw an exception.
Definition Phase.cpp:927
void setName(const string &nm)
Sets the string name for the phase.
Definition Phase.cpp:25
size_t m_mm
Number of elements.
Definition Phase.h:954
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition Phase.h:605
string elementName(size_t m) const
Name of the element with index m.
Definition Phase.cpp:49
double charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition Phase.h:563
void massFractionsToMoleFractions(const double *Y, double *X) const
Converts a mixture composition from mole fractions to mass fractions.
Definition Phase.cpp:969
double entropyElement298(size_t m) const
Entropy of the element in its standard state at 298 K and 1 bar.
Definition Phase.cpp:75
virtual void setMoleFractions_NoNorm(const double *const x)
Set the mole fractions to the specified values without normalizing.
Definition Phase.cpp:342
vector< double > m_speciesCharge
Vector of species charges.
Definition Phase.h:890
size_t addElement(const string &symbol, double weight=-12345.0, int atomicNumber=0, double entropy298=ENTROPY298_UNKNOWN, int elem_type=CT_ELEM_TYPE_ABSPOS)
Add an element.
Definition Phase.cpp:663
string name() const
Return the name of the phase.
Definition Phase.cpp:20
void clear()
Clear all cached values.
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.
Definition utilities.h:82
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition utilities.h:104
const double Faraday
Faraday constant [C/kmol].
Definition ct_defs.h:131
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:180
double getElementWeight(const string &ename)
Get the atomic weight of an element.
Definition Elements.cpp:251
const double Tiny
Small number to compare differences of mole fractions against.
Definition ct_defs.h:173
const double SmallNumber
smallest number to compare to zero.
Definition ct_defs.h:158
map< string, double > Composition
Map from string names to doubles.
Definition ct_defs.h:177
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...