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