28 void ThermoPhase::resetHf298(
size_t k) {
30 m_spthermo.resetHf298(k);
32 for (
size_t k = 0; k < nSpecies(); k++) {
33 m_spthermo.resetHf298(k);
39 int ThermoPhase::activityConvention()
const
44 int ThermoPhase::standardStateConvention()
const
46 return m_ssConvention;
49 Units ThermoPhase::standardConcentrationUnits()
const
52 return Units(1.0, 0, -
static_cast<double>(nDim()), 0, 0, 0, 1);
55 double ThermoPhase::logStandardConc(
size_t k)
const
57 return log(standardConcentration(k));
60 void ThermoPhase::getActivities(
double* a)
const
62 getActivityConcentrations(a);
63 for (
size_t k = 0; k < nSpecies(); k++) {
64 a[k] /= standardConcentration(k);
68 void ThermoPhase::getLnActivityCoefficients(
double* lnac)
const
70 getActivityCoefficients(lnac);
71 for (
size_t k = 0; k < m_kk; k++) {
72 lnac[k] = std::log(lnac[k]);
76 void ThermoPhase::getElectrochemPotentials(
double* mu)
const
78 getChemPotentials(mu);
79 double ve =
Faraday * electricPotential();
80 for (
size_t k = 0; k < m_kk; k++) {
81 mu[k] += ve*charge(k);
85 void ThermoPhase::setState_TPX(
double t,
double p,
const double* x)
91 void ThermoPhase::setState_TPX(
double t,
double p,
const Composition& x)
93 setMoleFractionsByName(x);
97 void ThermoPhase::setState_TPX(
double t,
double p,
const string& x)
99 setMoleFractionsByName(x);
103 void ThermoPhase::setState_TPY(
double t,
double p,
const double* y)
109 void ThermoPhase::setState_TPY(
double t,
double p,
const Composition& y)
111 setMassFractionsByName(y);
115 void ThermoPhase::setState_TPY(
double t,
double p,
const string& y)
117 setMassFractionsByName(y);
121 void ThermoPhase::setState_TP(
double t,
double p)
123 double tsave = temperature();
124 double dsave = density();
129 setState_TD(tsave, dsave);
134 void ThermoPhase::setState_HP(
double Htarget,
double p,
double rtol)
136 setState_HPorUV(Htarget, p, rtol,
false);
139 void ThermoPhase::setState_UV(
double u,
double v,
double rtol)
141 assertCompressible(
"setState_UV");
142 setState_HPorUV(u, v, rtol,
true);
145 void ThermoPhase::setState(
const AnyMap& input_state)
147 AnyMap state = input_state;
150 if (state.
hasKey(
"mass-fractions")) {
151 state[
"Y"] = state[
"mass-fractions"];
152 state.
erase(
"mass-fractions");
154 if (state.
hasKey(
"mole-fractions")) {
155 state[
"X"] = state[
"mole-fractions"];
156 state.
erase(
"mole-fractions");
158 if (state.
hasKey(
"temperature")) {
159 state[
"T"] = state[
"temperature"];
161 if (state.
hasKey(
"pressure")) {
162 state[
"P"] = state[
"pressure"];
164 if (state.
hasKey(
"enthalpy")) {
165 state[
"H"] = state[
"enthalpy"];
167 if (state.
hasKey(
"int-energy")) {
168 state[
"U"] = state[
"int-energy"];
170 if (state.
hasKey(
"internal-energy")) {
171 state[
"U"] = state[
"internal-energy"];
173 if (state.
hasKey(
"specific-volume")) {
174 state[
"V"] = state[
"specific-volume"];
176 if (state.
hasKey(
"entropy")) {
177 state[
"S"] = state[
"entropy"];
179 if (state.
hasKey(
"density")) {
180 state[
"D"] = state[
"density"];
182 if (state.
hasKey(
"vapor-fraction")) {
183 state[
"Q"] = state[
"vapor-fraction"];
188 if (state[
"X"].is<string>()) {
189 setMoleFractionsByName(state[
"X"].asString());
191 setMoleFractionsByName(state[
"X"].asMap<double>());
194 }
else if (state.
hasKey(
"Y")) {
195 if (state[
"Y"].is<string>()) {
196 setMassFractionsByName(state[
"Y"].asString());
198 setMassFractionsByName(state[
"Y"].asMap<double>());
203 if (state.
size() == 0) {
204 setState_TP(298.15,
OneAtm);
206 double T = state.
convert(
"T",
"K");
207 double P = state.
convert(
"P",
"Pa");
209 setState_TPQ(T, P, state[
"Q"].asDouble());
214 setState_TD(state.
convert(
"T",
"K"), state.
convert(
"D",
"kg/m^3"));
216 setState_TV(state.
convert(
"T",
"K"), state.
convert(
"V",
"m^3/kg"));
218 setState_HP(state.
convert(
"H",
"J/kg"), state.
convert(
"P",
"Pa"));
220 setState_UV(state.
convert(
"U",
"J/kg"), state.
convert(
"V",
"m^3/kg"));
222 setState_SP(state.
convert(
"S",
"J/kg/K"), state.
convert(
"P",
"Pa"));
224 setState_SV(state.
convert(
"S",
"J/kg/K"), state.
convert(
"V",
"m^3/kg"));
226 setState_ST(state.
convert(
"S",
"J/kg/K"), state.
convert(
"T",
"K"));
228 setState_PV(state.
convert(
"P",
"Pa"), state.
convert(
"V",
"m^3/kg"));
230 setState_UP(state.
convert(
"U",
"J/kg"), state.
convert(
"P",
"Pa"));
232 setState_VH(state.
convert(
"V",
"m^3/kg"), state.
convert(
"H",
"J/kg"));
236 setState_SH(state.
convert(
"S",
"J/kg/K"), state.
convert(
"H",
"J/kg"));
238 setState_DP(state.
convert(
"D",
"kg/m^3"), state.
convert(
"P",
"Pa"));
240 setState_Psat(state.
convert(
"P",
"Pa"), state[
"Q"].asDouble());
242 setState_Tsat(state.
convert(
"T",
"K"), state[
"Q"].asDouble());
243 }
else if (state.
hasKey(
"T")) {
245 }
else if (state.
hasKey(
"P")) {
246 setState_TP(298.15, state.
convert(
"P",
"Pa"));
249 "'state' did not specify a recognized set of properties.\n"
250 "Keys provided were: {}", input_state.
keys_str());
254 void ThermoPhase::setState_conditional_TP(
double t,
double p,
bool set_p)
262 void ThermoPhase::setState_HPorUV(
double Htarget,
double p,
263 double rtol,
bool doUV)
273 "Input specific volume is too small or negative. v = {}", v);
279 "Input pressure is too small or negative. p = {}", p);
283 double Tmax = maxTemp() + 0.1;
284 double Tmin = minTemp() - 0.1;
288 double Tnew = temperature();
292 }
else if (Tnew < Tmin) {
296 setState_conditional_TP(Tnew, p, !doUV);
299 double Hnew = (doUV) ? intEnergy_mass() : enthalpy_mass();
300 double Cpnew = (doUV) ? cv_mass() : cp_mass();
306 bool ignoreBounds =
false;
309 bool unstablePhase =
false;
312 double Tunstable = -1.0;
313 bool unstablePhaseNew =
false;
316 for (
int n = 0; n < 500; n++) {
321 unstablePhase =
true;
325 dt =
clip((Htarget - Hold)/cpd, -100.0, 100.0);
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);
337 }
else if (Htop > Htarget && Tnew > (0.75 * Ttop + 0.25 * Told)) {
338 dt = 0.75 * (Ttop - Told);
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) {
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) {
373 for (
int its = 0; its < 10; its++) {
375 if (Tnew < Told / 3.0) {
377 dt = -2.0 * Told / 3.0;
379 setState_conditional_TP(Tnew, p, !doUV);
381 Hnew = intEnergy_mass();
384 Hnew = enthalpy_mass();
388 unstablePhaseNew =
true;
391 unstablePhaseNew =
false;
394 if (unstablePhase ==
false && unstablePhaseNew ==
true) {
399 if (Hnew == Htarget) {
401 }
else if (Hnew > Htarget && (Htop < Htarget || Hnew < Htop)) {
404 }
else if (Hnew < Htarget && (Hbot > Htarget || Hnew > Hbot)) {
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) {
421 string ErrString =
"No convergence in 500 iterations\n";
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);
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);
442 ErrString += fmt::format(
443 "\t - The phase became unstable (Cp < 0) T_unstable_last = {}\n",
447 throw CanteraError(
"ThermoPhase::setState_HPorUV (UV)", ErrString);
449 throw CanteraError(
"ThermoPhase::setState_HPorUV (HP)", ErrString);
453 void ThermoPhase::setState_SP(
double Starget,
double p,
double rtol)
455 setState_SPorSV(Starget, p, rtol,
false);
458 void ThermoPhase::setState_SV(
double Starget,
double v,
double rtol)
460 assertCompressible(
"setState_SV");
461 setState_SPorSV(Starget, v, rtol,
true);
464 void ThermoPhase::setState_SPorSV(
double Starget,
double p,
465 double rtol,
bool doSV)
473 "Input specific volume is too small or negative. v = {}", v);
479 "Input pressure is too small or negative. p = {}", p);
483 double Tmax = maxTemp() + 0.1;
484 double Tmin = minTemp() - 0.1;
488 double Tnew = temperature();
492 }
else if (Tnew < Tmin) {
496 setState_conditional_TP(Tnew, p, !doSV);
499 double Snew = entropy_mass();
500 double Cpnew = (doSV) ? cv_mass() : cp_mass();
506 bool ignoreBounds =
false;
509 bool unstablePhase =
false;
510 double Tunstable = -1.0;
511 bool unstablePhaseNew =
false;
514 for (
int n = 0; n < 500; n++) {
519 unstablePhase =
true;
523 dt =
clip((Starget - Sold)*Told/cpd, -100.0, 100.0);
527 if ((dt > 0.0 && unstablePhase) || (dt <= 0.0 && !unstablePhase)) {
528 if (Sbot < Starget && Tnew < Tbot) {
529 dt = 0.75 * (Tbot - Told);
532 }
else if (Stop > Starget && Tnew > Ttop) {
533 dt = 0.75 * (Ttop - Told);
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) {
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) {
567 for (
int its = 0; its < 10; its++) {
569 setState_conditional_TP(Tnew, p, !doSV);
570 Cpnew = (doSV) ? cv_mass() : cp_mass();
571 Snew = entropy_mass();
573 unstablePhaseNew =
true;
576 unstablePhaseNew =
false;
579 if (unstablePhase ==
false && unstablePhaseNew ==
true) {
584 if (Snew == Starget) {
586 }
else if (Snew > Starget && (Stop < Starget || Snew < Stop)) {
589 }
else if (Snew < Starget && (Sbot > Starget || Snew > Sbot)) {
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) {
606 string ErrString =
"No convergence in 500 iterations\n";
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);
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);
627 ErrString += fmt::format(
"\t - The phase became unstable (Cp < 0) T_unstable_last = {}\n",
631 throw CanteraError(
"ThermoPhase::setState_SPorSV (SV)", ErrString);
633 throw CanteraError(
"ThermoPhase::setState_SPorSV (SP)", ErrString);
637 double ThermoPhase::o2Required(
const double* y)
const
640 size_t iC = elementIndex(
"C");
641 size_t iS = elementIndex(
"S");
642 size_t iH = elementIndex(
"H");
646 for (
size_t k = 0; k != m_kk; ++k) {
648 double x = y[k] / molecularWeights()[k];
650 o2req += x * nAtoms(k, iC);
653 o2req += x * nAtoms(k, iS);
656 o2req += x * 0.25 * nAtoms(k, iH);
661 "No composition specified");
666 double ThermoPhase::o2Present(
const double* y)
const
668 size_t iO = elementIndex(
"O");
671 for (
size_t k = 0; k != m_kk; ++k) {
673 o2pres += y[k] / molecularWeights()[k] * nAtoms(k, iO);
677 "No composition specified");
679 return 0.5 * o2pres / sum;
682 double ThermoPhase::stoichAirFuelRatio(
const Composition& fuelComp,
686 vector<double> fuel(getCompositionFromMap(fuelComp));
687 vector<double> ox(getCompositionFromMap(oxComp));
688 return stoichAirFuelRatio(fuel.data(), ox.data(), basis);
691 double ThermoPhase::stoichAirFuelRatio(
const string& fuelComp,
const string& oxComp,
694 return stoichAirFuelRatio(
700 double ThermoPhase::stoichAirFuelRatio(
const double* fuelComp,
701 const double* oxComp,
704 vector<double> fuel, ox;
705 if (basis == ThermoBasis::molar) {
708 moleFractionsToMassFractions(fuelComp, fuel.data());
709 moleFractionsToMassFractions(oxComp, ox.data());
710 fuelComp = fuel.data();
714 double o2_required_fuel = o2Required(fuelComp) - o2Present(fuelComp);
715 double o2_required_ox = o2Required(oxComp) - o2Present(oxComp);
717 if (o2_required_fuel < 0.0 || o2_required_ox > 0.0) {
719 "Fuel composition contains too much oxygen or "
720 "oxidizer contains not enough oxygen. "
721 "Fuel and oxidizer composition mixed up?");
724 if (o2_required_ox == 0.0) {
725 return std::numeric_limits<double>::infinity();
728 return o2_required_fuel / (-o2_required_ox);
731 void ThermoPhase::setEquivalenceRatio(
double phi,
const double* fuelComp,
736 "Equivalence ratio phi must be >= 0");
739 double p = pressure();
741 vector<double> fuel, ox;
742 if (basis == ThermoBasis::molar) {
745 moleFractionsToMassFractions(fuelComp, fuel.data());
746 moleFractionsToMassFractions(oxComp, ox.data());
747 fuelComp = fuel.data();
751 double AFR_st = stoichAirFuelRatio(fuelComp, oxComp, ThermoBasis::mass);
753 double sum_f = std::accumulate(fuelComp, fuelComp+m_kk, 0.0);
754 double sum_o = std::accumulate(oxComp, oxComp+m_kk, 0.0);
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;
761 setMassFractions(y.data());
765 void ThermoPhase::setEquivalenceRatio(
double phi,
const string& fuelComp,
768 setEquivalenceRatio(phi,
774 void ThermoPhase::setEquivalenceRatio(
double phi,
const Composition& fuelComp,
777 vector<double> fuel = getCompositionFromMap(fuelComp);
778 vector<double> ox = getCompositionFromMap(oxComp);
779 setEquivalenceRatio(phi, fuel.data(), ox.data(), basis);
782 double ThermoPhase::equivalenceRatio()
const
784 double o2_required = o2Required(massFractions());
785 double o2_present = o2Present(massFractions());
787 if (o2_present == 0.0) {
788 return std::numeric_limits<double>::infinity();
791 return o2_required / o2_present;
798 vector<double> fuel(getCompositionFromMap(fuelComp));
799 vector<double> ox(getCompositionFromMap(oxComp));
800 return equivalenceRatio(fuel.data(), ox.data(), basis);
803 double ThermoPhase::equivalenceRatio(
const string& fuelComp,
const string& oxComp,
806 return equivalenceRatio(
812 double ThermoPhase::equivalenceRatio(
const double* fuelComp,
813 const double* oxComp,
816 double Z = mixtureFraction(fuelComp, oxComp, basis);
823 return std::numeric_limits<double>::infinity();
826 vector<double> fuel, ox;
827 if (basis == ThermoBasis::molar) {
830 moleFractionsToMassFractions(fuelComp, fuel.data());
831 moleFractionsToMassFractions(oxComp, ox.data());
832 fuelComp = fuel.data();
836 double AFR_st = stoichAirFuelRatio(fuelComp, oxComp, ThermoBasis::mass);
838 return std::max(Z / (1.0 - Z) * AFR_st, 0.0);
841 void ThermoPhase::setMixtureFraction(
double mixFrac,
const Composition& fuelComp,
844 vector<double> fuel(getCompositionFromMap(fuelComp));
845 vector<double> ox(getCompositionFromMap(oxComp));
846 setMixtureFraction(mixFrac, fuel.data(), ox.data(), basis);
849 void ThermoPhase::setMixtureFraction(
double mixFrac,
const string& fuelComp,
852 setMixtureFraction(mixFrac,
858 void ThermoPhase::setMixtureFraction(
double mixFrac,
const double* fuelComp,
861 if (mixFrac < 0.0 || mixFrac > 1.0) {
863 "Mixture fraction must be between 0 and 1");
866 vector<double> fuel, ox;
867 if (basis == ThermoBasis::molar) {
870 moleFractionsToMassFractions(fuelComp, fuel.data());
871 moleFractionsToMassFractions(oxComp, ox.data());
872 fuelComp = fuel.data();
876 double sum_yf = std::accumulate(fuelComp, fuelComp+m_kk, 0.0);
877 double sum_yo = std::accumulate(oxComp, oxComp+m_kk, 0.0);
879 if (sum_yf == 0.0 || sum_yo == 0.0) {
881 "No fuel and/or oxidizer composition specified");
884 double p = pressure();
886 vector<double> y(m_kk);
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;
892 setMassFractions_NoNorm(y.data());
899 const string& element)
const
901 vector<double> fuel(getCompositionFromMap(fuelComp));
902 vector<double> ox(getCompositionFromMap(oxComp));
903 return mixtureFraction(fuel.data(), ox.data(), basis, element);
906 double ThermoPhase::mixtureFraction(
const string& fuelComp,
const string& oxComp,
909 return mixtureFraction(
915 double ThermoPhase::mixtureFraction(
const double* fuelComp,
const double* oxComp,
918 vector<double> fuel, ox;
919 if (basis == ThermoBasis::molar) {
922 moleFractionsToMassFractions(fuelComp, fuel.data());
923 moleFractionsToMassFractions(oxComp, ox.data());
924 fuelComp = fuel.data();
928 if (element ==
"Bilger")
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());
934 if (o2_required_fuel < 0.0 || o2_required_ox > 0.0) {
936 "Fuel composition contains too much oxygen or "
937 "oxidizer contains not enough oxygen. "
938 "Fuel and oxidizer composition mixed up?");
941 double denominator = o2_required_fuel - o2_required_ox;
943 if (denominator == 0.0) {
945 "Fuel and oxidizer have the same composition");
948 double Z = (o2_required_mix - o2_required_ox) / denominator;
950 return std::min(std::max(Z, 0.0), 1.0);
953 double sum_yf = std::accumulate(fuelComp, fuelComp+m_kk, 0.0);
954 double sum_yo = std::accumulate(oxComp, oxComp+m_kk, 0.0);
956 if (sum_yf == 0.0 || sum_yo == 0.0) {
958 "No fuel and/or oxidizer composition specified");
961 auto elementalFraction = [
this](
size_t m,
const double* y) {
963 for (
size_t k = 0; k != m_kk; ++k) {
964 Z_m += y[k] / molecularWeight(k) * nAtoms(k, m);
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());
974 if (Z_m_fuel == Z_m_ox) {
976 "Fuel and oxidizer have the same composition for element {}",
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);
995 void ThermoPhase::initThermoFile(
const string& inputFile,
const string&
id)
997 if (inputFile.empty()) {
1001 size_t dot = inputFile.find_last_of(
".");
1004 extension = inputFile.substr(
dot+1);
1006 if (extension ==
"xml" || extension ==
"cti") {
1008 "The CTI and XML formats are no longer supported.");
1011 AnyMap root = AnyMap::fromYamlFile(inputFile);
1012 auto& phase = root[
"phases"].getMapWhere(
"name",
id);
1016 void ThermoPhase::initThermo()
1019 if (!m_spthermo.ready(m_kk)) {
1021 "Missing species thermo data");
1025 void ThermoPhase::setState_TPQ(
double T,
double P,
double Q)
1027 if (T > critTemperature()) {
1028 if (P > critPressure() || Q == 1) {
1033 "Temperature ({}), pressure ({}) and vapor fraction ({}) "
1034 "are inconsistent, above the critical temperature.",
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)) {
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)",
1054 bool ThermoPhase::addSpecies(shared_ptr<Species> spec)
1056 if (!spec->thermo) {
1058 "Species {} has no thermo data", spec->name);
1060 bool added = Phase::addSpecies(spec);
1062 spec->thermo->validate(spec->name);
1063 m_spthermo.install_STIT(m_kk-1, spec->thermo);
1068 void ThermoPhase::modifySpecies(
size_t k, shared_ptr<Species> spec)
1070 if (!spec->thermo) {
1072 "Species {} has no thermo data", spec->name);
1074 Phase::modifySpecies(k, spec);
1075 if (speciesName(k) != spec->name) {
1077 "New species '{}' does not match existing species '{}' at index {}",
1078 spec->name, speciesName(k), k);
1080 spec->thermo->validate(spec->name);
1081 m_spthermo.modifySpecies(k, spec->thermo);
1084 void ThermoPhase::setParameters(
const AnyMap& phaseNode,
const AnyMap& rootNode)
1086 m_input = phaseNode;
1089 AnyMap ThermoPhase::parameters(
bool withInput)
const
1099 void ThermoPhase::getParameters(
AnyMap& phaseNode)
const
1101 phaseNode[
"name"] = name();
1102 phaseNode[
"thermo"] = ThermoFactory::factory()->canonicalize(type());
1104 for (
size_t i = 0; i < nElements(); i++) {
1108 phaseNode[
"species"] = speciesNames();
1111 auto stateVars = nativeState();
1112 if (stateVars.count(
"T")) {
1113 state[
"T"].setQuantity(temperature(),
"K");
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");
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);
1127 Y[speciesName(k)] = Yk;
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);
1137 X[speciesName(k)] = Xk;
1144 phaseNode[
"state"] = std::move(state);
1146 static bool reg = AnyMap::addOrderingRules(
"Phase", {{
"tail",
"state"}});
1148 phaseNode[
"__type__"] =
"Phase";
1157 AnyMap& ThermoPhase::input()
1162 void ThermoPhase::invalidateCache() {
1163 Phase::invalidateCache();
1167 void ThermoPhase::equilibrate(
const string& XY,
const string& solver,
1168 double rtol,
int max_steps,
int max_iter,
1169 int estimate_equil,
int log_level)
1171 if (solver ==
"auto" || solver ==
"element_potential") {
1172 vector<double> initial_state;
1173 saveState(initial_state);
1174 debuglog(
"Trying ChemEquil solver\n", log_level);
1179 int ret = E.
equilibrate(*
this, XY.c_str(), log_level-1);
1182 "ChemEquil solver failed. Return code: {}", ret);
1184 debuglog(
"ChemEquil solver succeeded\n", log_level);
1186 }
catch (std::exception& err) {
1187 debuglog(
"ChemEquil solver failed.\n", log_level);
1189 restoreState(initial_state);
1190 if (solver ==
"auto") {
1197 if (solver ==
"auto" || solver ==
"vcs" || solver ==
"gibbs") {
1201 M.
equilibrate(XY, solver, rtol, max_steps, max_iter,
1202 estimate_equil, log_level);
1206 if (solver !=
"auto") {
1208 "Invalid solver specified: '{}'", solver);
1212 void ThermoPhase::getdlnActCoeffdlnN(
const size_t ld,
double*
const dlnActCoeffdlnN)
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;
1222 void ThermoPhase::getdlnActCoeffdlnN_numderiv(
const size_t ld,
1223 double*
const dlnActCoeffdlnN)
1225 double deltaMoles_j = 0.0;
1226 double pres = pressure();
1229 vector<double> ActCoeff_Base(m_kk);
1230 getActivityCoefficients(ActCoeff_Base.data());
1231 vector<double> Xmol_Base(m_kk);
1232 getMoleFractions(Xmol_Base.data());
1235 vector<double> ActCoeff(m_kk);
1236 vector<double> Xmol(m_kk);
1237 double v_totalMoles = 1.0;
1238 double TMoles_base = v_totalMoles;
1241 for (
size_t j = 0; j < m_kk; j++) {
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;
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;
1256 Xmol[j] = (moles_j_base + deltaMoles_j) / v_totalMoles;
1259 setMoleFractions(Xmol.data());
1261 getActivityCoefficients(ActCoeff.data());
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);
1270 v_totalMoles = TMoles_base;
1273 setMoleFractions(Xmol_Base.data());
1277 string ThermoPhase::report(
bool show_thermo,
double threshold)
const
1279 if (type() ==
"none") {
1281 "Not implemented for thermo model 'none'");
1284 fmt::memory_buffer b;
1286 int name_width = 18;
1288 string blank_leader = fmt::format(
"{:{}}",
"", name_width);
1290 string string_property = fmt::format(
"{{:>{}}} {{}}\n", name_width);
1292 string one_property = fmt::format(
"{{:>{}}} {{:<.5g}} {{}}\n", name_width);
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"
1298 string Y_X_header = fmt::format(
1299 two_prop_header, blank_leader,
"mass frac. Y",
"mole frac. X"
1301 string two_prop_sep = fmt::format(
1302 "{} {:-^15} {:-^15}\n", blank_leader,
"",
""
1304 string two_property = fmt::format(
1305 "{{:>{}}} {{:15.5g}} {{:15.5g}} {{}}\n", name_width
1308 string three_prop_header = fmt::format(
1309 "{} {:^15} {:^15} {:^15}\n", blank_leader,
"mass frac. Y",
1310 "mole frac. X",
"chem. pot. / RT"
1312 string three_prop_sep = fmt::format(
1313 "{} {:-^15} {:-^15} {:-^15}\n", blank_leader,
"",
"",
""
1315 string three_property = fmt::format(
1316 "{{:>{}}} {{:15.5g}} {{:15.5g}} {{:15.5g}}\n", name_width
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");
1328 "mean mol. weight", meanMolecularWeight(),
"kg/kmol");
1330 double phi = electricPotential();
1332 fmt_append(b, one_property,
"potential", phi,
"V");
1335 fmt_append(b, string_property,
"phase of matter", phaseOfMatter());
1342 "enthalpy", enthalpy_mass(), enthalpy_mole(),
"J");
1344 "internal energy", intEnergy_mass(), intEnergy_mole(),
"J");
1346 "entropy", entropy_mass(), entropy_mole(),
"J/K");
1348 "Gibbs function", gibbs_mass(), gibbs_mole(),
"J");
1350 "heat capacity c_p", cp_mass(), cp_mole(),
"J/K");
1353 "heat capacity c_v", cv_mass(), cv_mole(),
"J/K");
1356 "heat capacity c_v",
"<not implemented>");
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]);
1367 double xMinor = 0.0;
1368 double yMinor = 0.0;
1373 for (
size_t k = 0; k < m_kk; k++) {
1374 if (abs(x[k]) >= threshold) {
1377 speciesName(k), y[k], x[k], mu[k]/RT());
1379 fmt_append(b, two_property, speciesName(k), y[k], x[k],
"");
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],
"");
1401 string minor = fmt::format(
"[{:+5d} minor]", nMinor);
1402 fmt_append(b, two_property, minor, yMinor, xMinor,
"");
1405 return to_string(b) + err.
what();
1407 return to_string(b);
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.
size_t size() const
Returns the number of elements in this map.
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
double convert(const string &key, const string &units) const
Convert the item stored by the given key to the units specified in units.
void setFlowStyle(bool flow=true)
Use "flow" style when outputting this AnyMap to YAML.
void erase(const string &key)
Erase the value held by key.
void update(const AnyMap &other, bool keepExisting=true)
Add items from other to this AnyMap.
string keys_str() const
Return a string listing the keys in this AnyMap, for use in error messages, for example.
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.
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.
double relTolerance
Relative tolerance.
int maxIterations
Maximum number of iterations.
A class for multiphase mixtures.
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.
An error indicating that an unimplemented function has been called.
A representation of the units associated with a dimensional quantity.
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....
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 MultiPhase object.
void debuglog(const string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
double dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
T clip(const T &value, const T &lower, const T &upper)
Clip value such that lower <= value <= upper.
const double Faraday
Faraday constant [C/kmol].
const double OneAtm
One atmosphere [Pa].
void setupPhase(ThermoPhase &thermo, const AnyMap &phaseNode, const AnyMap &rootNode)
Initialize a ThermoPhase object.
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
const vector< string > & elementNames()
Get a vector of the names of the elements defined in Cantera.
const int cAC_CONVENTION_MOLAR
Standard state uses the molar convention.
const double SmallNumber
smallest number to compare to zero.
ThermoBasis
Differentiate between mole fractions and mass fractions for input mixture composition.
map< string, double > Composition
Map from string names to doubles.
Contains declarations for string manipulation functions within Cantera.