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