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