Cantera  3.1.0a1
ThermoPhase.cpp
Go to the documentation of this file.
1 /**
2  * @file ThermoPhase.cpp
3  * Definition file for class ThermoPhase, the base class for phases with
4  * thermodynamic properties
5  * (see class @link Cantera::ThermoPhase ThermoPhase@endlink).
6  */
7 
8 // This file is part of Cantera. See License.txt in the top-level directory or
9 // at https://cantera.org/license.txt for license and copyright information.
10 
14 #include "cantera/thermo/Species.h"
18 
19 #include <iomanip>
20 #include <fstream>
21 #include <numeric>
22 
23 using namespace std;
24 
25 namespace Cantera
26 {
27 
28 void ThermoPhase::resetHf298(size_t k) {
29  if (k != npos) {
30  m_spthermo.resetHf298(k);
31  } else {
32  for (size_t k = 0; k < nSpecies(); k++) {
33  m_spthermo.resetHf298(k);
34  }
35  }
36  invalidateCache();
37 }
38 
39 int ThermoPhase::activityConvention() const
40 {
41  return cAC_CONVENTION_MOLAR;
42 }
43 
44 int ThermoPhase::standardStateConvention() const
45 {
46  return m_ssConvention;
47 }
48 
49 Units ThermoPhase::standardConcentrationUnits() const
50 {
51  // kmol/m^3 for bulk phases, kmol/m^2 for surface phases, etc.
52  return Units(1.0, 0, -static_cast<double>(nDim()), 0, 0, 0, 1);
53 }
54 
55 double ThermoPhase::logStandardConc(size_t k) const
56 {
57  return log(standardConcentration(k));
58 }
59 
60 void ThermoPhase::getActivities(double* a) const
61 {
62  getActivityConcentrations(a);
63  for (size_t k = 0; k < nSpecies(); k++) {
64  a[k] /= standardConcentration(k);
65  }
66 }
67 
68 void ThermoPhase::getLnActivityCoefficients(double* lnac) const
69 {
70  getActivityCoefficients(lnac);
71  for (size_t k = 0; k < m_kk; k++) {
72  lnac[k] = std::log(lnac[k]);
73  }
74 }
75 
76 void ThermoPhase::getElectrochemPotentials(double* mu) const
77 {
78  getChemPotentials(mu);
79  double ve = Faraday * electricPotential();
80  for (size_t k = 0; k < m_kk; k++) {
81  mu[k] += ve*charge(k);
82  }
83 }
84 
85 void ThermoPhase::setState_TPX(double t, double p, const double* x)
86 {
87  setMoleFractions(x);
88  setState_TP(t,p);
89 }
90 
91 void ThermoPhase::setState_TPX(double t, double p, const Composition& x)
92 {
93  setMoleFractionsByName(x);
94  setState_TP(t,p);
95 }
96 
97 void ThermoPhase::setState_TPX(double t, double p, const string& x)
98 {
99  setMoleFractionsByName(x);
100  setState_TP(t,p);
101 }
102 
103 void ThermoPhase::setState_TPY(double t, double p, const double* y)
104 {
105  setMassFractions(y);
106  setState_TP(t,p);
107 }
108 
109 void ThermoPhase::setState_TPY(double t, double p, const Composition& y)
110 {
111  setMassFractionsByName(y);
112  setState_TP(t,p);
113 }
114 
115 void ThermoPhase::setState_TPY(double t, double p, const string& y)
116 {
117  setMassFractionsByName(y);
118  setState_TP(t,p);
119 }
120 
121 void ThermoPhase::setState_TP(double t, double p)
122 {
123  double tsave = temperature();
124  double dsave = density();
125  try {
126  setTemperature(t);
127  setPressure(p);
128  } catch (CanteraError&) {
129  setState_TD(tsave, dsave);
130  throw;
131  }
132 }
133 
134 void ThermoPhase::setState_HP(double Htarget, double p, double rtol)
135 {
136  setState_HPorUV(Htarget, p, rtol, false);
137 }
138 
139 void ThermoPhase::setState_UV(double u, double v, double rtol)
140 {
141  assertCompressible("setState_UV");
142  setState_HPorUV(u, v, rtol, true);
143 }
144 
145 void ThermoPhase::setState(const AnyMap& input_state)
146 {
147  AnyMap state = input_state;
148 
149  // Remap allowable synonyms
150  if (state.hasKey("mass-fractions")) {
151  state["Y"] = state["mass-fractions"];
152  state.erase("mass-fractions");
153  }
154  if (state.hasKey("mole-fractions")) {
155  state["X"] = state["mole-fractions"];
156  state.erase("mole-fractions");
157  }
158  if (state.hasKey("temperature")) {
159  state["T"] = state["temperature"];
160  }
161  if (state.hasKey("pressure")) {
162  state["P"] = state["pressure"];
163  }
164  if (state.hasKey("enthalpy")) {
165  state["H"] = state["enthalpy"];
166  }
167  if (state.hasKey("int-energy")) {
168  state["U"] = state["int-energy"];
169  }
170  if (state.hasKey("internal-energy")) {
171  state["U"] = state["internal-energy"];
172  }
173  if (state.hasKey("specific-volume")) {
174  state["V"] = state["specific-volume"];
175  }
176  if (state.hasKey("entropy")) {
177  state["S"] = state["entropy"];
178  }
179  if (state.hasKey("density")) {
180  state["D"] = state["density"];
181  }
182  if (state.hasKey("vapor-fraction")) {
183  state["Q"] = state["vapor-fraction"];
184  }
185 
186  // Set composition
187  if (state.hasKey("X")) {
188  if (state["X"].is<string>()) {
189  setMoleFractionsByName(state["X"].asString());
190  } else {
191  setMoleFractionsByName(state["X"].asMap<double>());
192  }
193  state.erase("X");
194  } else if (state.hasKey("Y")) {
195  if (state["Y"].is<string>()) {
196  setMassFractionsByName(state["Y"].asString());
197  } else {
198  setMassFractionsByName(state["Y"].asMap<double>());
199  }
200  state.erase("Y");
201  }
202  // set thermodynamic state using whichever property set is found
203  if (state.size() == 0) {
204  setState_TP(298.15, OneAtm);
205  } else if (state.hasKey("T") && state.hasKey("P")) {
206  double T = state.convert("T", "K");
207  double P = state.convert("P", "Pa");
208  if (state.hasKey("Q")) {
209  setState_TPQ(T, P, state["Q"].asDouble());
210  } else {
211  setState_TP(T, P);
212  }
213  } else if (state.hasKey("T") && state.hasKey("D")) {
214  setState_TD(state.convert("T", "K"), state.convert("D", "kg/m^3"));
215  } else if (state.hasKey("T") && state.hasKey("V")) {
216  setState_TV(state.convert("T", "K"), state.convert("V", "m^3/kg"));
217  } else if (state.hasKey("H") && state.hasKey("P")) {
218  setState_HP(state.convert("H", "J/kg"), state.convert("P", "Pa"));
219  } else if (state.hasKey("U") && state.hasKey("V")) {
220  setState_UV(state.convert("U", "J/kg"), state.convert("V", "m^3/kg"));
221  } else if (state.hasKey("S") && state.hasKey("P")) {
222  setState_SP(state.convert("S", "J/kg/K"), state.convert("P", "Pa"));
223  } else if (state.hasKey("S") && state.hasKey("V")) {
224  setState_SV(state.convert("S", "J/kg/K"), state.convert("V", "m^3/kg"));
225  } else if (state.hasKey("S") && state.hasKey("T")) {
226  setState_ST(state.convert("S", "J/kg/K"), state.convert("T", "K"));
227  } else if (state.hasKey("P") && state.hasKey("V")) {
228  setState_PV(state.convert("P", "Pa"), state.convert("V", "m^3/kg"));
229  } else if (state.hasKey("U") && state.hasKey("P")) {
230  setState_UP(state.convert("U", "J/kg"), state.convert("P", "Pa"));
231  } else if (state.hasKey("V") && state.hasKey("H")) {
232  setState_VH(state.convert("V", "m^3/kg"), state.convert("H", "J/kg"));
233  } else if (state.hasKey("T") && state.hasKey("H")) {
234  setState_TH(state.convert("T", "K"), state.convert("H", "J/kg"));
235  } else if (state.hasKey("S") && state.hasKey("H")) {
236  setState_SH(state.convert("S", "J/kg/K"), state.convert("H", "J/kg"));
237  } else if (state.hasKey("D") && state.hasKey("P")) {
238  setState_DP(state.convert("D", "kg/m^3"), state.convert("P", "Pa"));
239  } else if (state.hasKey("P") && state.hasKey("Q")) {
240  setState_Psat(state.convert("P", "Pa"), state["Q"].asDouble());
241  } else if (state.hasKey("T") && state.hasKey("Q")) {
242  setState_Tsat(state.convert("T", "K"), state["Q"].asDouble());
243  } else if (state.hasKey("T")) {
244  setState_TP(state.convert("T", "K"), OneAtm);
245  } else if (state.hasKey("P")) {
246  setState_TP(298.15, state.convert("P", "Pa"));
247  } else {
248  throw CanteraError("ThermoPhase::setState",
249  "'state' did not specify a recognized set of properties.\n"
250  "Keys provided were: {}", input_state.keys_str());
251  }
252 }
253 
254 void ThermoPhase::setState_conditional_TP(double t, double p, bool set_p)
255 {
256  setTemperature(t);
257  if (set_p) {
258  setPressure(p);
259  }
260 }
261 
262 void ThermoPhase::setState_HPorUV(double Htarget, double p,
263  double rtol, bool doUV)
264 {
265  double dt;
266  double v = 0.0;
267 
268  // Assign the specific volume or pressure and make sure it's positive
269  if (doUV) {
270  double v = p;
271  if (v < 1.0E-300) {
272  throw CanteraError("ThermoPhase::setState_HPorUV (UV)",
273  "Input specific volume is too small or negative. v = {}", v);
274  }
275  setDensity(1.0/v);
276  } else {
277  if (p < 1.0E-300) {
278  throw CanteraError("ThermoPhase::setState_HPorUV (HP)",
279  "Input pressure is too small or negative. p = {}", p);
280  }
281  setPressure(p);
282  }
283  double Tmax = maxTemp() + 0.1;
284  double Tmin = minTemp() - 0.1;
285 
286  // Make sure we are within the temperature bounds at the start
287  // of the iteration
288  double Tnew = temperature();
289  double Tinit = Tnew;
290  if (Tnew > Tmax) {
291  Tnew = Tmax - 1.0;
292  } else if (Tnew < Tmin) {
293  Tnew = Tmin + 1.0;
294  }
295  if (Tnew != Tinit) {
296  setState_conditional_TP(Tnew, p, !doUV);
297  }
298 
299  double Hnew = (doUV) ? intEnergy_mass() : enthalpy_mass();
300  double Cpnew = (doUV) ? cv_mass() : cp_mass();
301  double Htop = Hnew;
302  double Ttop = Tnew;
303  double Hbot = Hnew;
304  double Tbot = Tnew;
305 
306  bool ignoreBounds = false;
307  // Unstable phases are those for which cp < 0.0. These are possible for
308  // cases where we have passed the spinodal curve.
309  bool unstablePhase = false;
310  // Counter indicating the last temperature point where the
311  // phase was unstable
312  double Tunstable = -1.0;
313  bool unstablePhaseNew = false;
314 
315  // Newton iteration
316  for (int n = 0; n < 500; n++) {
317  double Told = Tnew;
318  double Hold = Hnew;
319  double cpd = Cpnew;
320  if (cpd < 0.0) {
321  unstablePhase = true;
322  Tunstable = Tnew;
323  }
324  // limit step size to 100 K
325  dt = clip((Htarget - Hold)/cpd, -100.0, 100.0);
326 
327  // Calculate the new T
328  Tnew = Told + dt;
329 
330  // Limit the step size so that we are convergent This is the step that
331  // makes it different from a Newton's algorithm
332  if ((dt > 0.0 && unstablePhase) || (dt <= 0.0 && !unstablePhase)) {
333  if (Hbot < Htarget && Tnew < (0.75 * Tbot + 0.25 * Told)) {
334  dt = 0.75 * (Tbot - Told);
335  Tnew = Told + dt;
336  }
337  } else if (Htop > Htarget && Tnew > (0.75 * Ttop + 0.25 * Told)) {
338  dt = 0.75 * (Ttop - Told);
339  Tnew = Told + dt;
340  }
341 
342  // Check Max and Min values
343  if (Tnew > Tmax && !ignoreBounds) {
344  setState_conditional_TP(Tmax, p, !doUV);
345  double Hmax = (doUV) ? intEnergy_mass() : enthalpy_mass();
346  if (Hmax >= Htarget) {
347  if (Htop < Htarget) {
348  Ttop = Tmax;
349  Htop = Hmax;
350  }
351  } else {
352  Tnew = Tmax + 1.0;
353  ignoreBounds = true;
354  }
355  }
356  if (Tnew < Tmin && !ignoreBounds) {
357  setState_conditional_TP(Tmin, p, !doUV);
358  double Hmin = (doUV) ? intEnergy_mass() : enthalpy_mass();
359  if (Hmin <= Htarget) {
360  if (Hbot > Htarget) {
361  Tbot = Tmin;
362  Hbot = Hmin;
363  }
364  } else {
365  Tnew = Tmin - 1.0;
366  ignoreBounds = true;
367  }
368  }
369 
370  // Try to keep phase within its region of stability
371  // -> Could do a lot better if I calculate the
372  // spinodal value of H.
373  for (int its = 0; its < 10; its++) {
374  Tnew = Told + dt;
375  if (Tnew < Told / 3.0) {
376  Tnew = Told / 3.0;
377  dt = -2.0 * Told / 3.0;
378  }
379  setState_conditional_TP(Tnew, p, !doUV);
380  if (doUV) {
381  Hnew = intEnergy_mass();
382  Cpnew = cv_mass();
383  } else {
384  Hnew = enthalpy_mass();
385  Cpnew = cp_mass();
386  }
387  if (Cpnew < 0.0) {
388  unstablePhaseNew = true;
389  Tunstable = Tnew;
390  } else {
391  unstablePhaseNew = false;
392  break;
393  }
394  if (unstablePhase == false && unstablePhaseNew == true) {
395  dt *= 0.25;
396  }
397  }
398 
399  if (Hnew == Htarget) {
400  return;
401  } else if (Hnew > Htarget && (Htop < Htarget || Hnew < Htop)) {
402  Htop = Hnew;
403  Ttop = Tnew;
404  } else if (Hnew < Htarget && (Hbot > Htarget || Hnew > Hbot)) {
405  Hbot = Hnew;
406  Tbot = Tnew;
407  }
408  // Convergence in H
409  double Herr = Htarget - Hnew;
410  double acpd = std::max(fabs(cpd), 1.0E-5);
411  double denom = std::max(fabs(Htarget), acpd * Tnew);
412  double HConvErr = fabs((Herr)/denom);
413  if (HConvErr < rtol || fabs(dt/Tnew) < rtol) {
414  return;
415  }
416  }
417  // We are here when there hasn't been convergence
418 
419  // Formulate a detailed error message, since questions seem to arise often
420  // about the lack of convergence.
421  string ErrString = "No convergence in 500 iterations\n";
422  if (doUV) {
423  ErrString += fmt::format(
424  "\tTarget Internal Energy = {}\n"
425  "\tCurrent Specific Volume = {}\n"
426  "\tStarting Temperature = {}\n"
427  "\tCurrent Temperature = {}\n"
428  "\tCurrent Internal Energy = {}\n"
429  "\tCurrent Delta T = {}\n",
430  Htarget, v, Tinit, Tnew, Hnew, dt);
431  } else {
432  ErrString += fmt::format(
433  "\tTarget Enthalpy = {}\n"
434  "\tCurrent Pressure = {}\n"
435  "\tStarting Temperature = {}\n"
436  "\tCurrent Temperature = {}\n"
437  "\tCurrent Enthalpy = {}\n"
438  "\tCurrent Delta T = {}\n",
439  Htarget, p, Tinit, Tnew, Hnew, dt);
440  }
441  if (unstablePhase) {
442  ErrString += fmt::format(
443  "\t - The phase became unstable (Cp < 0) T_unstable_last = {}\n",
444  Tunstable);
445  }
446  if (doUV) {
447  throw CanteraError("ThermoPhase::setState_HPorUV (UV)", ErrString);
448  } else {
449  throw CanteraError("ThermoPhase::setState_HPorUV (HP)", ErrString);
450  }
451 }
452 
453 void ThermoPhase::setState_SP(double Starget, double p, double rtol)
454 {
455  setState_SPorSV(Starget, p, rtol, false);
456 }
457 
458 void ThermoPhase::setState_SV(double Starget, double v, double rtol)
459 {
460  assertCompressible("setState_SV");
461  setState_SPorSV(Starget, v, rtol, true);
462 }
463 
464 void ThermoPhase::setState_SPorSV(double Starget, double p,
465  double rtol, bool doSV)
466 {
467  double v = 0.0;
468  double dt;
469  if (doSV) {
470  v = p;
471  if (v < 1.0E-300) {
472  throw CanteraError("ThermoPhase::setState_SPorSV (SV)",
473  "Input specific volume is too small or negative. v = {}", v);
474  }
475  setDensity(1.0/v);
476  } else {
477  if (p < 1.0E-300) {
478  throw CanteraError("ThermoPhase::setState_SPorSV (SP)",
479  "Input pressure is too small or negative. p = {}", p);
480  }
481  setPressure(p);
482  }
483  double Tmax = maxTemp() + 0.1;
484  double Tmin = minTemp() - 0.1;
485 
486  // Make sure we are within the temperature bounds at the start
487  // of the iteration
488  double Tnew = temperature();
489  double Tinit = Tnew;
490  if (Tnew > Tmax) {
491  Tnew = Tmax - 1.0;
492  } else if (Tnew < Tmin) {
493  Tnew = Tmin + 1.0;
494  }
495  if (Tnew != Tinit) {
496  setState_conditional_TP(Tnew, p, !doSV);
497  }
498 
499  double Snew = entropy_mass();
500  double Cpnew = (doSV) ? cv_mass() : cp_mass();
501  double Stop = Snew;
502  double Ttop = Tnew;
503  double Sbot = Snew;
504  double Tbot = Tnew;
505 
506  bool ignoreBounds = false;
507  // Unstable phases are those for which Cp < 0.0. These are possible for
508  // cases where we have passed the spinodal curve.
509  bool unstablePhase = false;
510  double Tunstable = -1.0;
511  bool unstablePhaseNew = false;
512 
513  // Newton iteration
514  for (int n = 0; n < 500; n++) {
515  double Told = Tnew;
516  double Sold = Snew;
517  double cpd = Cpnew;
518  if (cpd < 0.0) {
519  unstablePhase = true;
520  Tunstable = Tnew;
521  }
522  // limit step size to 100 K
523  dt = clip((Starget - Sold)*Told/cpd, -100.0, 100.0);
524  Tnew = Told + dt;
525 
526  // Limit the step size so that we are convergent
527  if ((dt > 0.0 && unstablePhase) || (dt <= 0.0 && !unstablePhase)) {
528  if (Sbot < Starget && Tnew < Tbot) {
529  dt = 0.75 * (Tbot - Told);
530  Tnew = Told + dt;
531  }
532  } else if (Stop > Starget && Tnew > Ttop) {
533  dt = 0.75 * (Ttop - Told);
534  Tnew = Told + dt;
535  }
536 
537  // Check Max and Min values
538  if (Tnew > Tmax && !ignoreBounds) {
539  setState_conditional_TP(Tmax, p, !doSV);
540  double Smax = entropy_mass();
541  if (Smax >= Starget) {
542  if (Stop < Starget) {
543  Ttop = Tmax;
544  Stop = Smax;
545  }
546  } else {
547  Tnew = Tmax + 1.0;
548  ignoreBounds = true;
549  }
550  } else if (Tnew < Tmin && !ignoreBounds) {
551  setState_conditional_TP(Tmin, p, !doSV);
552  double Smin = entropy_mass();
553  if (Smin <= Starget) {
554  if (Sbot > Starget) {
555  Tbot = Tmin;
556  Sbot = Smin;
557  }
558  } else {
559  Tnew = Tmin - 1.0;
560  ignoreBounds = true;
561  }
562  }
563 
564  // Try to keep phase within its region of stability
565  // -> Could do a lot better if I calculate the
566  // spinodal value of H.
567  for (int its = 0; its < 10; its++) {
568  Tnew = Told + dt;
569  setState_conditional_TP(Tnew, p, !doSV);
570  Cpnew = (doSV) ? cv_mass() : cp_mass();
571  Snew = entropy_mass();
572  if (Cpnew < 0.0) {
573  unstablePhaseNew = true;
574  Tunstable = Tnew;
575  } else {
576  unstablePhaseNew = false;
577  break;
578  }
579  if (unstablePhase == false && unstablePhaseNew == true) {
580  dt *= 0.25;
581  }
582  }
583 
584  if (Snew == Starget) {
585  return;
586  } else if (Snew > Starget && (Stop < Starget || Snew < Stop)) {
587  Stop = Snew;
588  Ttop = Tnew;
589  } else if (Snew < Starget && (Sbot > Starget || Snew > Sbot)) {
590  Sbot = Snew;
591  Tbot = Tnew;
592  }
593  // Convergence in S
594  double Serr = Starget - Snew;
595  double acpd = std::max(fabs(cpd), 1.0E-5);
596  double denom = std::max(fabs(Starget), acpd * Tnew);
597  double SConvErr = fabs((Serr * Tnew)/denom);
598  if (SConvErr < rtol || fabs(dt/Tnew) < rtol) {
599  return;
600  }
601  }
602  // We are here when there hasn't been convergence
603 
604  // Formulate a detailed error message, since questions seem to arise often
605  // about the lack of convergence.
606  string ErrString = "No convergence in 500 iterations\n";
607  if (doSV) {
608  ErrString += fmt::format(
609  "\tTarget Entropy = {}\n"
610  "\tCurrent Specific Volume = {}\n"
611  "\tStarting Temperature = {}\n"
612  "\tCurrent Temperature = {}\n"
613  "\tCurrent Entropy = {}\n"
614  "\tCurrent Delta T = {}\n",
615  Starget, v, Tinit, Tnew, Snew, dt);
616  } else {
617  ErrString += fmt::format(
618  "\tTarget Entropy = {}\n"
619  "\tCurrent Pressure = {}\n"
620  "\tStarting Temperature = {}\n"
621  "\tCurrent Temperature = {}\n"
622  "\tCurrent Entropy = {}\n"
623  "\tCurrent Delta T = {}\n",
624  Starget, p, Tinit, Tnew, Snew, dt);
625  }
626  if (unstablePhase) {
627  ErrString += fmt::format("\t - The phase became unstable (Cp < 0) T_unstable_last = {}\n",
628  Tunstable);
629  }
630  if (doSV) {
631  throw CanteraError("ThermoPhase::setState_SPorSV (SV)", ErrString);
632  } else {
633  throw CanteraError("ThermoPhase::setState_SPorSV (SP)", ErrString);
634  }
635 }
636 
637 double ThermoPhase::o2Required(const double* y) const
638 {
639  // indices of fuel elements
640  size_t iC = elementIndex("C");
641  size_t iS = elementIndex("S");
642  size_t iH = elementIndex("H");
643 
644  double o2req = 0.0;
645  double sum = 0.0;
646  for (size_t k = 0; k != m_kk; ++k) {
647  sum += y[k];
648  double x = y[k] / molecularWeights()[k];
649  if (iC != npos) {
650  o2req += x * nAtoms(k, iC);
651  }
652  if (iS != npos) {
653  o2req += x * nAtoms(k, iS);
654  }
655  if (iH != npos) {
656  o2req += x * 0.25 * nAtoms(k, iH);
657  }
658  }
659  if (sum == 0.0) {
660  throw CanteraError("ThermoPhase::o2Required",
661  "No composition specified");
662  }
663  return o2req/sum;
664 }
665 
666 double ThermoPhase::o2Present(const double* y) const
667 {
668  size_t iO = elementIndex("O");
669  double o2pres = 0.0;
670  double sum = 0.0;
671  for (size_t k = 0; k != m_kk; ++k) {
672  sum += y[k];
673  o2pres += y[k] / molecularWeights()[k] * nAtoms(k, iO);
674  }
675  if (sum == 0.0) {
676  throw CanteraError("ThermoPhase::o2Present",
677  "No composition specified");
678  }
679  return 0.5 * o2pres / sum;
680 }
681 
682 double ThermoPhase::stoichAirFuelRatio(const Composition& fuelComp,
683  const Composition& oxComp,
684  ThermoBasis basis) const
685 {
686  vector<double> fuel(getCompositionFromMap(fuelComp));
687  vector<double> ox(getCompositionFromMap(oxComp));
688  return stoichAirFuelRatio(fuel.data(), ox.data(), basis);
689 }
690 
691 double ThermoPhase::stoichAirFuelRatio(const string& fuelComp, const string& oxComp,
692  ThermoBasis basis) const
693 {
694  return stoichAirFuelRatio(
695  parseCompString(fuelComp.find(":") != string::npos ? fuelComp : fuelComp+":1.0"),
696  parseCompString(oxComp.find(":") != string::npos ? oxComp : oxComp+":1.0"),
697  basis);
698 }
699 
700 double ThermoPhase::stoichAirFuelRatio(const double* fuelComp,
701  const double* oxComp,
702  ThermoBasis basis) const
703 {
704  vector<double> fuel, ox;
705  if (basis == ThermoBasis::molar) { // convert input compositions to mass fractions
706  fuel.resize(m_kk);
707  ox.resize(m_kk);
708  moleFractionsToMassFractions(fuelComp, fuel.data());
709  moleFractionsToMassFractions(oxComp, ox.data());
710  fuelComp = fuel.data();
711  oxComp = ox.data();
712  }
713 
714  double o2_required_fuel = o2Required(fuelComp) - o2Present(fuelComp);
715  double o2_required_ox = o2Required(oxComp) - o2Present(oxComp);
716 
717  if (o2_required_fuel < 0.0 || o2_required_ox > 0.0) {
718  throw CanteraError("ThermoPhase::stoichAirFuelRatio",
719  "Fuel composition contains too much oxygen or "
720  "oxidizer contains not enough oxygen. "
721  "Fuel and oxidizer composition mixed up?");
722  }
723 
724  if (o2_required_ox == 0.0) {
725  return std::numeric_limits<double>::infinity();
726  }
727 
728  return o2_required_fuel / (-o2_required_ox);
729 }
730 
731 void ThermoPhase::setEquivalenceRatio(double phi, const double* fuelComp,
732  const double* oxComp, ThermoBasis basis)
733 {
734  if (phi < 0.0) {
735  throw CanteraError("ThermoPhase::setEquivalenceRatio",
736  "Equivalence ratio phi must be >= 0");
737  }
738 
739  double p = pressure();
740 
741  vector<double> fuel, ox;
742  if (basis == ThermoBasis::molar) { // convert input compositions to mass fractions
743  fuel.resize(m_kk);
744  ox.resize(m_kk);
745  moleFractionsToMassFractions(fuelComp, fuel.data());
746  moleFractionsToMassFractions(oxComp, ox.data());
747  fuelComp = fuel.data();
748  oxComp = ox.data();
749  }
750 
751  double AFR_st = stoichAirFuelRatio(fuelComp, oxComp, ThermoBasis::mass);
752 
753  double sum_f = std::accumulate(fuelComp, fuelComp+m_kk, 0.0);
754  double sum_o = std::accumulate(oxComp, oxComp+m_kk, 0.0);
755 
756  vector<double> y(m_kk);
757  for (size_t k = 0; k != m_kk; ++k) {
758  y[k] = phi * fuelComp[k]/sum_f + AFR_st * oxComp[k]/sum_o;
759  }
760 
761  setMassFractions(y.data());
762  setPressure(p);
763 }
764 
765 void ThermoPhase::setEquivalenceRatio(double phi, const string& fuelComp,
766  const string& oxComp, ThermoBasis basis)
767 {
768  setEquivalenceRatio(phi,
769  parseCompString(fuelComp.find(":") != string::npos ? fuelComp : fuelComp+":1.0"),
770  parseCompString(oxComp.find(":") != string::npos ? oxComp : oxComp+":1.0"),
771  basis);
772 }
773 
774 void ThermoPhase::setEquivalenceRatio(double phi, const Composition& fuelComp,
775  const Composition& oxComp, ThermoBasis basis)
776 {
777  vector<double> fuel = getCompositionFromMap(fuelComp);
778  vector<double> ox = getCompositionFromMap(oxComp);
779  setEquivalenceRatio(phi, fuel.data(), ox.data(), basis);
780 }
781 
782 double ThermoPhase::equivalenceRatio() const
783 {
784  double o2_required = o2Required(massFractions());
785  double o2_present = o2Present(massFractions());
786 
787  if (o2_present == 0.0) { // pure fuel
788  return std::numeric_limits<double>::infinity();
789  }
790 
791  return o2_required / o2_present;
792 }
793 
794 double ThermoPhase::equivalenceRatio(const Composition& fuelComp,
795  const Composition& oxComp,
796  ThermoBasis basis) const
797 {
798  vector<double> fuel(getCompositionFromMap(fuelComp));
799  vector<double> ox(getCompositionFromMap(oxComp));
800  return equivalenceRatio(fuel.data(), ox.data(), basis);
801 }
802 
803 double ThermoPhase::equivalenceRatio(const string& fuelComp, const string& oxComp,
804  ThermoBasis basis) const
805 {
806  return equivalenceRatio(
807  parseCompString(fuelComp.find(":") != string::npos ? fuelComp : fuelComp+":1.0"),
808  parseCompString(oxComp.find(":") != string::npos ? oxComp : oxComp+":1.0"),
809  basis);
810 }
811 
812 double ThermoPhase::equivalenceRatio(const double* fuelComp,
813  const double* oxComp,
814  ThermoBasis basis) const
815 {
816  double Z = mixtureFraction(fuelComp, oxComp, basis);
817 
818  if (Z == 0.0) {
819  return 0.0; // pure oxidizer
820  }
821 
822  if (Z == 1.0) {
823  return std::numeric_limits<double>::infinity(); // pure fuel
824  }
825 
826  vector<double> fuel, ox;
827  if (basis == ThermoBasis::molar) { // convert input compositions to mass fractions
828  fuel.resize(m_kk);
829  ox.resize(m_kk);
830  moleFractionsToMassFractions(fuelComp, fuel.data());
831  moleFractionsToMassFractions(oxComp, ox.data());
832  fuelComp = fuel.data();
833  oxComp = ox.data();
834  }
835 
836  double AFR_st = stoichAirFuelRatio(fuelComp, oxComp, ThermoBasis::mass);
837 
838  return std::max(Z / (1.0 - Z) * AFR_st, 0.0);
839 }
840 
841 void ThermoPhase::setMixtureFraction(double mixFrac, const Composition& fuelComp,
842  const Composition& oxComp, ThermoBasis basis)
843 {
844  vector<double> fuel(getCompositionFromMap(fuelComp));
845  vector<double> ox(getCompositionFromMap(oxComp));
846  setMixtureFraction(mixFrac, fuel.data(), ox.data(), basis);
847 }
848 
849 void ThermoPhase::setMixtureFraction(double mixFrac, const string& fuelComp,
850  const string& oxComp, ThermoBasis basis)
851 {
852  setMixtureFraction(mixFrac,
853  parseCompString(fuelComp.find(":") != string::npos ? fuelComp : fuelComp+":1.0"),
854  parseCompString(oxComp.find(":") != string::npos ? oxComp : oxComp+":1.0"),
855  basis);
856 }
857 
858 void ThermoPhase::setMixtureFraction(double mixFrac, const double* fuelComp,
859  const double* oxComp, ThermoBasis basis)
860 {
861  if (mixFrac < 0.0 || mixFrac > 1.0) {
862  throw CanteraError("ThermoPhase::setMixtureFraction",
863  "Mixture fraction must be between 0 and 1");
864  }
865 
866  vector<double> fuel, ox;
867  if (basis == ThermoBasis::molar) { // convert input compositions to mass fractions
868  fuel.resize(m_kk);
869  ox.resize(m_kk);
870  moleFractionsToMassFractions(fuelComp, fuel.data());
871  moleFractionsToMassFractions(oxComp, ox.data());
872  fuelComp = fuel.data();
873  oxComp = ox.data();
874  }
875 
876  double sum_yf = std::accumulate(fuelComp, fuelComp+m_kk, 0.0);
877  double sum_yo = std::accumulate(oxComp, oxComp+m_kk, 0.0);
878 
879  if (sum_yf == 0.0 || sum_yo == 0.0) {
880  throw CanteraError("ThermoPhase::setMixtureFraction",
881  "No fuel and/or oxidizer composition specified");
882  }
883 
884  double p = pressure();
885 
886  vector<double> y(m_kk);
887 
888  for (size_t k = 0; k != m_kk; ++k) {
889  y[k] = mixFrac * fuelComp[k]/sum_yf + (1.0-mixFrac) * oxComp[k]/sum_yo;
890  }
891 
892  setMassFractions_NoNorm(y.data());
893  setPressure(p);
894 }
895 
896 double ThermoPhase::mixtureFraction(const Composition& fuelComp,
897  const Composition& oxComp,
898  ThermoBasis basis,
899  const string& element) const
900 {
901  vector<double> fuel(getCompositionFromMap(fuelComp));
902  vector<double> ox(getCompositionFromMap(oxComp));
903  return mixtureFraction(fuel.data(), ox.data(), basis, element);
904 }
905 
906 double ThermoPhase::mixtureFraction(const string& fuelComp, const string& oxComp,
907  ThermoBasis basis, const string& element) const
908 {
909  return mixtureFraction(
910  parseCompString(fuelComp.find(":") != string::npos ? fuelComp : fuelComp+":1.0"),
911  parseCompString(oxComp.find(":") != string::npos ? oxComp : oxComp+":1.0"),
912  basis, element);
913 }
914 
915 double ThermoPhase::mixtureFraction(const double* fuelComp, const double* oxComp,
916  ThermoBasis basis, const string& element) const
917 {
918  vector<double> fuel, ox;
919  if (basis == ThermoBasis::molar) { // convert input compositions to mass fractions
920  fuel.resize(m_kk);
921  ox.resize(m_kk);
922  moleFractionsToMassFractions(fuelComp, fuel.data());
923  moleFractionsToMassFractions(oxComp, ox.data());
924  fuelComp = fuel.data();
925  oxComp = ox.data();
926  }
927 
928  if (element == "Bilger") // compute the mixture fraction based on the Bilger mixture fraction
929  {
930  double o2_required_fuel = o2Required(fuelComp) - o2Present(fuelComp);
931  double o2_required_ox = o2Required(oxComp) - o2Present(oxComp);
932  double o2_required_mix = o2Required(massFractions()) - o2Present(massFractions());
933 
934  if (o2_required_fuel < 0.0 || o2_required_ox > 0.0) {
935  throw CanteraError("ThermoPhase::mixtureFraction",
936  "Fuel composition contains too much oxygen or "
937  "oxidizer contains not enough oxygen. "
938  "Fuel and oxidizer composition mixed up?");
939  }
940 
941  double denominator = o2_required_fuel - o2_required_ox;
942 
943  if (denominator == 0.0) {
944  throw CanteraError("ThermoPhase::mixtureFraction",
945  "Fuel and oxidizer have the same composition");
946  }
947 
948  double Z = (o2_required_mix - o2_required_ox) / denominator;
949 
950  return std::min(std::max(Z, 0.0), 1.0);
951  } else {
952  // compute the mixture fraction from a single element
953  double sum_yf = std::accumulate(fuelComp, fuelComp+m_kk, 0.0);
954  double sum_yo = std::accumulate(oxComp, oxComp+m_kk, 0.0);
955 
956  if (sum_yf == 0.0 || sum_yo == 0.0) {
957  throw CanteraError("ThermoPhase::mixtureFraction",
958  "No fuel and/or oxidizer composition specified");
959  }
960 
961  auto elementalFraction = [this](size_t m, const double* y) {
962  double Z_m = 0.0;
963  for (size_t k = 0; k != m_kk; ++k) {
964  Z_m += y[k] / molecularWeight(k) * nAtoms(k, m);
965  }
966  return Z_m;
967  };
968 
969  size_t m = elementIndex(element);
970  double Z_m_fuel = elementalFraction(m, fuelComp)/sum_yf;
971  double Z_m_ox = elementalFraction(m, oxComp)/sum_yo;
972  double Z_m_mix = elementalFraction(m, massFractions());
973 
974  if (Z_m_fuel == Z_m_ox) {
975  throw CanteraError("ThermoPhase::mixtureFraction",
976  "Fuel and oxidizer have the same composition for element {}",
977  element);
978  }
979  double Z = (Z_m_mix - Z_m_ox) / (Z_m_fuel - Z_m_ox);
980  return std::min(std::max(Z, 0.0), 1.0);
981  }
982 }
983 
984 MultiSpeciesThermo& ThermoPhase::speciesThermo(int k)
985 {
986  return m_spthermo;
987 }
988 
989 const MultiSpeciesThermo& ThermoPhase::speciesThermo(int k) const
990 {
991  return m_spthermo;
992 }
993 
994 
995 void ThermoPhase::initThermoFile(const string& inputFile, const string& id)
996 {
997  if (inputFile.empty()) {
998  // No input file specified - nothing to set up
999  return;
1000  }
1001  size_t dot = inputFile.find_last_of(".");
1002  string extension;
1003  if (dot != npos) {
1004  extension = inputFile.substr(dot+1);
1005  }
1006  if (extension == "xml" || extension == "cti") {
1007  throw CanteraError("ThermoPhase::initThermoFile",
1008  "The CTI and XML formats are no longer supported.");
1009  }
1010 
1011  AnyMap root = AnyMap::fromYamlFile(inputFile);
1012  auto& phase = root["phases"].getMapWhere("name", id);
1013  setupPhase(*this, phase, root);
1014 }
1015 
1016 void ThermoPhase::initThermo()
1017 {
1018  // Check to see that all of the species thermo objects have been initialized
1019  if (!m_spthermo.ready(m_kk)) {
1020  throw CanteraError("ThermoPhase::initThermo()",
1021  "Missing species thermo data");
1022  }
1023 }
1024 
1025 void ThermoPhase::setState_TPQ(double T, double P, double Q)
1026 {
1027  if (T > critTemperature()) {
1028  if (P > critPressure() || Q == 1) {
1029  setState_TP(T, P);
1030  return;
1031  } else {
1032  throw CanteraError("ThermoPhase::setState_TPQ",
1033  "Temperature ({}), pressure ({}) and vapor fraction ({}) "
1034  "are inconsistent, above the critical temperature.",
1035  T, P, Q);
1036  }
1037  }
1038 
1039  double Psat = satPressure(T);
1040  if (std::abs(Psat / P - 1) < 1e-6) {
1041  setState_Tsat(T, Q);
1042  } else if ((Q == 0 && P >= Psat) || (Q == 1 && P <= Psat)) {
1043  setState_TP(T, P);
1044  } else {
1045  throw CanteraError("ThermoPhase::setState_TPQ",
1046  "Temperature ({}), pressure ({}) and vapor fraction ({}) "
1047  "are inconsistent.\nPsat at this T: {}\n"
1048  "Consider specifying the state using two fully independent "
1049  "properties (for example, temperature and density)",
1050  T, P, Q, Psat);
1051  }
1052 }
1053 
1054 bool ThermoPhase::addSpecies(shared_ptr<Species> spec)
1055 {
1056  if (!spec->thermo) {
1057  throw CanteraError("ThermoPhase::addSpecies",
1058  "Species {} has no thermo data", spec->name);
1059  }
1060  bool added = Phase::addSpecies(spec);
1061  if (added) {
1062  spec->thermo->validate(spec->name);
1063  m_spthermo.install_STIT(m_kk-1, spec->thermo);
1064  }
1065  return added;
1066 }
1067 
1068 void ThermoPhase::modifySpecies(size_t k, shared_ptr<Species> spec)
1069 {
1070  if (!spec->thermo) {
1071  throw CanteraError("ThermoPhase::modifySpecies",
1072  "Species {} has no thermo data", spec->name);
1073  }
1074  Phase::modifySpecies(k, spec);
1075  if (speciesName(k) != spec->name) {
1076  throw CanteraError("ThermoPhase::modifySpecies",
1077  "New species '{}' does not match existing species '{}' at index {}",
1078  spec->name, speciesName(k), k);
1079  }
1080  spec->thermo->validate(spec->name);
1081  m_spthermo.modifySpecies(k, spec->thermo);
1082 }
1083 
1084 void ThermoPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode)
1085 {
1086  m_input = phaseNode;
1087 }
1088 
1089 AnyMap ThermoPhase::parameters(bool withInput) const
1090 {
1091  AnyMap out;
1092  getParameters(out);
1093  if (withInput) {
1094  out.update(m_input);
1095  }
1096  return out;
1097 }
1098 
1099 void ThermoPhase::getParameters(AnyMap& phaseNode) const
1100 {
1101  phaseNode["name"] = name();
1102  phaseNode["thermo"] = ThermoFactory::factory()->canonicalize(type());
1103  vector<string> elementNames;
1104  for (size_t i = 0; i < nElements(); i++) {
1105  elementNames.push_back(elementName(i));
1106  }
1107  phaseNode["elements"] = elementNames;
1108  phaseNode["species"] = speciesNames();
1109 
1110  AnyMap state;
1111  auto stateVars = nativeState();
1112  if (stateVars.count("T")) {
1113  state["T"].setQuantity(temperature(), "K");
1114  }
1115 
1116  if (stateVars.count("D")) {
1117  state["density"].setQuantity(density(), "kg/m^3");
1118  } else if (stateVars.count("P")) {
1119  state["P"].setQuantity(pressure(), "Pa");
1120  }
1121 
1122  if (stateVars.count("Y")) {
1123  map<string, double> Y;
1124  for (size_t k = 0; k < m_kk; k++) {
1125  double Yk = massFraction(k);
1126  if (Yk > 0) {
1127  Y[speciesName(k)] = Yk;
1128  }
1129  }
1130  state["Y"] = Y;
1131  state["Y"].setFlowStyle();
1132  } else if (stateVars.count("X")) {
1133  map<string, double> X;
1134  for (size_t k = 0; k < m_kk; k++) {
1135  double Xk = moleFraction(k);
1136  if (Xk > 0) {
1137  X[speciesName(k)] = Xk;
1138  }
1139  }
1140  state["X"] = X;
1141  state["X"].setFlowStyle();
1142  }
1143 
1144  phaseNode["state"] = std::move(state);
1145 
1146  static bool reg = AnyMap::addOrderingRules("Phase", {{"tail", "state"}});
1147  if (reg) {
1148  phaseNode["__type__"] = "Phase";
1149  }
1150 }
1151 
1152 const AnyMap& ThermoPhase::input() const
1153 {
1154  return m_input;
1155 }
1156 
1157 AnyMap& ThermoPhase::input()
1158 {
1159  return m_input;
1160 }
1161 
1162 void ThermoPhase::invalidateCache() {
1163  Phase::invalidateCache();
1164  m_tlast += 0.1234;
1165 }
1166 
1167 void ThermoPhase::equilibrate(const string& XY, const string& solver,
1168  double rtol, int max_steps, int max_iter,
1169  int estimate_equil, int log_level)
1170 {
1171  if (solver == "auto" || solver == "element_potential") {
1172  vector<double> initial_state;
1173  saveState(initial_state);
1174  debuglog("Trying ChemEquil solver\n", log_level);
1175  try {
1176  ChemEquil E;
1177  E.options.maxIterations = max_steps;
1178  E.options.relTolerance = rtol;
1179  int ret = E.equilibrate(*this, XY.c_str(), log_level-1);
1180  if (ret < 0) {
1181  throw CanteraError("ThermoPhase::equilibrate",
1182  "ChemEquil solver failed. Return code: {}", ret);
1183  }
1184  debuglog("ChemEquil solver succeeded\n", log_level);
1185  return;
1186  } catch (std::exception& err) {
1187  debuglog("ChemEquil solver failed.\n", log_level);
1188  debuglog(err.what(), log_level);
1189  restoreState(initial_state);
1190  if (solver == "auto") {
1191  } else {
1192  throw;
1193  }
1194  }
1195  }
1196 
1197  if (solver == "auto" || solver == "vcs" || solver == "gibbs") {
1198  MultiPhase M;
1199  M.addPhase(this, 1.0);
1200  M.init();
1201  M.equilibrate(XY, solver, rtol, max_steps, max_iter,
1202  estimate_equil, log_level);
1203  return;
1204  }
1205 
1206  if (solver != "auto") {
1207  throw CanteraError("ThermoPhase::equilibrate",
1208  "Invalid solver specified: '{}'", solver);
1209  }
1210 }
1211 
1212 void ThermoPhase::getdlnActCoeffdlnN(const size_t ld, double* const dlnActCoeffdlnN)
1213 {
1214  for (size_t m = 0; m < m_kk; m++) {
1215  for (size_t k = 0; k < m_kk; k++) {
1216  dlnActCoeffdlnN[ld * k + m] = 0.0;
1217  }
1218  }
1219  return;
1220 }
1221 
1222 void ThermoPhase::getdlnActCoeffdlnN_numderiv(const size_t ld,
1223  double* const dlnActCoeffdlnN)
1224 {
1225  double deltaMoles_j = 0.0;
1226  double pres = pressure();
1227 
1228  // Evaluate the current base activity coefficients if necessary
1229  vector<double> ActCoeff_Base(m_kk);
1230  getActivityCoefficients(ActCoeff_Base.data());
1231  vector<double> Xmol_Base(m_kk);
1232  getMoleFractions(Xmol_Base.data());
1233 
1234  // Make copies of ActCoeff and Xmol_ for use in taking differences
1235  vector<double> ActCoeff(m_kk);
1236  vector<double> Xmol(m_kk);
1237  double v_totalMoles = 1.0;
1238  double TMoles_base = v_totalMoles;
1239 
1240  // Loop over the columns species to be deltad
1241  for (size_t j = 0; j < m_kk; j++) {
1242  // Calculate a value for the delta moles of species j
1243  // -> Note Xmol_[] and Tmoles are always positive or zero quantities.
1244  // -> experience has shown that you always need to make the deltas
1245  // greater than needed to change the other mole fractions in order
1246  // to capture some effects.
1247  double moles_j_base = v_totalMoles * Xmol_Base[j];
1248  deltaMoles_j = 1.0E-7 * moles_j_base + v_totalMoles * 1.0E-13 + 1.0E-150;
1249 
1250  // Now, update the total moles in the phase and all of the mole
1251  // fractions based on this.
1252  v_totalMoles = TMoles_base + deltaMoles_j;
1253  for (size_t k = 0; k < m_kk; k++) {
1254  Xmol[k] = Xmol_Base[k] * TMoles_base / v_totalMoles;
1255  }
1256  Xmol[j] = (moles_j_base + deltaMoles_j) / v_totalMoles;
1257 
1258  // Go get new values for the activity coefficients.
1259  setMoleFractions(Xmol.data());
1260  setPressure(pres);
1261  getActivityCoefficients(ActCoeff.data());
1262 
1263  // Calculate the column of the matrix
1264  double* const lnActCoeffCol = dlnActCoeffdlnN + ld * j;
1265  for (size_t k = 0; k < m_kk; k++) {
1266  lnActCoeffCol[k] = (2*moles_j_base + deltaMoles_j) *(ActCoeff[k] - ActCoeff_Base[k]) /
1267  ((ActCoeff[k] + ActCoeff_Base[k]) * deltaMoles_j);
1268  }
1269  // Revert to the base case Xmol_, v_totalMoles
1270  v_totalMoles = TMoles_base;
1271  Xmol = Xmol_Base;
1272  }
1273  setMoleFractions(Xmol_Base.data());
1274  setPressure(pres);
1275 }
1276 
1277 string ThermoPhase::report(bool show_thermo, double threshold) const
1278 {
1279  if (type() == "none") {
1280  throw NotImplementedError("ThermoPhase::report",
1281  "Not implemented for thermo model 'none'");
1282  }
1283 
1284  fmt::memory_buffer b;
1285  // This is the width of the first column of names in the report.
1286  int name_width = 18;
1287 
1288  string blank_leader = fmt::format("{:{}}", "", name_width);
1289 
1290  string string_property = fmt::format("{{:>{}}} {{}}\n", name_width);
1291 
1292  string one_property = fmt::format("{{:>{}}} {{:<.5g}} {{}}\n", name_width);
1293 
1294  string two_prop_header = "{} {:^15} {:^15}\n";
1295  string kg_kmol_header = fmt::format(
1296  two_prop_header, blank_leader, "1 kg", "1 kmol"
1297  );
1298  string Y_X_header = fmt::format(
1299  two_prop_header, blank_leader, "mass frac. Y", "mole frac. X"
1300  );
1301  string two_prop_sep = fmt::format(
1302  "{} {:-^15} {:-^15}\n", blank_leader, "", ""
1303  );
1304  string two_property = fmt::format(
1305  "{{:>{}}} {{:15.5g}} {{:15.5g}} {{}}\n", name_width
1306  );
1307 
1308  string three_prop_header = fmt::format(
1309  "{} {:^15} {:^15} {:^15}\n", blank_leader, "mass frac. Y",
1310  "mole frac. X", "chem. pot. / RT"
1311  );
1312  string three_prop_sep = fmt::format(
1313  "{} {:-^15} {:-^15} {:-^15}\n", blank_leader, "", "", ""
1314  );
1315  string three_property = fmt::format(
1316  "{{:>{}}} {{:15.5g}} {{:15.5g}} {{:15.5g}}\n", name_width
1317  );
1318 
1319  try {
1320  if (name() != "") {
1321  fmt_append(b, "\n {}:\n", name());
1322  }
1323  fmt_append(b, "\n");
1324  fmt_append(b, one_property, "temperature", temperature(), "K");
1325  fmt_append(b, one_property, "pressure", pressure(), "Pa");
1326  fmt_append(b, one_property, "density", density(), "kg/m^3");
1327  fmt_append(b, one_property,
1328  "mean mol. weight", meanMolecularWeight(), "kg/kmol");
1329 
1330  double phi = electricPotential();
1331  if (phi != 0.0) {
1332  fmt_append(b, one_property, "potential", phi, "V");
1333  }
1334 
1335  fmt_append(b, string_property, "phase of matter", phaseOfMatter());
1336 
1337  if (show_thermo) {
1338  fmt_append(b, "\n");
1339  fmt_append(b, kg_kmol_header);
1340  fmt_append(b, two_prop_sep);
1341  fmt_append(b, two_property,
1342  "enthalpy", enthalpy_mass(), enthalpy_mole(), "J");
1343  fmt_append(b, two_property,
1344  "internal energy", intEnergy_mass(), intEnergy_mole(), "J");
1345  fmt_append(b, two_property,
1346  "entropy", entropy_mass(), entropy_mole(), "J/K");
1347  fmt_append(b, two_property,
1348  "Gibbs function", gibbs_mass(), gibbs_mole(), "J");
1349  fmt_append(b, two_property,
1350  "heat capacity c_p", cp_mass(), cp_mole(), "J/K");
1351  try {
1352  fmt_append(b, two_property,
1353  "heat capacity c_v", cv_mass(), cv_mole(), "J/K");
1354  } catch (NotImplementedError&) {
1355  fmt_append(b, string_property,
1356  "heat capacity c_v", "<not implemented>");
1357  }
1358  }
1359 
1360  vector<double> x(m_kk);
1361  vector<double> y(m_kk);
1362  vector<double> mu(m_kk);
1363  getMoleFractions(&x[0]);
1364  getMassFractions(&y[0]);
1365  getChemPotentials(&mu[0]);
1366  int nMinor = 0;
1367  double xMinor = 0.0;
1368  double yMinor = 0.0;
1369  fmt_append(b, "\n");
1370  if (show_thermo) {
1371  fmt_append(b, three_prop_header);
1372  fmt_append(b, three_prop_sep);
1373  for (size_t k = 0; k < m_kk; k++) {
1374  if (abs(x[k]) >= threshold) {
1375  if (abs(x[k]) > SmallNumber) {
1376  fmt_append(b, three_property,
1377  speciesName(k), y[k], x[k], mu[k]/RT());
1378  } else {
1379  fmt_append(b, two_property, speciesName(k), y[k], x[k], "");
1380  }
1381  } else {
1382  nMinor++;
1383  xMinor += x[k];
1384  yMinor += y[k];
1385  }
1386  }
1387  } else {
1388  fmt_append(b, Y_X_header);
1389  fmt_append(b, two_prop_sep);
1390  for (size_t k = 0; k < m_kk; k++) {
1391  if (abs(x[k]) >= threshold) {
1392  fmt_append(b, two_property, speciesName(k), y[k], x[k], "");
1393  } else {
1394  nMinor++;
1395  xMinor += x[k];
1396  yMinor += y[k];
1397  }
1398  }
1399  }
1400  if (nMinor) {
1401  string minor = fmt::format("[{:+5d} minor]", nMinor);
1402  fmt_append(b, two_property, minor, yMinor, xMinor, "");
1403  }
1404  } catch (CanteraError& err) {
1405  return to_string(b) + err.what();
1406  }
1407  return to_string(b);
1408 }
1409 
1410 }
Chemical equilibrium.
Headers for the MultiPhase object that is used to set up multiphase equilibrium problems (see Chemica...
Pure Virtual Base class for individual species reference state thermodynamic managers and text for th...
Declaration for class Cantera::Species.
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
A map of string keys to values whose type can vary at runtime.
Definition: AnyMap.h:427
size_t size() const
Returns the number of elements in this map.
Definition: AnyMap.h:622
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition: AnyMap.cpp:1423
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
void setFlowStyle(bool flow=true)
Use "flow" style when outputting this AnyMap to YAML.
Definition: AnyMap.cpp:1726
void erase(const string &key)
Erase the value held by key.
Definition: AnyMap.cpp:1428
void update(const AnyMap &other, bool keepExisting=true)
Add items from other to this AnyMap.
Definition: AnyMap.cpp:1438
string keys_str() const
Return a string listing the keys in this AnyMap, for use in error messages, for example.
Definition: AnyMap.cpp:1447
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
const char * what() const override
Get a description of the error.
Class ChemEquil implements a chemical equilibrium solver for single-phase solutions.
Definition: ChemEquil.h:79
int equilibrate(ThermoPhase &s, const char *XY, int loglevel=0)
Equilibrate a phase, holding the elemental composition fixed at the initial value found within the Th...
Definition: ChemEquil.cpp:293
EquilOpt options
Options controlling how the calculation is carried out.
Definition: ChemEquil.h:127
double relTolerance
Relative tolerance.
Definition: ChemEquil.h:29
int maxIterations
Maximum number of iterations.
Definition: ChemEquil.h:31
A class for multiphase mixtures.
Definition: MultiPhase.h:57
void init()
Process phases and build atomic composition array.
Definition: MultiPhase.cpp:106
void addPhase(ThermoPhase *p, double moles)
Add a phase to the mixture.
Definition: MultiPhase.cpp:39
A species thermodynamic property manager for a phase.
An error indicating that an unimplemented function has been called.
Definition: ctexceptions.h:195
A representation of the units associated with a dimensional quantity.
Definition: Units.h:35
void fmt_append(fmt::memory_buffer &b, Args... args)
Versions 6.2.0 and 6.2.1 of fmtlib do not include this define before they include windows....
Definition: fmt.h:29
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
void equilibrate(const string &XY, const string &solver="auto", double rtol=1e-9, int max_steps=50000, int max_iter=100, int estimate_equil=0, int log_level=0)
Equilibrate a MultiPhase object.
Definition: MultiPhase.cpp:630
void debuglog(const string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Definition: global.h:158
double dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
Definition: utilities.h:82
T clip(const T &value, const T &lower, const T &upper)
Clip value such that lower <= value <= upper.
Definition: global.h:329
const double Faraday
Faraday constant [C/kmol].
Definition: ct_defs.h:131
const double OneAtm
One atmosphere [Pa].
Definition: ct_defs.h:96
void setupPhase(ThermoPhase &thermo, const AnyMap &phaseNode, const AnyMap &rootNode)
Initialize a ThermoPhase object.
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
const vector< string > & elementNames()
Get a vector of the names of the elements defined in Cantera.
Definition: Elements.cpp:227
const int cAC_CONVENTION_MOLAR
Standard state uses the molar convention.
Definition: ThermoPhase.h:135
const double SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:158
ThermoBasis
Differentiate between mole fractions and mass fractions for input mixture composition.
Definition: ThermoPhase.h:154
map< string, double > Composition
Map from string names to doubles.
Definition: ct_defs.h:177
Contains declarations for string manipulation functions within Cantera.