Cantera  3.1.0a1
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 
9 #include "cantera/thermo/Phase.h"
10 #include "cantera/base/utilities.h"
12 #include "cantera/thermo/Species.h"
14 
15 using namespace std;
16 
17 namespace Cantera
18 {
19 
20 string Phase::name() const
21 {
22  return m_name;
23 }
24 
25 void Phase::setName(const string& name)
26 {
27  m_name = name;
28 }
29 
30 size_t Phase::nElements() const
31 {
32  return m_mm;
33 }
34 
35 void 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 
42 void Phase::checkElementArraySize(size_t mm) const
43 {
44  if (m_mm > mm) {
45  throw ArraySizeError("Phase::checkElementArraySize", mm, m_mm);
46  }
47 }
48 
49 string Phase::elementName(size_t m) const
50 {
51  checkElementIndex(m);
52  return m_elementNames[m];
53 }
54 
55 size_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 
65 const vector<string>& Phase::elementNames() const
66 {
67  return m_elementNames;
68 }
69 
70 double Phase::atomicWeight(size_t m) const
71 {
72  return m_atomicWeights[m];
73 }
74 
75 double Phase::entropyElement298(size_t m) const
76 {
77  checkElementIndex(m);
78  return m_entropy298[m];
79 }
80 
81 const vector<double>& Phase::atomicWeights() const
82 {
83  return m_atomicWeights;
84 }
85 
86 int Phase::atomicNumber(size_t m) const
87 {
88  return m_atomicNumbers[m];
89 }
90 
91 int Phase::elementType(size_t m) const
92 {
93  return m_elem_type[m];
94 }
95 
96 int 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 
103 double Phase::nAtoms(size_t k, size_t m) const
104 {
105  checkElementIndex(m);
106  checkSpeciesIndex(k);
107  return m_speciesComp[m_mm * k + m];
108 }
109 
110 size_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 
129 size_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 
142 string Phase::speciesName(size_t k) const
143 {
144  checkSpeciesIndex(k);
145  return m_speciesNames[k];
146 }
147 
148 const vector<string>& Phase::speciesNames() const
149 {
150  return m_speciesNames;
151 }
152 
153 void Phase::checkSpeciesIndex(size_t k) const
154 {
155  if (k >= m_kk) {
156  throw IndexError("Phase::checkSpeciesIndex", "species", k, m_kk-1);
157  }
158 }
159 
160 void Phase::checkSpeciesArraySize(size_t kk) const
161 {
162  if (m_kk > kk) {
163  throw ArraySizeError("Phase::checkSpeciesArraySize", kk, m_kk);
164  }
165 }
166 
167 map<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 
184 string 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 
197 vector<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 
215 vector<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 
228 size_t Phase::stateSize() const {
229  if (isPure()) {
230  return 2;
231  } else {
232  return nSpecies() + 2;
233  }
234 }
235 
236 void Phase::saveState(vector<double>& state) const
237 {
238  state.resize(stateSize());
239  saveState(state.size(), &state[0]);
240 }
241 
242 void Phase::saveState(size_t lenstate, double* state) const
243 {
244  auto native = nativeState();
245 
246  // function assumes default definition of nativeState
247  state[native.at("T")] = temperature();
248  if (isCompressible()) {
249  state[native.at("D")] = density();
250  } else {
251  state[native.at("P")] = pressure();
252  }
253  if (native.count("X")) {
254  getMoleFractions(state + native["X"]);
255  } else if (native.count("Y")) {
256  getMassFractions(state + native["Y"]);
257  }
258 }
259 
260 void Phase::restoreState(const vector<double>& state)
261 {
262  restoreState(state.size(),&state[0]);
263 }
264 
265 void Phase::restoreState(size_t lenstate, const double* state)
266 {
267  size_t ls = stateSize();
268  if (lenstate < ls) {
269  throw ArraySizeError("Phase::restoreState",
270  lenstate, ls);
271  }
272 
273  auto native = nativeState();
274  setTemperature(state[native.at("T")]);
275  if (isCompressible()) {
276  setDensity(state[native.at("D")]);
277  } else {
278  setPressure(state[native.at("P")]);
279  }
280 
281  if (native.count("X")) {
282  setMoleFractions_NoNorm(state + native["X"]);
283  } else if (native.count("Y")) {
284  setMassFractions_NoNorm(state + native["Y"]);
285  }
286  compositionChanged();
287 }
288 
289 void Phase::setMoleFractions(const double* const x)
290 {
291  // Use m_y as a temporary work vector for the non-negative mole fractions
292  double norm = 0.0;
293  // sum is calculated below as the unnormalized molecular weight
294  double sum = 0;
295  for (size_t k = 0; k < m_kk; k++) {
296  double xk = std::max(x[k], 0.0); // Ignore negative mole fractions
297  m_y[k] = xk;
298  norm += xk;
299  sum += m_molwts[k] * xk;
300  }
301 
302  // Set m_ym to the normalized mole fractions divided by the normalized mean
303  // molecular weight:
304  // m_ym_k = X_k / (sum_k X_k M_k)
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;
308  }
309 
310  // Now set m_y to the normalized mass fractions:
311  // m_y = X_k M_k / (sum_k X_k M_k)
312  for (size_t k=0; k < m_kk; k++) {
313  m_y[k] = m_ym[k] * m_molwts[k];
314  }
315 
316  // Calculate the normalized molecular weight
317  m_mmw = sum/norm;
318  compositionChanged();
319 }
320 
321 void Phase::setMoleFractions_NoNorm(const double* const x)
322 {
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();
328 }
329 
330 void Phase::setMoleFractionsByName(const Composition& xMap)
331 {
332  vector<double> mf = getCompositionFromMap(xMap);
333  setMoleFractions(mf.data());
334 }
335 
336 void Phase::setMoleFractionsByName(const string& x)
337 {
338  setMoleFractionsByName(parseCompString(x));
339 }
340 
341 void Phase::setMassFractions(const double* const y)
342 {
343  for (size_t k = 0; k < m_kk; k++) {
344  m_y[k] = std::max(y[k], 0.0); // Ignore negative mass fractions
345  }
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);
348 
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();
353 }
354 
355 void Phase::setMassFractions_NoNorm(const double* const y)
356 {
357  double sum = 0.0;
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);
362  m_mmw = 1.0/sum;
363  compositionChanged();
364 }
365 
366 void Phase::setMassFractionsByName(const Composition& yMap)
367 {
368  vector<double> mf = getCompositionFromMap(yMap);
369  setMassFractions(mf.data());
370 }
371 
372 void Phase::setMassFractionsByName(const string& y)
373 {
374  setMassFractionsByName(parseCompString(y));
375 }
376 
377 void Phase::setState_TD(double t, double rho)
378 {
379  setTemperature(t);
380  setDensity(rho);
381 }
382 
383 double Phase::molecularWeight(size_t k) const
384 {
385  checkSpeciesIndex(k);
386  return m_molwts[k];
387 }
388 
389 void Phase::getMolecularWeights(double* weights) const
390 {
391  const vector<double>& mw = molecularWeights();
392  copy(mw.begin(), mw.end(), weights);
393 }
394 
395 const vector<double>& Phase::molecularWeights() const
396 {
397  return m_molwts;
398 }
399 
400 const vector<double>& Phase::inverseMolecularWeights() const
401 {
402  return m_rmolwts;
403 }
404 
405 void Phase::getCharges(double* charges) const
406 {
407  copy(m_speciesCharge.begin(), m_speciesCharge.end(), charges);
408 }
409 
410 Composition Phase::getMoleFractionsByName(double threshold) const
411 {
412  Composition comp;
413  for (size_t k = 0; k < m_kk; k++) {
414  double x = moleFraction(k);
415  if (x > threshold) {
416  comp[speciesName(k)] = x;
417  }
418  }
419  return comp;
420 }
421 
422 Composition Phase::getMassFractionsByName(double threshold) const
423 {
424  Composition comp;
425  for (size_t k = 0; k < m_kk; k++) {
426  double x = massFraction(k);
427  if (x > threshold) {
428  comp[speciesName(k)] = x;
429  }
430  }
431  return comp;
432 }
433 
434 void Phase::getMoleFractions(double* const x) const
435 {
436  scale(m_ym.begin(), m_ym.end(), x, m_mmw);
437 }
438 
439 double Phase::moleFraction(size_t k) const
440 {
441  checkSpeciesIndex(k);
442  return m_ym[k] * m_mmw;
443 }
444 
445 double Phase::moleFraction(const string& nameSpec) const
446 {
447  size_t iloc = speciesIndex(nameSpec);
448  if (iloc != npos) {
449  return moleFraction(iloc);
450  } else {
451  return 0.0;
452  }
453 }
454 
455 double Phase::massFraction(size_t k) const
456 {
457  checkSpeciesIndex(k);
458  return m_y[k];
459 }
460 
461 double Phase::massFraction(const string& nameSpec) const
462 {
463  size_t iloc = speciesIndex(nameSpec);
464  if (iloc != npos) {
465  return massFractions()[iloc];
466  } else {
467  return 0.0;
468  }
469 }
470 
471 void Phase::getMassFractions(double* const y) const
472 {
473  copy(m_y.begin(), m_y.end(), y);
474 }
475 
476 double Phase::concentration(const size_t k) const
477 {
478  checkSpeciesIndex(k);
479  return m_y[k] * m_dens * m_rmolwts[k];
480 }
481 
482 void Phase::getConcentrations(double* const c) const
483 {
484  scale(m_ym.begin(), m_ym.end(), c, m_dens);
485 }
486 
487 void Phase::setConcentrations(const double* const conc)
488 {
489  assertCompressible("setConcentrations");
490 
491  // Use m_y as temporary storage for non-negative concentrations
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); // Ignore negative concentrations
495  m_y[k] = ck;
496  sum += ck * m_molwts[k];
497  norm += ck;
498  }
499  m_mmw = sum/norm;
500  setDensity(sum);
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]; // m_y is now the mass fraction
505  }
506  compositionChanged();
507 }
508 
509 void Phase::setConcentrationsNoNorm(const double* const conc)
510 {
511  assertCompressible("setConcentrationsNoNorm");
512 
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];
516  norm += conc[k];
517  }
518  m_mmw = sum/norm;
519  setDensity(sum);
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];
524  }
525  compositionChanged();
526 }
527 
528 void Phase::setMolesNoTruncate(const double* const N)
529 {
530  // get total moles
531  copy(N, N + m_kk, m_ym.begin());
532  double totalMoles = accumulate(m_ym.begin(), m_ym.end(), 0.0);
533  // get total mass
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);
537  // mean molecular weight
538  m_mmw = totalMass/totalMoles;
539  // mass fractions
540  scale(m_y.begin(), m_y.end(), m_y.begin(), 1/totalMass);
541  // moles fractions/m_mmw
542  scale(m_ym.begin(), m_ym.end(), m_ym.begin(), 1/(m_mmw * totalMoles));
543  // composition has changed
544  compositionChanged();
545 }
546 
547 double Phase::elementalMassFraction(const size_t m) const
548 {
549  checkElementIndex(m);
550  double Z_m = 0.0;
551  for (size_t k = 0; k != m_kk; ++k) {
552  Z_m += nAtoms(k, m) * atomicWeight(m) / molecularWeight(k)
553  * massFraction(k);
554  }
555  return Z_m;
556 }
557 
558 double Phase::elementalMoleFraction(const size_t m) const
559 {
560  checkElementIndex(m);
561  double denom = 0;
562  for (size_t k = 0; k < m_kk; k++) {
563  double atoms = 0;
564  for (size_t j = 0; j < nElements(); j++) {
565  atoms += nAtoms(k, j);
566  }
567  denom += atoms * moleFraction(k);
568  }
569  double numerator = 0.0;
570  for (size_t k = 0; k != m_kk; ++k) {
571  numerator += nAtoms(k, m) * moleFraction(k);
572  }
573  return numerator / denom;
574 }
575 
576 double Phase::molarDensity() const
577 {
578  return density()/meanMolecularWeight();
579 }
580 
581 double Phase::molarVolume() const
582 {
583  return 1.0/molarDensity();
584 }
585 
586 void Phase::setDensity(const double density_)
587 {
588  assertCompressible("setDensity");
589  if (density_ > 0.0) {
590  m_dens = density_;
591  } else {
592  throw CanteraError("Phase::setDensity",
593  "density must be positive. density = {}", density_);
594  }
595 }
596 
597 void Phase::assignDensity(const double density_)
598 {
599  if (density_ > 0.0) {
600  m_dens = density_;
601  } else {
602  throw CanteraError("Phase::assignDensity",
603  "density must be positive. density = {}", density_);
604  }
605 }
606 
607 double Phase::chargeDensity() const
608 {
609  double cdens = 0.0;
610  for (size_t k = 0; k < m_kk; k++) {
611  cdens += charge(k)*moleFraction(k);
612  }
613  return cdens * Faraday;
614 }
615 
616 double Phase::mean_X(const double* const Q) const
617 {
618  return m_mmw*std::inner_product(m_ym.begin(), m_ym.end(), Q, 0.0);
619 }
620 
621 double Phase::mean_X(const vector<double>& Q) const
622 {
623  return m_mmw*std::inner_product(m_ym.begin(), m_ym.end(), Q.begin(), 0.0);
624 }
625 
626 double Phase::sum_xlogx() const
627 {
628  double sumxlogx = 0;
629  for (size_t k = 0; k < m_kk; k++) {
630  sumxlogx += m_ym[k] * std::log(std::max(m_ym[k], SmallNumber));
631  }
632  return m_mmw * sumxlogx + std::log(m_mmw);
633 }
634 
635 size_t Phase::addElement(const string& symbol, double weight, int atomic_number,
636  double entropy298, int elem_type)
637 {
638  // Look up the atomic weight if not given
639  if (weight == 0.0) {
640  try {
641  weight = getElementWeight(symbol);
642  } catch (CanteraError&) {
643  // assume this is just a custom element with zero atomic weight
644  }
645  } else if (weight == -12345.0) {
646  weight = getElementWeight(symbol);
647  }
648 
649  // Try to look up the standard entropy if not given. Fail silently.
650  if (entropy298 == ENTROPY298_UNKNOWN) {
651  try {
652  const static AnyMap db = AnyMap::fromYamlFile(
653  "element-standard-entropies.yaml");
654  const AnyMap& elem = db["elements"].getMapWhere("symbol", symbol);
655  entropy298 = elem.convert("entropy298", "J/kmol/K", ENTROPY298_UNKNOWN);
656  } catch (CanteraError&) {
657  }
658  }
659 
660  // Check for duplicates
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) {
665  throw CanteraError("Phase::addElement",
666  "Duplicate elements ({}) have different weights", symbol);
667  } else {
668  // Ignore attempt to add duplicate element with the same weight
669  return m;
670  }
671  }
672 
673  // Add the new element
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);
678  if (symbol == "E") {
679  m_elem_type.push_back(CT_ELEM_TYPE_ELECTRONCHARGE);
680  } else {
681  m_elem_type.push_back(elem_type);
682  }
683  m_mm++;
684 
685  // Update species compositions
686  if (m_kk) {
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];
693  }
694  m_speciesComp[k * (m_mm) + (m_mm-1)] = 0.0;
695  }
696  }
697 
698  return m_mm-1;
699 }
700 
701 bool Phase::addSpecies(shared_ptr<Species> spec) {
702  if (m_species.find(spec->name) != m_species.end()) {
703  throw CanteraError("Phase::addSpecies",
704  "Phase '{}' already contains a species named '{}'.",
705  m_name, spec->name);
706  }
707 
708  vector<double> comp(nElements());
709  for (const auto& [eName, stoich] : spec->composition) {
710  size_t m = elementIndex(eName);
711  if (m == npos) { // Element doesn't exist in this phase
712  switch (m_undefinedElementBehavior) {
713  case UndefElement::ignore:
714  return false;
715 
716  case UndefElement::add:
717  addElement(eName);
718  comp.resize(nElements());
719  m = elementIndex(eName);
720  break;
721 
722  case UndefElement::error:
723  default:
724  throw CanteraError("Phase::addSpecies",
725  "Species '{}' contains an undefined element '{}'.",
726  spec->name, eName);
727  }
728  }
729  comp[m] = stoich;
730  }
731 
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) {
739  if (ecomp != 0.0) {
740  throw CanteraError("Phase::addSpecies",
741  "Input charge and element E compositions differ "
742  "for species " + spec->name);
743  } else {
744  // Just fix up the element E composition based on the input
745  // species charge
746  comp[eindex] = -spec->charge;
747  }
748  }
749  } else {
750  addElement("E", 0.000545, 0, 0.0, CT_ELEM_TYPE_ELECTRONCHARGE);
751  ne = nElements();
752  eindex = elementIndex("E");
753  comp.resize(ne);
754  comp[ne - 1] = - spec->charge;
755  }
756  }
757 
758  double wt = 0.0;
759  for (size_t m = 0; m < ne; m++) {
760  wt += comp[m] * aw[m];
761  }
762 
763  // Some surface phases may define species representing empty sites
764  // that have zero molecular weight. Give them a very small molecular
765  // weight to avoid dividing by zero.
766  wt = std::max(wt, Tiny);
767 
768  spec->setMolecularWeight(wt);
769 
770  m_molwts.push_back(wt);
771  m_rmolwts.push_back(1.0/wt);
772 
773  for (size_t m = 0; m < ne; m++) {
774  m_speciesComp.push_back(comp[m]);
775  }
776 
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);
781 
782  string nLower = toLowerCopy(spec->name);
783  if (m_speciesLower.find(nLower) == m_speciesLower.end()) {
784  m_speciesLower[nLower] = m_kk;
785  } else {
786  m_speciesLower[nLower] = npos;
787  }
788  m_kk++;
789 
790  // Ensure that the Phase has a valid mass fraction vector that sums to
791  // one. We will assume that species 0 has a mass fraction of 1.0 and mass
792  // fraction of all other species is 0.0.
793  if (m_kk == 1) {
794  m_y.push_back(1.0);
795  m_ym.push_back(m_rmolwts[0]);
796  m_mmw = 1.0 / m_ym[0];
797  } else {
798  m_y.push_back(0.0);
799  m_ym.push_back(0.0);
800  }
801  invalidateCache();
802  return true;
803 }
804 
805 void Phase::modifySpecies(size_t k, shared_ptr<Species> spec)
806 {
807  if (speciesName(k) != spec->name) {
808  throw CanteraError("Phase::modifySpecies",
809  "New species name '{}' does not match existing name '{}'",
810  spec->name, speciesName(k));
811  }
812  const shared_ptr<Species>& old = m_species[spec->name];
813  if (spec->composition != old->composition) {
814  throw CanteraError("Phase::modifySpecies",
815  "New composition for '{}' does not match existing composition",
816  spec->name);
817  }
818  m_species[spec->name] = spec;
819  invalidateCache();
820 }
821 
822 void Phase::addSpeciesAlias(const string& name, const string& alias)
823 {
824  if (speciesIndex(alias) != npos) {
825  throw CanteraError("Phase::addSpeciesAlias",
826  "Invalid alias '{}': species already exists", alias);
827  }
828  size_t k = speciesIndex(name);
829  if (k != npos) {
830  m_speciesIndices[alias] = k;
831  } else {
832  throw CanteraError("Phase::addSpeciesAlias",
833  "Unable to add alias '{}' "
834  "(original species '{}' not found).", alias, name);
835  }
836 }
837 
838 vector<string> Phase::findIsomers(const Composition& compMap) const
839 {
840  vector<string> isomerNames;
841 
842  for (const auto& [name, species] : m_species) {
843  if (species->composition == compMap) {
844  isomerNames.emplace_back(name);
845  }
846  }
847 
848  return isomerNames;
849 }
850 
851 vector<string> Phase::findIsomers(const string& comp) const
852 {
853  return findIsomers(parseCompString(comp));
854 }
855 
856 shared_ptr<Species> Phase::species(const string& name) const
857 {
858  size_t k = speciesIndex(name);
859  if (k != npos) {
860  return m_species.at(speciesName(k));
861  } else {
862  throw CanteraError("Phase::species",
863  "Unknown species '{}'", name);
864  }
865 }
866 
867 shared_ptr<Species> Phase::species(size_t k) const
868 {
869  checkSpeciesIndex(k);
870  return m_species.at(m_speciesNames[k]);
871 }
872 
873 void Phase::ignoreUndefinedElements() {
874  m_undefinedElementBehavior = UndefElement::ignore;
875 }
876 
877 void Phase::addUndefinedElements() {
878  m_undefinedElementBehavior = UndefElement::add;
879 }
880 
881 void Phase::throwUndefinedElements() {
882  m_undefinedElementBehavior = UndefElement::error;
883 }
884 
885 bool Phase::ready() const
886 {
887  return (m_kk > 0);
888 }
889 
890 void Phase::invalidateCache() {
891  m_stateNum++;
892  m_cache.clear();
893 }
894 
895 void Phase::setMolecularWeight(const int k, const double mw)
896 {
897  m_molwts[k] = mw;
898  m_rmolwts[k] = 1.0/mw;
899 
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);
903 }
904 
905 void Phase::compositionChanged() {
906  m_stateNum++;
907 }
908 
909 vector<double> Phase::getCompositionFromMap(const Composition& comp) const
910 {
911  vector<double> X(m_kk);
912  for (const auto& [name, value] : comp) {
913  size_t loc = speciesIndex(name);
914  if (loc == npos) {
915  throw CanteraError("Phase::getCompositionFromMap",
916  "Unknown species '{}'", name);
917  }
918  X[loc] = value;
919  }
920  return X;
921 }
922 
923 void Phase::massFractionsToMoleFractions(const double* Y, double* X) const
924 {
925  double rmmw = 0.0;
926  for (size_t k = 0; k != m_kk; ++k) {
927  rmmw += Y[k]/m_molwts[k];
928  }
929  if (rmmw == 0.0) {
930  throw CanteraError("Phase::massFractionsToMoleFractions",
931  "no input composition given");
932  }
933  for (size_t k = 0; k != m_kk; ++k) {
934  X[k] = Y[k]/(rmmw*m_molwts[k]);
935  }
936 }
937 
938 void Phase::moleFractionsToMassFractions(const double* X, double* Y) const
939 {
940  double mmw = dot(X, X+m_kk, m_molwts.data());
941  if (mmw == 0.0) {
942  throw CanteraError("Phase::moleFractionsToMassFractions",
943  "no input composition given");
944  }
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;
948  }
949 }
950 
951 } // 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
Array size error.
Definition: ctexceptions.h:135
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
An array index is out of range.
Definition: ctexceptions.h:165
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.
Definition: stringUtils.cpp:58
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: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 vector< string > & elementNames()
Get a vector of the names of the elements defined in Cantera.
Definition: Elements.cpp:227
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...