Cantera  3.0.0
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-1);
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
110void Phase::getAtoms(size_t k, double* atomArray) const
111{
112 warn_deprecated("Phase::getAtoms", "To be removed after Cantera 3.0. "
113 "Use 'nAtoms(species_index, element_index)' instead.");
114 for (size_t m = 0; m < m_mm; m++) {
115 atomArray[m] = (double) m_speciesComp[m_mm * k + m];
116 }
117}
118
119size_t Phase::findSpeciesLower(const string& name) const
120{
121 string nLower = toLowerCopy(name);
122 size_t loc = npos;
123 auto it = m_speciesLower.find(nLower);
124 if (it == m_speciesLower.end()) {
125 return npos;
126 } else {
127 loc = it->second;
128 if (loc == npos) {
129 throw CanteraError("Phase::findSpeciesLower",
130 "Lowercase species name '{}' is not unique. "
131 "Set Phase::caseSensitiveSpecies to true to "
132 "enforce case sensitive species names", nLower);
133 }
134 }
135 return loc;
136}
137
138size_t Phase::speciesIndex(const string& nameStr) const
139{
140 size_t loc = npos;
141
142 auto it = m_speciesIndices.find(nameStr);
143 if (it != m_speciesIndices.end()) {
144 return it->second;
145 } else if (!m_caseSensitiveSpecies) {
146 loc = findSpeciesLower(nameStr);
147 }
148 return loc;
149}
150
151string Phase::speciesName(size_t k) const
152{
154 return m_speciesNames[k];
155}
156
157const vector<string>& Phase::speciesNames() const
158{
159 return m_speciesNames;
160}
161
162void Phase::checkSpeciesIndex(size_t k) const
163{
164 if (k >= m_kk) {
165 throw IndexError("Phase::checkSpeciesIndex", "species", k, m_kk-1);
166 }
167}
168
169void Phase::checkSpeciesArraySize(size_t kk) const
170{
171 if (m_kk > kk) {
172 throw ArraySizeError("Phase::checkSpeciesArraySize", kk, m_kk);
173 }
174}
175
176string Phase::speciesSPName(int k) const
177{
178 warn_deprecated("Phase::speciesSPName", "To be removed after Cantera 3.0");
179 return m_name + ":" + speciesName(k);
180}
181
182map<string, size_t> Phase::nativeState() const
183{
184 if (isPure()) {
185 if (isCompressible()) {
186 return { {"T", 0}, {"D", 1} };
187 } else {
188 return { {"T", 0}, {"P", 1} };
189 }
190 } else {
191 if (isCompressible()) {
192 return { {"T", 0}, {"D", 1}, {"Y", 2} };
193 } else {
194 return { {"T", 0}, {"P", 1}, {"Y", 2} };
195 }
196 }
197}
198
199string Phase::nativeMode() const
200{
201 map<size_t, string> states; // reverse lookup
202 for (auto& [name, value] : nativeState()) {
203 states[value] = name;
204 }
205 string out;
206 for (auto& [value, name] : states) {
207 out += name;
208 }
209 return out;
210}
211
212vector<string> Phase::fullStates() const
213{
214 if (isPure()) {
215 if (isCompressible()) {
216 return {"TD", "TP", "UV", "DP", "HP", "SP", "SV"};
217 } else {
218 return {"TP", "HP", "SP"};
219 }
220 } else {
221 if (isCompressible()) {
222 return {"TDX", "TDY", "TPX", "TPY", "UVX", "UVY", "DPX", "DPY",
223 "HPX", "HPY", "SPX", "SPY", "SVX", "SVY"};
224 } else {
225 return {"TPX", "TPY", "HPX", "HPY", "SPX", "SPY"};
226 }
227 }
228}
229
230vector<string> Phase::partialStates() const
231{
232 if (isPure()) {
233 return {};
234 } else {
235 if (isCompressible()) {
236 return {"TD", "TP", "UV", "DP", "HP", "SP", "SV"};
237 } else {
238 return {"TP", "HP", "SP"};
239 }
240 }
241}
242
243size_t Phase::stateSize() const {
244 if (isPure()) {
245 return 2;
246 } else {
247 return nSpecies() + 2;
248 }
249}
250
251void Phase::saveState(vector<double>& state) const
252{
253 state.resize(stateSize());
254 saveState(state.size(), &state[0]);
255}
256
257void Phase::saveState(size_t lenstate, double* state) const
258{
259 auto native = nativeState();
260
261 // function assumes default definition of nativeState
262 state[native.at("T")] = temperature();
263 if (isCompressible()) {
264 state[native.at("D")] = density();
265 } else {
266 state[native.at("P")] = pressure();
267 }
268 if (native.count("X")) {
269 getMoleFractions(state + native["X"]);
270 } else if (native.count("Y")) {
271 getMassFractions(state + native["Y"]);
272 }
273}
274
275void Phase::restoreState(const vector<double>& state)
276{
277 restoreState(state.size(),&state[0]);
278}
279
280void Phase::restoreState(size_t lenstate, const double* state)
281{
282 size_t ls = stateSize();
283 if (lenstate < ls) {
284 throw ArraySizeError("Phase::restoreState",
285 lenstate, ls);
286 }
287
288 auto native = nativeState();
289 setTemperature(state[native.at("T")]);
290 if (isCompressible()) {
291 setDensity(state[native.at("D")]);
292 } else {
293 setPressure(state[native.at("P")]);
294 }
295
296 if (native.count("X")) {
297 setMoleFractions_NoNorm(state + native["X"]);
298 } else if (native.count("Y")) {
299 setMassFractions_NoNorm(state + native["Y"]);
300 }
302}
303
304void Phase::setMoleFractions(const double* const x)
305{
306 // Use m_y as a temporary work vector for the non-negative mole fractions
307 double norm = 0.0;
308 // sum is calculated below as the unnormalized molecular weight
309 double sum = 0;
310 for (size_t k = 0; k < m_kk; k++) {
311 double xk = std::max(x[k], 0.0); // Ignore negative mole fractions
312 m_y[k] = xk;
313 norm += xk;
314 sum += m_molwts[k] * xk;
315 }
316
317 // Set m_ym to the normalized mole fractions divided by the normalized mean
318 // molecular weight:
319 // m_ym_k = X_k / (sum_k X_k M_k)
320 const double invSum = 1.0/sum;
321 for (size_t k=0; k < m_kk; k++) {
322 m_ym[k] = m_y[k]*invSum;
323 }
324
325 // Now set m_y to the normalized mass fractions:
326 // m_y = X_k M_k / (sum_k X_k M_k)
327 for (size_t k=0; k < m_kk; k++) {
328 m_y[k] = m_ym[k] * m_molwts[k];
329 }
330
331 // Calculate the normalized molecular weight
332 m_mmw = sum/norm;
334}
335
336void Phase::setMoleFractions_NoNorm(const double* const x)
337{
338 m_mmw = dot(x, x + m_kk, m_molwts.begin());
339 scale(x, x + m_kk, m_ym.begin(), 1.0/m_mmw);
340 transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(),
341 m_y.begin(), multiplies<double>());
343}
344
346{
347 vector<double> mf = getCompositionFromMap(xMap);
348 setMoleFractions(mf.data());
349}
350
351void Phase::setMoleFractionsByName(const string& x)
352{
354}
355
356void Phase::setMassFractions(const double* const y)
357{
358 for (size_t k = 0; k < m_kk; k++) {
359 m_y[k] = std::max(y[k], 0.0); // Ignore negative mass fractions
360 }
361 double norm = accumulate(m_y.begin(), m_y.end(), 0.0);
362 scale(m_y.begin(), m_y.end(), m_y.begin(), 1.0/norm);
363
364 transform(m_y.begin(), m_y.end(), m_rmolwts.begin(),
365 m_ym.begin(), multiplies<double>());
366 m_mmw = 1.0 / accumulate(m_ym.begin(), m_ym.end(), 0.0);
368}
369
370void Phase::setMassFractions_NoNorm(const double* const y)
371{
372 double sum = 0.0;
373 copy(y, y + m_kk, m_y.begin());
374 transform(m_y.begin(), m_y.end(), m_rmolwts.begin(), m_ym.begin(),
375 multiplies<double>());
376 sum = accumulate(m_ym.begin(), m_ym.end(), 0.0);
377 m_mmw = 1.0/sum;
379}
380
382{
383 vector<double> mf = getCompositionFromMap(yMap);
384 setMassFractions(mf.data());
385}
386
387void Phase::setMassFractionsByName(const string& y)
388{
390}
391
392void Phase::setState_TRX(double t, double dens, const double* x)
393{
394 warn_deprecated("Phase::setState_TRX",
395 "To be removed after Cantera 3.0. Replaceable by calls to "
396 "setMoleFractions and setState_TD.");
398 setState_TD(t, dens);
399}
400
401void Phase::setState_TNX(double t, double n, const double* x)
402{
403 warn_deprecated("Phase::setState_TNX", "To be removed after Cantera 3.0. "
404 "Use 'setMoleFractions' and 'setState_TD' instead.");
408}
409
410void Phase::setState_TRX(double t, double dens, const Composition& x)
411{
412 warn_deprecated("Phase::setState_TRX",
413 "To be removed after Cantera 3.0. Replaceable by calls to "
414 "setMoleFractionsByName and setState_TD.");
416 setState_TD(t, dens);
417}
418
419void Phase::setState_TRY(double t, double dens, const double* y)
420{
421 warn_deprecated("Phase::setState_TRY",
422 "To be removed after Cantera 3.0. Replaceable by calls to "
423 "setMassFractions and setState_TD.");
425 setState_TD(t, dens);
426}
427
428void Phase::setState_TRY(double t, double dens, const Composition& y)
429{
430 warn_deprecated("Phase::setState_TRY",
431 "To be removed after Cantera 3.0. Replaceable by calls to "
432 "setMassFractionsByName and setState_TD.");
434 setState_TD(t, dens);
435}
436
437void Phase::setState_TR(double t, double rho)
438{
439 warn_deprecated("Phase::setState_TR",
440 "To be removed after Cantera 3.0. Renamed to setState_TD.");
441 setState_TD(t, rho);
442}
443
444void Phase::setState_TD(double t, double rho)
445{
447 setDensity(rho);
448}
449
450void Phase::setState_TX(double t, double* x)
451{
452 warn_deprecated("Phase::setState_TX", "To be removed after Cantera 3.0. "
453 "Use calls to 'setTemperature' and 'setMoleFractions' instead.");
456}
457
458void Phase::setState_TY(double t, double* y)
459{
460 warn_deprecated("Phase::setState_TY", "To be removed after Cantera 3.0. "
461 "Use calls to 'setTemperature' and 'setMassFractions' instead.");
464}
465
466void Phase::setState_RX(double rho, double* x)
467{
468 warn_deprecated("Phase::setState_RX", "To be removed after Cantera 3.0. "
469 "Use calls to 'setDensity' and 'setMoleFractions' instead.");
471 setDensity(rho);
472}
473
474void Phase::setState_RY(double rho, double* y)
475{
476 warn_deprecated("Phase::setState_RY", "To be removed after Cantera 3.0. "
477 "Use calls to 'setDensity' and 'setMassFractions' instead.");
479 setDensity(rho);
480}
481
482double Phase::molecularWeight(size_t k) const
483{
485 return m_molwts[k];
486}
487
488void Phase::getMolecularWeights(vector<double>& weights) const
489{
490 warn_deprecated("Phase::getMolecularWeights(vector<double>&)", "To be removed after "
491 "Cantera 3.0. Use 'getMolecularWeights(vec.data())' instead.");
492 weights = molecularWeights();
493}
494
495void Phase::getMolecularWeights(double* weights) const
496{
497 const vector<double>& mw = molecularWeights();
498 copy(mw.begin(), mw.end(), weights);
499}
500
501const vector<double>& Phase::molecularWeights() const
502{
503 return m_molwts;
504}
505
506const vector<double>& Phase::inverseMolecularWeights() const
507{
508 return m_rmolwts;
509}
510
511void Phase::getCharges(double* charges) const
512{
513 copy(m_speciesCharge.begin(), m_speciesCharge.end(), charges);
514}
515
517{
518 Composition comp;
519 for (size_t k = 0; k < m_kk; k++) {
520 double x = moleFraction(k);
521 if (x > threshold) {
522 comp[speciesName(k)] = x;
523 }
524 }
525 return comp;
526}
527
529{
530 Composition comp;
531 for (size_t k = 0; k < m_kk; k++) {
532 double x = massFraction(k);
533 if (x > threshold) {
534 comp[speciesName(k)] = x;
535 }
536 }
537 return comp;
538}
539
540void Phase::getMoleFractions(double* const x) const
541{
542 scale(m_ym.begin(), m_ym.end(), x, m_mmw);
543}
544
545double Phase::moleFraction(size_t k) const
546{
548 return m_ym[k] * m_mmw;
549}
550
551double Phase::moleFraction(const string& nameSpec) const
552{
553 size_t iloc = speciesIndex(nameSpec);
554 if (iloc != npos) {
555 return moleFraction(iloc);
556 } else {
557 return 0.0;
558 }
559}
560
561const double* Phase::moleFractdivMMW() const
562{
563 warn_deprecated("Phase::moleFractdivMMW", "To be removed after Cantera 3.0. "
564 "Generally replaceable by 'getMoleFractions' and 'meanMolecularWeight'.");
565 return &m_ym[0];
566}
567
568double Phase::massFraction(size_t k) const
569{
571 return m_y[k];
572}
573
574double Phase::massFraction(const string& nameSpec) const
575{
576 size_t iloc = speciesIndex(nameSpec);
577 if (iloc != npos) {
578 return massFractions()[iloc];
579 } else {
580 return 0.0;
581 }
582}
583
584void Phase::getMassFractions(double* const y) const
585{
586 copy(m_y.begin(), m_y.end(), y);
587}
588
589double Phase::concentration(const size_t k) const
590{
592 return m_y[k] * m_dens * m_rmolwts[k];
593}
594
595void Phase::getConcentrations(double* const c) const
596{
597 scale(m_ym.begin(), m_ym.end(), c, m_dens);
598}
599
600void Phase::setConcentrations(const double* const conc)
601{
602 assertCompressible("setConcentrations");
603
604 // Use m_y as temporary storage for non-negative concentrations
605 double sum = 0.0, norm = 0.0;
606 for (size_t k = 0; k != m_kk; ++k) {
607 double ck = std::max(conc[k], 0.0); // Ignore negative concentrations
608 m_y[k] = ck;
609 sum += ck * m_molwts[k];
610 norm += ck;
611 }
612 m_mmw = sum/norm;
613 setDensity(sum);
614 double rsum = 1.0/sum;
615 for (size_t k = 0; k != m_kk; ++k) {
616 m_ym[k] = m_y[k] * rsum;
617 m_y[k] = m_ym[k] * m_molwts[k]; // m_y is now the mass fraction
618 }
620}
621
622void Phase::setConcentrationsNoNorm(const double* const conc)
623{
624 assertCompressible("setConcentrationsNoNorm");
625
626 double sum = 0.0, norm = 0.0;
627 for (size_t k = 0; k != m_kk; ++k) {
628 sum += conc[k] * m_molwts[k];
629 norm += conc[k];
630 }
631 m_mmw = sum/norm;
632 setDensity(sum);
633 double rsum = 1.0/sum;
634 for (size_t k = 0; k != m_kk; ++k) {
635 m_ym[k] = conc[k] * rsum;
636 m_y[k] = m_ym[k] * m_molwts[k];
637 }
639}
640
641void Phase::setMolesNoTruncate(const double* const N)
642{
643 // get total moles
644 copy(N, N + m_kk, m_ym.begin());
645 double totalMoles = accumulate(m_ym.begin(), m_ym.end(), 0.0);
646 // get total mass
647 copy(N, N + m_kk, m_y.begin());
648 transform(m_y.begin(), m_y.end(), m_molwts.begin(), m_y.begin(), multiplies<double>());
649 double totalMass = accumulate(m_y.begin(), m_y.end(), 0.0);
650 // mean molecular weight
651 m_mmw = totalMass/totalMoles;
652 // mass fractions
653 scale(m_y.begin(), m_y.end(), m_y.begin(), 1/totalMass);
654 // moles fractions/m_mmw
655 scale(m_ym.begin(), m_ym.end(), m_ym.begin(), 1/(m_mmw * totalMoles));
656 // composition has changed
658}
659
660double Phase::elementalMassFraction(const size_t m) const
661{
663 double Z_m = 0.0;
664 for (size_t k = 0; k != m_kk; ++k) {
665 Z_m += nAtoms(k, m) * atomicWeight(m) / molecularWeight(k)
666 * massFraction(k);
667 }
668 return Z_m;
669}
670
671double Phase::elementalMoleFraction(const size_t m) const
672{
674 double denom = 0;
675 for (size_t k = 0; k < m_kk; k++) {
676 double atoms = 0;
677 for (size_t j = 0; j < nElements(); j++) {
678 atoms += nAtoms(k, j);
679 }
680 denom += atoms * moleFraction(k);
681 }
682 double numerator = 0.0;
683 for (size_t k = 0; k != m_kk; ++k) {
684 numerator += nAtoms(k, m) * moleFraction(k);
685 }
686 return numerator / denom;
687}
688
690{
691 return density()/meanMolecularWeight();
692}
693
694void Phase::setMolarDensity(const double molar_density)
695{
696 warn_deprecated("Phase::setMolarDensity", "To be removed after Cantera 3.0. "
697 "Use 'setDensity(molar_density * meanMolecularWeight())' instead.");
698 assertCompressible("setMolarDensity");
699 m_dens = molar_density*meanMolecularWeight();
700}
701
702double Phase::molarVolume() const
703{
704 return 1.0/molarDensity();
705}
706
707void Phase::setDensity(const double density_)
708{
709 assertCompressible("setDensity");
710 if (density_ > 0.0) {
711 m_dens = density_;
712 } else {
713 throw CanteraError("Phase::setDensity",
714 "density must be positive. density = {}", density_);
715 }
716}
717
718void Phase::assignDensity(const double density_)
719{
720 if (density_ > 0.0) {
721 m_dens = density_;
722 } else {
723 throw CanteraError("Phase::assignDensity",
724 "density must be positive. density = {}", density_);
725 }
726}
727
729{
730 double cdens = 0.0;
731 for (size_t k = 0; k < m_kk; k++) {
732 cdens += charge(k)*moleFraction(k);
733 }
734 return cdens * Faraday;
735}
736
737double Phase::mean_X(const double* const Q) const
738{
739 return m_mmw*std::inner_product(m_ym.begin(), m_ym.end(), Q, 0.0);
740}
741
742double Phase::mean_X(const vector<double>& Q) const
743{
744 return m_mmw*std::inner_product(m_ym.begin(), m_ym.end(), Q.begin(), 0.0);
745}
746
747double Phase::sum_xlogx() const
748{
749 double sumxlogx = 0;
750 for (size_t k = 0; k < m_kk; k++) {
751 sumxlogx += m_ym[k] * std::log(std::max(m_ym[k], SmallNumber));
752 }
753 return m_mmw * sumxlogx + std::log(m_mmw);
754}
755
756size_t Phase::addElement(const string& symbol, double weight, int atomic_number,
757 double entropy298, int elem_type)
758{
759 // Look up the atomic weight if not given
760 if (weight == 0.0) {
761 try {
762 weight = getElementWeight(symbol);
763 } catch (CanteraError&) {
764 // assume this is just a custom element with zero atomic weight
765 }
766 } else if (weight == -12345.0) {
767 weight = getElementWeight(symbol);
768 }
769
770 // Try to look up the standard entropy if not given. Fail silently.
771 if (entropy298 == ENTROPY298_UNKNOWN) {
772 try {
773 const static AnyMap db = AnyMap::fromYamlFile(
774 "element-standard-entropies.yaml");
775 const AnyMap& elem = db["elements"].getMapWhere("symbol", symbol);
776 entropy298 = elem.convert("entropy298", "J/kmol/K", ENTROPY298_UNKNOWN);
777 } catch (CanteraError&) {
778 }
779 }
780
781 // Check for duplicates
782 auto iter = find(m_elementNames.begin(), m_elementNames.end(), symbol);
783 if (iter != m_elementNames.end()) {
784 size_t m = iter - m_elementNames.begin();
785 if (m_atomicWeights[m] != weight) {
786 throw CanteraError("Phase::addElement",
787 "Duplicate elements ({}) have different weights", symbol);
788 } else {
789 // Ignore attempt to add duplicate element with the same weight
790 return m;
791 }
792 }
793
794 // Add the new element
795 m_atomicWeights.push_back(weight);
796 m_elementNames.push_back(symbol);
797 m_atomicNumbers.push_back(atomic_number);
798 m_entropy298.push_back(entropy298);
799 if (symbol == "E") {
801 } else {
802 m_elem_type.push_back(elem_type);
803 }
805
806 // Update species compositions
807 if (m_kk) {
808 vector<double> old(m_speciesComp);
809 m_speciesComp.resize(m_kk*m_mm, 0.0);
810 for (size_t k = 0; k < m_kk; k++) {
811 size_t m_old = m_mm - 1;
812 for (size_t m = 0; m < m_old; m++) {
813 m_speciesComp[k * m_mm + m] = old[k * (m_old) + m];
814 }
815 m_speciesComp[k * (m_mm) + (m_mm-1)] = 0.0;
816 }
817 }
818
819 return m_mm-1;
820}
821
822bool Phase::addSpecies(shared_ptr<Species> spec) {
823 if (m_species.find(spec->name) != m_species.end()) {
824 throw CanteraError("Phase::addSpecies",
825 "Phase '{}' already contains a species named '{}'.",
826 m_name, spec->name);
827 }
828
829 vector<double> comp(nElements());
830 for (const auto& [eName, stoich] : spec->composition) {
831 size_t m = elementIndex(eName);
832 if (m == npos) { // Element doesn't exist in this phase
834 case UndefElement::ignore:
835 return false;
836
837 case UndefElement::add:
838 addElement(eName);
839 comp.resize(nElements());
840 m = elementIndex(eName);
841 break;
842
843 case UndefElement::error:
844 default:
845 throw CanteraError("Phase::addSpecies",
846 "Species '{}' contains an undefined element '{}'.",
847 spec->name, eName);
848 }
849 }
850 comp[m] = stoich;
851 }
852
853 size_t ne = nElements();
854 const vector<double>& aw = atomicWeights();
855 if (spec->charge != 0.0) {
856 size_t eindex = elementIndex("E");
857 if (eindex != npos) {
858 double ecomp = comp[eindex];
859 if (fabs(spec->charge + ecomp) > 0.001) {
860 if (ecomp != 0.0) {
861 throw CanteraError("Phase::addSpecies",
862 "Input charge and element E compositions differ "
863 "for species " + spec->name);
864 } else {
865 // Just fix up the element E composition based on the input
866 // species charge
867 comp[eindex] = -spec->charge;
868 }
869 }
870 } else {
871 addElement("E", 0.000545, 0, 0.0, CT_ELEM_TYPE_ELECTRONCHARGE);
872 ne = nElements();
873 eindex = elementIndex("E");
874 comp.resize(ne);
875 comp[ne - 1] = - spec->charge;
876 }
877 }
878
879 double wt = 0.0;
880 for (size_t m = 0; m < ne; m++) {
881 wt += comp[m] * aw[m];
882 }
883
884 // Some surface phases may define species representing empty sites
885 // that have zero molecular weight. Give them a very small molecular
886 // weight to avoid dividing by zero.
887 wt = std::max(wt, Tiny);
888
889 spec->setMolecularWeight(wt);
890
891 m_molwts.push_back(wt);
892 m_rmolwts.push_back(1.0/wt);
893
894 for (size_t m = 0; m < ne; m++) {
895 m_speciesComp.push_back(comp[m]);
896 }
897
898 m_speciesNames.push_back(spec->name);
899 m_species[spec->name] = spec;
900 m_speciesIndices[spec->name] = m_kk;
901 m_speciesCharge.push_back(spec->charge);
902
903 string nLower = toLowerCopy(spec->name);
904 if (m_speciesLower.find(nLower) == m_speciesLower.end()) {
905 m_speciesLower[nLower] = m_kk;
906 } else {
907 m_speciesLower[nLower] = npos;
908 }
909 m_kk++;
910
911 // Ensure that the Phase has a valid mass fraction vector that sums to
912 // one. We will assume that species 0 has a mass fraction of 1.0 and mass
913 // fraction of all other species is 0.0.
914 if (m_kk == 1) {
915 m_y.push_back(1.0);
916 m_ym.push_back(m_rmolwts[0]);
917 m_mmw = 1.0 / m_ym[0];
918 } else {
919 m_y.push_back(0.0);
920 m_ym.push_back(0.0);
921 }
923 return true;
924}
925
926void Phase::modifySpecies(size_t k, shared_ptr<Species> spec)
927{
928 if (speciesName(k) != spec->name) {
929 throw CanteraError("Phase::modifySpecies",
930 "New species name '{}' does not match existing name '{}'",
931 spec->name, speciesName(k));
932 }
933 const shared_ptr<Species>& old = m_species[spec->name];
934 if (spec->composition != old->composition) {
935 throw CanteraError("Phase::modifySpecies",
936 "New composition for '{}' does not match existing composition",
937 spec->name);
938 }
939 m_species[spec->name] = spec;
941}
942
943void Phase::addSpeciesAlias(const string& name, const string& alias)
944{
945 if (speciesIndex(alias) != npos) {
946 throw CanteraError("Phase::addSpeciesAlias",
947 "Invalid alias '{}': species already exists", alias);
948 }
949 size_t k = speciesIndex(name);
950 if (k != npos) {
951 m_speciesIndices[alias] = k;
952 } else {
953 throw CanteraError("Phase::addSpeciesAlias",
954 "Unable to add alias '{}' "
955 "(original species '{}' not found).", alias, name);
956 }
957}
958
959vector<string> Phase::findIsomers(const Composition& compMap) const
960{
961 vector<string> isomerNames;
962
963 for (const auto& [name, species] : m_species) {
964 if (species->composition == compMap) {
965 isomerNames.emplace_back(name);
966 }
967 }
968
969 return isomerNames;
970}
971
972vector<string> Phase::findIsomers(const string& comp) const
973{
974 return findIsomers(parseCompString(comp));
975}
976
977shared_ptr<Species> Phase::species(const string& name) const
978{
979 size_t k = speciesIndex(name);
980 if (k != npos) {
981 return m_species.at(speciesName(k));
982 } else {
983 throw CanteraError("Phase::species",
984 "Unknown species '{}'", name);
985 }
986}
987
988shared_ptr<Species> Phase::species(size_t k) const
989{
991 return m_species.at(m_speciesNames[k]);
992}
993
995 m_undefinedElementBehavior = UndefElement::ignore;
996}
997
999 m_undefinedElementBehavior = UndefElement::add;
1000}
1001
1003 m_undefinedElementBehavior = UndefElement::error;
1004}
1005
1006bool Phase::ready() const
1007{
1008 return (m_kk > 0);
1009}
1010
1012 m_stateNum++;
1013 m_cache.clear();
1014}
1015
1016void Phase::setMolecularWeight(const int k, const double mw)
1017{
1018 m_molwts[k] = mw;
1019 m_rmolwts[k] = 1.0/mw;
1020
1021 transform(m_y.begin(), m_y.end(), m_rmolwts.begin(), m_ym.begin(),
1022 multiplies<double>());
1023 m_mmw = 1.0 / accumulate(m_ym.begin(), m_ym.end(), 0.0);
1024}
1025
1027 m_stateNum++;
1028}
1029
1030vector<double> Phase::getCompositionFromMap(const Composition& comp) const
1031{
1032 vector<double> X(m_kk);
1033 for (const auto& [name, value] : comp) {
1034 size_t loc = speciesIndex(name);
1035 if (loc == npos) {
1036 throw CanteraError("Phase::getCompositionFromMap",
1037 "Unknown species '{}'", name);
1038 }
1039 X[loc] = value;
1040 }
1041 return X;
1042}
1043
1044void Phase::massFractionsToMoleFractions(const double* Y, double* X) const
1045{
1046 double rmmw = 0.0;
1047 for (size_t k = 0; k != m_kk; ++k) {
1048 rmmw += Y[k]/m_molwts[k];
1049 }
1050 if (rmmw == 0.0) {
1051 throw CanteraError("Phase::massFractionsToMoleFractions",
1052 "no input composition given");
1053 }
1054 for (size_t k = 0; k != m_kk; ++k) {
1055 X[k] = Y[k]/(rmmw*m_molwts[k]);
1056 }
1057}
1058
1059void Phase::moleFractionsToMassFractions(const double* X, double* Y) const
1060{
1061 double mmw = dot(X, X+m_kk, m_molwts.data());
1062 if (mmw == 0.0) {
1063 throw CanteraError("Phase::moleFractionsToMassFractions",
1064 "no input composition given");
1065 }
1066 double rmmw = 1.0/mmw;
1067 for (size_t k = 0; k != m_kk; ++k) {
1068 Y[k] = X[k]*m_molwts[k]*rmmw;
1069 }
1070}
1071
1072} // 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:427
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:1535
static AnyMap fromYamlFile(const string &name, const string &parent_name="")
Create an AnyMap from a YAML file.
Definition AnyMap.cpp:1771
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:230
void getCharges(double *charges) const
Copy the vector of species charges into array charges.
Definition Phase.cpp:511
virtual void getConcentrations(double *const c) const
Get the species concentrations (kmol/m^3).
Definition Phase.cpp:595
map< string, size_t > m_speciesLower
Map of lower-case species names to indices.
Definition Phase.h:1015
double massFraction(size_t k) const
Return the mass fraction of a single species.
Definition Phase.cpp:568
virtual double molarDensity() const
Molar density (kmol/m^3).
Definition Phase.cpp:689
void assignDensity(const double density_)
Set the internally stored constant density (kg/m^3) of the phase.
Definition Phase.cpp:718
virtual bool addSpecies(shared_ptr< Species > spec)
Add a Species to this Phase.
Definition Phase.cpp:822
virtual void setMoleFractions(const double *const x)
Set the mole fractions to the specified values.
Definition Phase.cpp:304
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:212
void assertCompressible(const string &setter) const
Ensure that phase is compressible.
Definition Phase.h:907
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range.
Definition Phase.cpp:162
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
Definition Phase.cpp:275
vector< double > m_speciesComp
Atomic composition of the species.
Definition Phase.h:956
ValueCache m_cache
Cached for saved calculations within each ThermoPhase.
Definition Phase.h:927
void getMolecularWeights(vector< double > &weights) const
Copy the vector of molecular weights into vector weights.
Definition Phase.cpp:488
const double * moleFractdivMMW() const
Returns a const pointer to the start of the moleFraction/MW array.
Definition Phase.cpp:561
vector< string > m_speciesNames
Vector of the species names.
Definition Phase.h:1009
void setState_RY(double rho, double *y)
Set the density (kg/m^3) and mass fractions.
Definition Phase.cpp:474
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:245
bool m_caseSensitiveSpecies
Flag determining whether case sensitive species names are enforced.
Definition Phase.h:966
vector< string > m_elementNames
element names
Definition Phase.h:1020
void ignoreUndefinedElements()
Set behavior when adding a species containing undefined elements to just skip the species.
Definition Phase.cpp:994
UndefElement::behavior m_undefinedElementBehavior
Flag determining behavior when adding species with an undefined element.
Definition Phase.h:963
virtual void setMassFractions_NoNorm(const double *const y)
Set the mass fractions to the specified values without normalizing.
Definition Phase.cpp:370
virtual map< string, size_t > nativeState() const
Return a map of properties defining the native state of a substance.
Definition Phase.cpp:182
double chargeDensity() const
Charge density [C/m^3].
Definition Phase.cpp:728
void addUndefinedElements()
Set behavior when adding a species containing undefined elements to add those elements to the phase.
Definition Phase.cpp:998
virtual void setConcentrationsNoNorm(const double *const conc)
Set the concentrations without ignoring negative concentrations.
Definition Phase.cpp:622
vector< int > m_atomicNumbers
element atomic numbers
Definition Phase.h:1019
size_t m_kk
Number of species in the phase.
Definition Phase.h:947
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:926
void setState_TRY(double t, double dens, const double *y)
Set the internally stored temperature (K), density, and mass fractions.
Definition Phase.cpp:419
double m_mmw
mean molecular weight of the mixture (kg kmol-1)
Definition Phase.h:987
double elementalMoleFraction(const size_t m) const
Elemental mole fraction of element m.
Definition Phase.cpp:671
void setState_TR(double t, double rho)
Set the internally stored temperature (K) and density (kg/m^3)
Definition Phase.cpp:437
void setState_TD(double t, double rho)
Set the internally stored temperature (K) and density (kg/m^3)
Definition Phase.cpp:444
vector< double > m_rmolwts
inverse of species molecular weights (kmol kg-1)
Definition Phase.h:1002
double temperature() const
Temperature (K).
Definition Phase.h:662
string speciesSPName(int k) const
Returns the expanded species name of a species, including the phase name This is guaranteed to be uni...
Definition Phase.cpp:176
virtual void setPressure(double p)
Set the internally stored pressure (Pa) at constant temperature and composition.
Definition Phase.h:721
virtual bool isCompressible() const
Return whether phase represents a compressible substance.
Definition Phase.h:271
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:760
void moleFractionsToMassFractions(const double *X, double *Y) const
Converts a mixture composition from mass fractions to mole fractions.
Definition Phase.cpp:1059
virtual void setConcentrations(const double *const conc)
Set the concentrations to the specified values within the phase.
Definition Phase.cpp:600
void setState_RX(double rho, double *x)
Set the density (kg/m^3) and mole fractions.
Definition Phase.cpp:466
virtual void setMolarDensity(const double molarDensity)
Set the internally stored molar density (kmol/m^3) of the phase.
Definition Phase.cpp:694
void saveState(vector< double > &state) const
Save the current internal state of the phase.
Definition Phase.cpp:251
void setState_TX(double t, double *x)
Set the internally stored temperature (K) and mole fractions.
Definition Phase.cpp:450
Composition getMoleFractionsByName(double threshold=0.0) const
Get the mole fractions by name.
Definition Phase.cpp:516
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:589
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:381
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:151
map< string, size_t > m_speciesIndices
Map of species names to indices.
Definition Phase.h:1012
virtual void setDensity(const double density_)
Set the internally stored density (kg/m^3) of the phase.
Definition Phase.cpp:707
Composition getMassFractionsByName(double threshold=0.0) const
Get the mass fractions by name.
Definition Phase.cpp:528
void setState_TRX(double t, double dens, const double *x)
Set the internally stored temperature (K), density, and mole fractions.
Definition Phase.cpp:392
virtual size_t stateSize() const
Return size of vector defining internal state of the phase.
Definition Phase.cpp:243
string nativeMode() const
Return string acronym representing the native state of a Phase.
Definition Phase.cpp:199
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:1030
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:119
vector< double > m_molwts
species molecular weights (kg kmol-1)
Definition Phase.h:1000
virtual vector< string > findIsomers(const Composition &compMap) const
Return a vector with isomers names matching a given composition map.
Definition Phase.cpp:959
const vector< double > & inverseMolecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition Phase.cpp:506
virtual bool isPure() const
Return whether phase represents a pure (single species) substance.
Definition Phase.h:261
vector< double > m_y
Mass fractions of the species.
Definition Phase.h:998
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition Phase.cpp:540
void setMoleFractionsByName(const Composition &xMap)
Set the species mole fractions by name.
Definition Phase.cpp:345
const double * massFractions() const
Return a const pointer to the mass fraction array.
Definition Phase.h:535
vector< int > m_elem_type
Vector of element types.
Definition Phase.h:1021
void setState_TNX(double t, double n, const double *x)
Set the internally stored temperature (K), molar density (kmol/m^3), and mole fractions.
Definition Phase.cpp:401
double sum_xlogx() const
Evaluate .
Definition Phase.cpp:747
string m_name
Name of the phase.
Definition Phase.h:977
const vector< double > & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition Phase.cpp:501
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition Phase.cpp:138
double moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition Phase.cpp:545
double m_dens
Density (kg m-3).
Definition Phase.h:985
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:687
virtual void compositionChanged()
Apply changes to the state which are needed after the composition changes.
Definition Phase.cpp:1026
vector< double > m_atomicWeights
element atomic weights (kg kmol-1)
Definition Phase.h:1018
void checkSpeciesArraySize(size_t kk) const
Check that an array size is at least nSpecies().
Definition Phase.cpp:169
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:943
virtual void setMolesNoTruncate(const double *const N)
Set the state of the object with moles in [kmol].
Definition Phase.cpp:641
virtual void setTemperature(double temp)
Set the internally stored temperature of the phase (K).
Definition Phase.h:728
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:1016
double mean_X(const double *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition Phase.cpp:737
vector< double > m_entropy298
Entropy at 298.15 K and 1 bar of stable state pure elements (J kmol-1)
Definition Phase.h:1024
vector< double > m_ym
m_ym[k] = mole fraction of species k divided by the mean molecular weight of mixture.
Definition Phase.h:991
virtual void setMassFractions(const double *const y)
Set the mass fractions to the specified values and normalize them.
Definition Phase.cpp:356
const vector< string > & speciesNames() const
Return a const reference to the vector of species names.
Definition Phase.cpp:157
virtual bool ready() const
Returns a bool indicating whether the object is ready for use.
Definition Phase.cpp:1006
double molecularWeight(size_t k) const
Molecular weight of species k.
Definition Phase.cpp:482
double elementalMassFraction(const size_t m) const
Elemental mass fraction of element m.
Definition Phase.cpp:660
shared_ptr< Species > species(const string &name) const
Return the Species object for the named species.
Definition Phase.cpp:977
virtual double molarVolume() const
Molar volume (m^3/kmol).
Definition Phase.cpp:702
virtual void invalidateCache()
Invalidate any cached values which are normally updated only when a change in state is detected.
Definition Phase.cpp:1011
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:584
void setState_TY(double t, double *y)
Set the internally stored temperature (K) and mass fractions.
Definition Phase.cpp:458
void getAtoms(size_t k, double *atomArray) const
Get a vector containing the atomic composition of species k.
Definition Phase.cpp:110
int m_stateNum
State Change variable.
Definition Phase.h:1006
void throwUndefinedElements()
Set the behavior when adding a species containing undefined elements to throw an exception.
Definition Phase.cpp:1002
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:1017
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition Phase.h:680
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:638
void massFractionsToMoleFractions(const double *Y, double *X) const
Converts a mixture composition from mole fractions to mass fractions.
Definition Phase.cpp:1044
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:336
vector< double > m_speciesCharge
Vector of species charges. length m_kk.
Definition Phase.h:958
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:756
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:564
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:195
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
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Definition AnyMap.cpp:1926
map< string, double > Composition
Map from string names to doubles.
Definition ct_defs.h:184
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...