22 string flag = string(xy);
25 }
else if (flag ==
"TV") {
27 }
else if (flag ==
"HP") {
29 }
else if (flag ==
"UV") {
31 }
else if (flag ==
"SP") {
33 }
else if (flag ==
"SV") {
35 }
else if (flag ==
"UP") {
38 throw CanteraError(
"_equilflag",
"unknown property pair "+flag);
61 m_jwork1.resize(
m_mm+2);
62 m_jwork2.resize(
m_mm+2);
63 m_startSoln.resize(
m_mm+1);
68 m_orderVectorElements.resize(
m_mm);
70 for (
size_t m = 0; m <
m_mm; m++) {
71 m_orderVectorElements[m] = m;
73 m_orderVectorSpecies.resize(
m_kk);
74 for (
size_t k = 0; k <
m_kk; k++) {
75 m_orderVectorSpecies[k] = k;
80 for (
size_t m = 0; m <
m_mm; m++) {
81 for (
size_t k = 0; k <
m_kk; k++) {
89 if (mneg !=
npos && mneg != m) {
91 "negative atom numbers allowed for only one element");
99 "species {} has {} atoms of element {}, "
100 "but this element is not an electron.",
109 for (
size_t k = 0; k <
m_kk; k++) {
110 for (
size_t m = 0; m <
m_mm; m++) {
117 const vector<double>& lambda_RT,
double t)
120 fill(m_mu_RT.begin(), m_mu_RT.end(), 0.0);
121 for (
size_t k = 0; k <
m_kk; k++) {
122 for (
size_t m = 0; m <
m_mm; m++) {
123 m_mu_RT[k] += lambda_RT[m]*
nAtoms(k,m);
146 for (
size_t m = 0; m <
m_mm; m++) {
148 for (
size_t k = 0; k <
m_kk; k++) {
152 "negative mole fraction for {}: {}",
161 for (
size_t m = 0; m <
m_mm; m++) {
172 e.setInitialMixMoles(loglevel-1);
177 m_component[m] = e.componentIndex(m);
185 writelog(
"setInitialMoles: Estimated Mole Fractions\n");
188 for (
size_t k = 0; k <
m_kk; k++) {
192 writelog(
" Element_Name ElementGoal ElementMF\n");
193 for (
size_t m = 0; m <
m_mm; m++) {
202 vector<double>& elMolesGoal,
int loglevel)
204 vector<double> b(
m_mm, -999.0);
205 vector<double> mu_RT(
m_kk, 0.0);
206 vector<double> xMF_est(
m_kk, 0.0);
209 for (
size_t n = 0; n < s.
nSpecies(); n++) {
210 xMF_est[n] = std::max(xMF_est[n], 1e-20);
218 int usedZeroedSpecies = 0;
219 vector<double> formRxnMatrix;
221 &mp, m_orderVectorSpecies,
222 m_orderVectorElements, formRxnMatrix);
225 size_t k = m_orderVectorSpecies[m];
227 xMF_est[k] = std::max(xMF_est[k], 1e-8);
233 m_orderVectorSpecies, m_orderVectorElements);
236 scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(),
241 size_t isp = m_component[m];
247 for (
size_t n = 0; n < s.
nSpecies(); n++) {
255 aa(m,n) =
nAtoms(m_component[m], m_orderVectorElements[n]);
257 b[m] = mu_RT[m_component[m]];
267 lambda_RT[m_orderVectorElements[m]] = b[m];
270 lambda_RT[m_orderVectorElements[m]] = 0.0;
274 writelog(
" id CompSpecies ChemPot EstChemPot Diff\n");
276 size_t isp = m_component[m];
278 for (
size_t n = 0; n <
m_mm; n++) {
279 tmp +=
nAtoms(isp, n) * lambda_RT[n];
281 writelogf(
"%3d %16s %10.5g %10.5g %10.5g\n",
282 m, s.
speciesName(isp), mu_RT[isp], tmp, tmp - mu_RT[isp]);
286 for (
size_t m = 0; m <
m_mm; m++) {
298 return equilibrate(s, XY, elMolesGoal, loglevel-1);
302 vector<double>& elMolesGoal,
int loglevel)
305 bool tempFixed =
true;
307 vector<double> state;
314 "Input ThermoPhase is incompatible with initialization");
356 "illegal property pair '{}'", XYstr);
363 throw CanteraError(
"ChemEquil::equilibrate",
"Specified temperature"
364 " ({} K) outside valid range of {} K to {} K\n",
371 double xval = m_p1(s);
372 double yval = m_p2(s);
375 size_t nvar = mm + 1;
377 vector<double> x(nvar, -102.0);
378 vector<double> res_trial(nvar, 0.0);
386 size_t m = m_orderVectorElements[im];
387 if (elMolesGoal[m] > tmp) {
389 tmp = elMolesGoal[m];
394 "Element Abundance Vector is zeroed");
400 vector<double> xmm(
m_kk, 0.0);
401 for (
size_t k = 0; k <
m_kk; k++) {
410 double tmaxPhase = s.
maxTemp();
411 double tminPhase = s.
minTemp();
414 double tmin = std::max(s.
temperature(), tminPhase);
415 if (tmin > tmaxPhase) {
416 tmin = tmaxPhase - 20;
418 double tmax = std::min(tmin + 10., tmaxPhase);
419 if (tmax < tminPhase) {
420 tmax = tminPhase + 20;
423 double slope, phigh, plow, pval, dt;
438 double t0 = 0.5*(tmin + tmax);
442 for (
int it = 0; it < 10; it++) {
461 slope = (phigh - plow)/(tmax - tmin);
462 dt = (xval - pval)/slope;
465 if (fabs(dt) < 50.0) {
468 dt =
clip(dt, -200.0, 200.0);
469 if ((t0 + dt) < tminPhase) {
470 dt = 0.5*((t0) + tminPhase) - t0;
472 if ((t0 + dt) > tmaxPhase) {
473 dt = 0.5*((t0) + tmaxPhase) - t0;
477 if (t0 <= tminPhase || t0 >= tmaxPhase || t0 < 100.0) {
478 throw CanteraError(
"ChemEquil::equilibrate",
"T out of bounds");
510 vector<double> above(nvar);
511 vector<double> below(nvar);
512 for (
size_t m = 0; m < mm; m++) {
522 above[mm] = log(s.
maxTemp() + 25.0);
523 below[mm] = log(s.
minTemp() - 25.0);
525 vector<double> grad(nvar, 0.0);
526 vector<double> oldx(nvar, 0.0);
527 vector<double> oldresid(nvar, 0.0);
532 double f = 0.5*
dot(res_trial.begin(), res_trial.end(), res_trial.begin());
535 double deltax = (xx - xval)/xval;
536 double deltay = (yy - yval)/yval;
537 bool passThis =
true;
538 for (
size_t m = 0; m < nvar; m++) {
559 if (fabs(res_trial[m]) > tval) {
568 adjustEloc(s, elMolesGoal);
574 "Temperature ({} K) outside valid range of {} K "
582 f = 0.5*
dot(res_trial.begin(), res_trial.end(), res_trial.begin());
585 equilJacobian(s, x, elMolesGoal, jac, xval, yval);
588 writelogf(
"Jacobian matrix %d:\n", iter);
589 for (
size_t m = 0; m <=
m_mm; m++) {
591 for (
size_t n = 0; n <=
m_mm; n++) {
599 }
else if (m == m_skip) {
604 writelog(
" = - ({:10.5g})\n", res_trial[m]);
610 scale(res_trial.begin(), res_trial.end(), res_trial.begin(), -1.0);
618 "Jacobian is singular. \nTry adding more species, "
619 "changing the elemental composition slightly, \nor removing "
626 for (
size_t m = 0; m < nvar; m++) {
627 double newval = x[m] + res_trial[m];
628 if (newval > above[m]) {
630 std::min(fctr,0.8*(above[m] - x[m])/(newval - x[m])));
631 }
else if (newval < below[m]) {
632 if (m <
m_mm && (m != m_skip)) {
634 if (x[m] < below[m] + 50.) {
635 res_trial[m] = below[m] - x[m];
638 fctr = std::min(fctr, 0.8*(x[m] - below[m])/(x[m] - newval));
642 if (m == mm && fabs(res_trial[mm]) > 0.2) {
643 fctr = std::min(fctr, 0.2/fabs(res_trial[mm]));
646 if (fctr != 1.0 && loglevel > 0) {
648 "Soln Damping because of bounds: %g", fctr);
652 scale(res_trial.begin(), res_trial.end(), res_trial.begin(), fctr);
654 if (!
dampStep(s, oldx, oldf, grad, res_trial,
655 x, f, elMolesGoal , xval, yval)) {
660 "Cannot find an acceptable Newton damping coefficient.");
675 vector<double>& grad, vector<double>& step, vector<double>& x,
676 double& f, vector<double>& elmols,
double xval,
double yval)
681 for (
size_t m = 0; m <
m_mm; m++) {
683 if (step[m] > 1.25) {
684 damp = std::min(damp, 1.25 /step[m]);
686 if (step[m] < -1.25) {
687 damp = std::min(damp, -1.25 / step[m]);
690 if (step[m] > 0.75) {
691 damp = std::min(damp, 0.75 /step[m]);
693 if (step[m] < -0.75) {
694 damp = std::min(damp, -0.75 / step[m]);
700 for (
size_t m = 0; m < x.size(); m++) {
701 x[m] = oldx[m] + damp * step[m];
704 writelogf(
"Solution Unknowns: damp = %g\n", damp);
706 for (
size_t m = 0; m <
m_mm; m++) {
707 writelogf(
" % -10.5g % -10.5g % -10.5g\n", x[m], oldx[m], step[m]);
714 const vector<double>& elmFracGoal, vector<double>& resid,
715 double xval,
double yval,
int loglevel)
721 for (
size_t n = 0; n <
m_mm; n++) {
722 size_t m = m_orderVectorElements[n];
725 resid[m] = x[m] + 1000.0;
731 if (elmFracGoal[m] < 1.0E-10 || elmFrac[m] < 1.0E-10 || m ==
m_eloc) {
732 resid[m] = elmFracGoal[m] - elmFrac[m];
734 resid[m] = log((1.0 + elmFracGoal[m]) / (1.0 + elmFrac[m]));
739 if (loglevel > 0 && !m_doResPerturb) {
740 writelog(
"Residual: ElFracGoal ElFracCurrent Resid\n");
741 for (
size_t n = 0; n <
m_mm; n++) {
742 writelogf(
" % -14.7E % -14.7E % -10.5E\n",
743 elmFracGoal[n], elmFrac[n], resid[n]);
749 resid[
m_mm] = xx/xval - 1.0;
750 resid[m_skip] = yy/yval - 1.0;
752 if (loglevel > 0 && !m_doResPerturb) {
754 writelogf(
" XX : % -14.7E % -14.7E % -10.5E\n", xval, xx, resid[
m_mm]);
755 writelogf(
" YY(%1d): % -14.7E % -14.7E % -10.5E\n", m_skip, yval, yy, resid[m_skip]);
759void ChemEquil::equilJacobian(
ThermoPhase& s, vector<double>& x,
761 double xval,
double yval,
int loglevel)
763 vector<double>& r0 = m_jwork1;
764 vector<double>& r1 = m_jwork2;
765 size_t len = x.size();
768 double atol = 1.e-10;
772 m_doResPerturb =
false;
773 for (
size_t n = 0; n <
len; n++) {
775 double dx = std::max(atol, fabs(xsave) * 1.0E-7);
784 for (
size_t m = 0; m < x.size(); m++) {
785 jac(m, n) = (r1[m] - r0[m])*rdx;
789 m_doResPerturb =
false;
793 const vector<double>& Xmol_i_calc,
794 vector<double>& eMolesCalc, vector<double>& n_i_calc,
795 double pressureConst)
797 double n_t_calc = 0.0;
801 vector<double> actCoeff(
m_kk, 1.0);
806 for (
size_t k = 0; k <
m_kk; k++) {
807 double tmp = - (
m_muSS_RT[k] + log(actCoeff[k]));
808 for (
size_t m = 0; m <
m_mm; m++) {
809 tmp +=
nAtoms(k,m) * x[m];
811 tmp = std::min(tmp, 100.0);
815 n_i_calc[k] = n_t * exp(tmp);
817 n_t_calc += n_i_calc[k];
819 for (
size_t m = 0; m <
m_mm; m++) {
821 for (
size_t k = 0; k <
m_kk; k++) {
822 eMolesCalc[m] +=
nAtoms(k,m) * n_i_calc[k];
829 vector<double>& elMoles)
833 vector<double> state;
835 bool modifiedMatrix =
false;
839 vector<double> b(neq, 0.0);
840 vector<double> n_i(
m_kk,0.0);
841 vector<double> n_i_calc(
m_kk,0.0);
842 vector<double> actCoeff(
m_kk, 1.0);
846 double pressureConst = s.
pressure();
847 vector<double> Xmol_i_calc = n_i;
849 vector<double> x_old(
m_mm+1, 0.0);
850 vector<double> resid(
m_mm+1, 0.0);
851 vector<int> lumpSum(
m_mm+1, 0);
857 vector<double> eMolesCalc(
m_mm, 0.0);
858 vector<double> eMolesFix(
m_mm, 0.0);
859 double elMolesTotal = 0.0;
860 for (
size_t m = 0; m <
m_mm; m++) {
861 elMolesTotal += elMoles[m];
862 for (
size_t k = 0; k <
m_kk; k++) {
863 eMolesFix[m] +=
nAtoms(k,m) * n_i[k];
867 for (
size_t m = 0; m <
m_mm; m++) {
868 if (elMoles[m] > 1.0E-70) {
869 x[m] =
clip(x[m], -100.0, 50.0);
871 x[m] =
clip(x[m], -1000.0, 50.0);
876 double nAtomsMax = 1.0;
880 for (
size_t k = 0; k <
m_kk; k++) {
881 double tmp = - (
m_muSS_RT[k] + log(actCoeff[k]));
883 for (
size_t m = 0; m <
m_mm; m++) {
887 nAtomsMax = std::max(nAtomsMax, sum2);
897 writelog(
"estimateEP_Brinkley::\n\n");
900 writelog(
"Initial mole numbers and mu_SS:\n");
901 writelog(
" Name MoleNum mu_SS actCoeff\n");
902 for (
size_t k = 0; k <
m_kk; k++) {
906 writelogf(
"Initial n_t = %10.5g\n", n_t);
907 writelog(
"Comparison of Goal Element Abundance with Initial Guess:\n");
908 writelog(
" eName eCurrent eGoal\n");
909 for (
size_t m = 0; m <
m_mm; m++) {
914 for (
size_t m = 0; m <
m_mm; m++) {
923 for (
size_t m = 0; m <
m_mm; m++) {
929 writelogf(
"START ITERATION %d:\n", iter);
932 double n_t_calc =
calcEmoles(s, x, n_t, Xmol_i_calc, eMolesCalc, n_i_calc,
935 for (
size_t k = 0; k <
m_kk; k++) {
936 Xmol_i_calc[k] = n_i_calc[k]/n_t_calc;
940 writelog(
" Species: Calculated_Moles Calculated_Mole_Fraction\n");
941 for (
size_t k = 0; k <
m_kk; k++) {
945 writelogf(
"%15s: %10.5g\n",
"Total Molar Sum", n_t_calc);
946 writelogf(
"(iter %d) element moles bal: Goal Calculated\n", iter);
947 for (
size_t m = 0; m <
m_mm; m++) {
953 bool normalStep =
true;
956 for (
size_t m = 0; m <
m_mm; m++) {
957 if (elMoles[m] > 0.001 * elMolesTotal) {
958 if (eMolesCalc[m] > 1000. * elMoles[m]) {
962 if (1000 * eMolesCalc[m] < elMoles[m]) {
969 writelogf(
" NOTE: iter(%d) Doing an abnormal step due to row %d\n", iter, iM);
974 for (
size_t im = 0; im <
m_mm; im++) {
975 size_t m = m_orderVectorElements[im];
977 if (im <
m_nComponents && elMoles[m] > 0.001 * elMolesTotal) {
978 if (eMolesCalc[m] > 1000. * elMoles[m]) {
982 if (1000 * eMolesCalc[m] < elMoles[m]) {
988 if (n_t < (elMolesTotal / nAtomsMax)) {
989 if (resid[
m_mm] < 0.0) {
992 }
else if (n_t > elMolesTotal) {
993 resid[
m_mm] = std::min(resid[
m_mm], 0.0);
1018 for (
size_t m = 0; m <
m_mm; m++) {
1022 double nCutoff = 1.0E-9 * n_t_calc;
1024 writelog(
" Lump Sum Elements Calculation: \n");
1026 for (
size_t m = 0; m <
m_mm; m++) {
1028 size_t kMSp2 =
npos;
1029 for (
size_t k = 0; k <
m_kk; k++) {
1030 if (n_i_calc[k] > nCutoff && fabs(
nAtoms(k,m)) > 0.001) {
1033 double factor = fabs(
nAtoms(kMSp,m) /
nAtoms(kMSp2,m));
1034 for (
size_t n = 0; n <
m_mm; n++) {
1035 if (fabs(factor *
nAtoms(kMSp2,n) -
nAtoms(kMSp,n)) > 1.0E-8) {
1052 for (
size_t im = 0; im <
m_mm; im++) {
1053 size_t m = m_orderVectorElements[im];
1055 for (
size_t n = 0; n <
m_mm; n++) {
1057 for (
size_t k = 0; k <
m_kk; k++) {
1061 a1(m,
m_mm) = eMolesCalc[m];
1062 a1(
m_mm, m) = eMolesCalc[m];
1064 for (
size_t n = 0; n <=
m_mm; n++) {
1075 for (
size_t im = 0; im <
m_mm; im++) {
1076 size_t m = m_orderVectorElements[im];
1078 resid[m] = elMoles[m] - eMolesCalc[m];
1097 for (
size_t m = 0; m <
m_mm; m++) {
1098 if (a1(m,m) < 1.0E-50) {
1100 writelogf(
" NOTE: Diagonalizing the analytical Jac row %d\n", m);
1102 for (
size_t n = 0; n <
m_mm; n++) {
1106 if (resid[m] > 0.0) {
1108 }
else if (resid[m] < 0.0) {
1116 resid[
m_mm] = n_t - n_t_calc;
1120 for (
size_t m = 0; m <=
m_mm; m++) {
1122 for (
size_t n = 0; n <=
m_mm; n++) {
1129 sum += pow(resid[
m_mm] /(n_t + 1.0E-15), 2);
1131 writelogf(
"(it %d) Convergence = %g\n", iter, sum);
1144 for (
size_t m = 0; m <=
m_mm; m++) {
1146 for (
size_t n = 0; n <=
m_mm; n++) {
1147 tmp += fabs(a1(m,n));
1149 if (m <
m_mm && tmp < 1.0E-30) {
1151 writelogf(
" NOTE: Diagonalizing row %d\n", m);
1153 for (
size_t n = 0; n <=
m_mm; n++) {
1161 for (
size_t n = 0; n <=
m_mm; n++) {
1169 for (
size_t m = 0; m <=
m_mm; m++) {
1171 for (
size_t n = 0; n <=
m_mm; n++) {
1197 modifiedMatrix =
false;
1198 for (
size_t m = 0; m <
m_mm; m++) {
1199 size_t sameAsRow =
npos;
1200 for (
size_t im = 0; im < m; im++) {
1201 bool theSame =
true;
1202 for (
size_t n = 0; n <
m_mm; n++) {
1203 if (fabs(a1(m,n) - a1(im,n)) > 1.0E-7) {
1212 if (sameAsRow !=
npos || lumpSum[m]) {
1215 writelogf(
"Lump summing row %d, due to rank deficiency analysis\n", m);
1216 }
else if (sameAsRow !=
npos) {
1217 writelogf(
"Identified that rows %d and %d are the same\n", m, sameAsRow);
1220 modifiedMatrix =
true;
1221 for (
size_t n = 0; n <
m_mm; n++) {
1223 a1(m,m) += fabs(a1(m,n));
1231 writelog(
"Row Summed, MODIFIED Matrix:\n");
1232 for (
size_t m = 0; m <=
m_mm; m++) {
1234 for (
size_t n = 0; n <=
m_mm; n++) {
1246 "Jacobian is singular. \nTry adding more species, "
1247 "changing the elemental composition slightly, \nor removing "
1254 for (
size_t m = 0; m <
m_mm; m++) {
1255 if (resid[m] > 1.0) {
1256 beta = std::min(beta, 1.0 / resid[m]);
1258 if (resid[m] < -1.0) {
1259 beta = std::min(beta, -1.0 / resid[m]);
1263 writelogf(
"(it %d) Beta = %g\n", iter, beta);
1267 for (
size_t m = 0; m <
m_mm; m++) {
1268 x[m] += beta * resid[m];
1270 n_t *= exp(beta * resid[
m_mm]);
1273 writelogf(
"(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter);
1274 for (
size_t m = 0; m <
m_mm; m++) {
1275 writelogf(
" %5s %10.5g %10.5g %10.5g\n",
1286 writelogf(
" ChemEquil::estimateEP_Brinkley() SUCCESS: equilibrium found at T = %g, Pres = %g\n",
1289 writelogf(
" ChemEquil::estimateEP_Brinkley() FAILURE: equilibrium not found at T = %g, Pres = %g\n",
1297void ChemEquil::adjustEloc(
ThermoPhase& s, vector<double>& elMolesGoal)
1302 if (fabs(elMolesGoal[
m_eloc]) > 1.0E-20) {
1306 size_t maxPosEloc =
npos;
1307 size_t maxNegEloc =
npos;
1308 double maxPosVal = -1.0;
1309 double maxNegVal = -1.0;
1311 for (
size_t k = 0; k <
m_kk; k++) {
1323 double sumPos = 0.0;
1324 double sumNeg = 0.0;
1325 for (
size_t k = 0; k <
m_kk; k++) {
1335 if (sumPos >= sumNeg) {
1336 if (sumPos <= 0.0) {
1339 double factor = (elMolesGoal[
m_eloc] + sumNeg) / sumPos;
1340 if (
m_loglevel > 0 && factor < 0.9999999999) {
1341 writelogf(
"adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n",
1345 for (
size_t k = 0; k <
m_kk; k++) {
1351 double factor = (-elMolesGoal[
m_eloc] + sumPos) / sumNeg;
1352 if (
m_loglevel > 0 && factor < 0.9999999999) {
1353 writelogf(
"adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n",
1357 for (
size_t k = 0; k <
m_kk; k++) {
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
vector< double > & data()
Return a reference to the data vector.
Base class for exceptions thrown by Cantera classes.
virtual string getMessage() const
Method overridden by derived classes to format the error message.
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...
size_t m_kk
number of species in the phase
int m_loglevel
Verbosity of printed output.
size_t m_nComponents
This is equal to the rank of the stoichiometric coefficient matrix when it is computed.
int setInitialMoles(ThermoPhase &s, vector< double > &elMoleGoal, int loglevel=0)
Estimate the initial mole numbers.
ThermoPhase * m_phase
Pointer to the ThermoPhase object used to initialize this object.
double m_elementTotalSum
Current value of the sum of the element abundances given the current element potentials.
void update(const ThermoPhase &s)
Update internally stored state information.
void setToEquilState(ThermoPhase &s, const vector< double > &x, double t)
Set mixture to an equilibrium state consistent with specified element potentials and temperature.
size_t m_eloc
Index of the element id corresponding to the electric charge of each species.
int estimateElementPotentials(ThermoPhase &s, vector< double > &lambda, vector< double > &elMolesGoal, int loglevel=0)
Generate a starting estimate for the element potentials.
int estimateEP_Brinkley(ThermoPhase &s, vector< double > &lambda, vector< double > &elMoles)
Do a calculation of the element potentials using the Brinkley method, p.
void initialize(ThermoPhase &s)
Prepare for equilibrium calculations.
EquilOpt options
Options controlling how the calculation is carried out.
vector< double > m_molefractions
Current value of the mole fractions in the single phase. length = m_kk.
int dampStep(ThermoPhase &s, vector< double > &oldx, double oldf, vector< double > &grad, vector< double > &step, vector< double > &x, double &f, vector< double > &elmols, double xval, double yval)
Find an acceptable step size and take it.
double calcEmoles(ThermoPhase &s, vector< double > &x, const double &n_t, const vector< double > &Xmol_i_calc, vector< double > &eMolesCalc, vector< double > &n_i_calc, double pressureConst)
Given a vector of dimensionless element abundances, this routine calculates the moles of the elements...
vector< double > m_comp
Storage of the element compositions. natom(k,m) = m_comp[k*m_mm+ m];.
double nAtoms(size_t k, size_t m) const
number of atoms of element m in species k.
void equilResidual(ThermoPhase &s, const vector< double > &x, const vector< double > &elmtotal, vector< double > &resid, double xval, double yval, int loglevel=0)
Evaluates the residual vector F, of length m_mm.
double m_elemFracCutoff
element fractional cutoff, below which the element will be zeroed.
vector< double > m_muSS_RT
Dimensionless values of the Gibbs free energy for the standard state of each species,...
vector< double > m_elementmolefracs
Current value of the element mole fractions.
size_t m_mm
number of elements in the phase
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
double absElemTol
Abs Tol in element number.
int iterations
Iteration counter.
double relTolerance
Relative tolerance.
int maxIterations
Maximum number of iterations.
Multiphase chemical equilibrium solver.
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.
virtual void setMoleFractions(const double *const x)
Set the mole fractions to the specified values.
void restoreState(const vector< double > &state)
Restore a state saved on a previous call to saveState.
size_t nSpecies() const
Returns the number of species in the phase.
double temperature() const
Temperature (K).
virtual void setPressure(double p)
Set the internally stored pressure (Pa) at constant temperature and composition.
void saveState(vector< double > &state) const
Save the current internal state of the phase.
double atomicWeight(size_t m) const
Atomic weight of element m.
string speciesName(size_t k) const
Name of the species with index k.
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
double moleFraction(size_t k) const
Return the mole fraction of a single species.
virtual double density() const
Density (kg/m^3).
double nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
virtual void setTemperature(double temp)
Set the internally stored temperature of the phase (K).
size_t nElements() const
Number of elements.
virtual double pressure() const
Return the thermodynamic pressure (Pa).
string elementName(size_t m) const
Name of the element with index m.
Base class for a phase with thermodynamic properties.
virtual double minTemp(size_t k=npos) const
Minimum temperature for which the thermodynamic data for the species or phase are valid.
virtual double maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
virtual void getGibbs_RT(double *grt) const
Get the nondimensional Gibbs functions for the species in their standard states at the current T and ...
virtual void getActivityCoefficients(double *ac) const
Get the array of non-dimensional molar-based activity coefficients at the current solution temperatur...
double entropy_mass() const
Specific entropy. Units: J/kg/K.
virtual void getChemPotentials(double *mu) const
Get the species chemical potentials. Units: J/kmol.
double intEnergy_mass() const
Specific internal energy. Units: J/kg.
virtual double refPressure() const
Returns the reference pressure in Pa.
double enthalpy_mass() const
Specific enthalpy. Units: J/kg.
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
virtual void setToEquilState(const double *mu_RT)
This method is used by the ChemEquil equilibrium solver.
void ElemRearrange(size_t nComponents, const vector< double > &elementAbundances, MultiPhase *mphase, vector< size_t > &orderVectorSpecies, vector< size_t > &orderVectorElements)
Handles the potential rearrangement of the constraint equations represented by the Formula Matrix.
size_t BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, MultiPhase *mphase, vector< size_t > &orderVectorSpecies, vector< size_t > &orderVectorElements, vector< double > &formRxnMatrix)
Choose the optimum basis of species for the equilibrium calculations.
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
U len(const T &container)
Get the size of a container, cast to a signed integer type.
double dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
T clip(const T &value, const T &lower, const T &upper)
Clip value such that lower <= value <= upper.
const double GasConstant
Universal Gas Constant [J/kmol/K].
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
int _equilflag(const char *xy)
map property strings to integers
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...