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