23 string flag = string(xy);
26 }
else if (flag ==
"TV") {
28 }
else if (flag ==
"HP") {
30 }
else if (flag ==
"UV") {
32 }
else if (flag ==
"SP") {
34 }
else if (flag ==
"SV") {
36 }
else if (flag ==
"UP") {
39 throw CanteraError(
"_equilflag",
"unknown property pair "+flag);
44 ChemEquil::ChemEquil() : m_skip(
npos), m_elementTotalSum(1.0),
46 m_elemFracCutoff(1.0E-100),
52 m_elementTotalSum(1.0),
54 m_elemFracCutoff(1.0E-100),
60 ChemEquil::~ChemEquil()
78 m_jwork1.resize(
m_mm+2);
79 m_jwork2.resize(
m_mm+2);
80 m_startSoln.resize(
m_mm+1);
85 m_orderVectorElements.resize(
m_mm);
87 for (
size_t m = 0; m <
m_mm; m++) {
88 m_orderVectorElements[m] = m;
90 m_orderVectorSpecies.resize(
m_kk);
91 for (
size_t k = 0; k <
m_kk; k++) {
92 m_orderVectorSpecies[k] = k;
97 for (
size_t m = 0; m <
m_mm; m++) {
98 for (
size_t k = 0; k <
m_kk; k++) {
102 if (s.
nAtoms(k,m) < 0.0) {
106 if (mneg !=
npos && mneg != m) {
108 "negative atom numbers allowed for only one element");
115 writelog(
"WARNING: species {} has {} atoms of element {}," 116 " but this element is not an electron.\n",
125 for (
size_t k = 0; k <
m_kk; k++) {
126 for (
size_t m = 0; m <
m_mm; m++) {
133 const vector_fp& lambda_RT, doublereal t)
136 fill(m_mu_RT.begin(), m_mu_RT.end(), 0.0);
137 for (
size_t k = 0; k <
m_kk; k++) {
138 for (
size_t m = 0; m <
m_mm; m++) {
139 m_mu_RT[k] += lambda_RT[m]*
nAtoms(k,m);
162 for (
size_t m = 0; m <
m_mm; m++) {
164 for (
size_t k = 0; k <
m_kk; k++) {
168 "negative mole fraction for {}: {}",
177 for (
size_t m = 0; m <
m_mm; m++) {
189 e.setInitialMixMoles(loglevel-1);
194 m_component[m] = e.componentIndex(m);
202 writelog(
"setInitialMoles: Estimated Mole Fractions\n");
205 for (
size_t k = 0; k <
m_kk; k++) {
209 writelog(
" Element_Name ElementGoal ElementMF\n");
210 for (
size_t m = 0; m <
m_mm; m++) {
226 for (
size_t n = 0; n < s.
nSpecies(); n++) {
227 xMF_est[n] = std::max(xMF_est[n], 1e-20);
235 int usedZeroedSpecies = 0;
238 &mp, m_orderVectorSpecies,
239 m_orderVectorElements, formRxnMatrix);
242 size_t k = m_orderVectorSpecies[m];
244 xMF_est[k] = std::max(xMF_est[k], 1e-8);
250 m_orderVectorSpecies, m_orderVectorElements);
253 scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(),
258 size_t isp = m_component[m];
264 for (
size_t n = 0; n < s.
nSpecies(); n++) {
272 aa(m,n) =
nAtoms(m_component[m], m_orderVectorElements[n]);
274 b[m] = mu_RT[m_component[m]];
282 lambda_RT[m_orderVectorElements[m]] = b[m];
285 lambda_RT[m_orderVectorElements[m]] = 0.0;
289 writelog(
" id CompSpecies ChemPot EstChemPot Diff\n");
291 size_t isp = m_component[m];
293 for (
size_t n = 0; n <
m_mm; n++) {
294 tmp +=
nAtoms(isp, n) * lambda_RT[n];
296 writelogf(
"%3d %16s %10.5g %10.5g %10.5g\n",
297 m, s.
speciesName(isp), mu_RT[isp], tmp, tmp - mu_RT[isp]);
301 for (
size_t m = 0; m <
m_mm; m++) {
313 return equilibrate(s, XY, elMolesGoal, loglevel-1);
320 bool tempFixed =
true;
329 "Input ThermoPhase is incompatible with initialization");
370 throw CanteraError(
"equilibrate",
"illegal property pair.");
377 throw CanteraError(
"ChemEquil::equilibrate",
"Specified temperature" 378 " ({} K) outside valid range of {} K to {} K\n",
385 double xval = m_p1(s);
386 double yval = m_p2(s);
389 size_t nvar = mm + 1;
400 size_t m = m_orderVectorElements[im];
401 if (elMolesGoal[m] > tmp) {
403 tmp = elMolesGoal[m];
408 "Element Abundance Vector is zeroed");
415 for (
size_t k = 0; k <
m_kk; k++) {
424 doublereal tmaxPhase = s.
maxTemp();
425 doublereal tminPhase = s.
minTemp();
428 doublereal tmin = std::max(s.
temperature(), tminPhase);
429 if (tmin > tmaxPhase) {
430 tmin = tmaxPhase - 20;
432 doublereal tmax = std::min(tmin + 10., tmaxPhase);
433 if (tmax < tminPhase) {
434 tmax = tminPhase + 20;
437 doublereal slope, phigh, plow, pval, dt;
452 doublereal t0 = 0.5*(tmin + tmax);
456 for (
int it = 0; it < 10; it++) {
475 slope = (phigh - plow)/(tmax - tmin);
476 dt = (xval - pval)/slope;
479 if (fabs(dt) < 50.0) {
482 dt =
clip(dt, -200.0, 200.0);
483 if ((t0 + dt) < tminPhase) {
484 dt = 0.5*((t0) + tminPhase) - t0;
486 if ((t0 + dt) > tmaxPhase) {
487 dt = 0.5*((t0) + tmaxPhase) - t0;
491 if (t0 <= tminPhase || t0 >= tmaxPhase || t0 < 100.0) {
492 throw CanteraError(
"ChemEquil::equilibrate",
"T out of bounds");
526 for (
size_t m = 0; m < mm; m++) {
536 above[mm] = log(s.
maxTemp() + 25.0);
537 below[mm] = log(s.
minTemp() - 25.0);
546 double f = 0.5*
dot(res_trial.begin(), res_trial.end(), res_trial.begin());
549 double deltax = (xx - xval)/xval;
550 double deltay = (yy - yval)/yval;
551 bool passThis =
true;
552 for (
size_t m = 0; m < nvar; m++) {
573 if (fabs(res_trial[m]) > tval) {
580 for (
size_t m = 0; m <
m_mm; m++) {
585 adjustEloc(s, elMolesGoal);
590 writelog(
"Warning: Temperature ({} K) outside valid range of " 599 f = 0.5*
dot(res_trial.begin(), res_trial.end(), res_trial.begin());
602 equilJacobian(s, x, elMolesGoal, jac, xval, yval);
605 writelogf(
"Jacobian matrix %d:\n", iter);
606 for (
size_t m = 0; m <=
m_mm; m++) {
608 for (
size_t n = 0; n <=
m_mm; n++) {
616 }
else if (m == m_skip) {
621 writelog(
" = - ({:10.5g})\n", res_trial[m]);
627 scale(res_trial.begin(), res_trial.end(), res_trial.begin(), -1.0);
635 "Jacobian is singular. \nTry adding more species, " 636 "changing the elemental composition slightly, \nor removing " 643 for (
size_t m = 0; m < nvar; m++) {
644 double newval = x[m] + res_trial[m];
645 if (newval > above[m]) {
647 std::min(fctr,0.8*(above[m] - x[m])/(newval - x[m])));
648 }
else if (newval < below[m]) {
649 if (m <
m_mm && (m != m_skip)) {
651 if (x[m] < below[m] + 50.) {
652 res_trial[m] = below[m] - x[m];
655 fctr = std::min(fctr, 0.8*(x[m] - below[m])/(x[m] - newval));
659 if (m == mm && fabs(res_trial[mm]) > 0.2) {
660 fctr = std::min(fctr, 0.2/fabs(res_trial[mm]));
663 if (fctr != 1.0 && loglevel > 0) {
664 writelogf(
"WARNING Soln Damping because of bounds: %g\n", fctr);
668 scale(res_trial.begin(), res_trial.end(), res_trial.begin(), fctr);
670 if (!
dampStep(s, oldx, oldf, grad, res_trial,
671 x, f, elMolesGoal , xval, yval)) {
676 "Cannot find an acceptable Newton damping coefficient.");
692 double& f,
vector_fp& elmols,
double xval,
double yval)
697 for (
size_t m = 0; m <
m_mm; m++) {
699 if (step[m] > 1.25) {
700 damp = std::min(damp, 1.25 /step[m]);
702 if (step[m] < -1.25) {
703 damp = std::min(damp, -1.25 / step[m]);
706 if (step[m] > 0.75) {
707 damp = std::min(damp, 0.75 /step[m]);
709 if (step[m] < -0.75) {
710 damp = std::min(damp, -0.75 / step[m]);
716 for (
size_t m = 0; m < x.size(); m++) {
717 x[m] = oldx[m] + damp * step[m];
720 writelogf(
"Solution Unknowns: damp = %g\n", damp);
722 for (
size_t m = 0; m <
m_mm; m++) {
723 writelogf(
" % -10.5g % -10.5g % -10.5g\n", x[m], oldx[m], step[m]);
731 doublereal xval, doublereal yval,
int loglevel)
737 for (
size_t n = 0; n <
m_mm; n++) {
738 size_t m = m_orderVectorElements[n];
741 resid[m] = x[m] + 1000.0;
747 if (elmFracGoal[m] < 1.0E-10 || elmFrac[m] < 1.0E-10 || m ==
m_eloc) {
748 resid[m] = elmFracGoal[m] - elmFrac[m];
750 resid[m] = log((1.0 + elmFracGoal[m]) / (1.0 + elmFrac[m]));
755 if (loglevel > 0 && !m_doResPerturb) {
756 writelog(
"Residual: ElFracGoal ElFracCurrent Resid\n");
757 for (
size_t n = 0; n <
m_mm; n++) {
758 writelogf(
" % -14.7E % -14.7E % -10.5E\n",
759 elmFracGoal[n], elmFrac[n], resid[n]);
765 resid[
m_mm] = xx/xval - 1.0;
766 resid[m_skip] = yy/yval - 1.0;
768 if (loglevel > 0 && !m_doResPerturb) {
770 writelogf(
" XX : % -14.7E % -14.7E % -10.5E\n", xval, xx, resid[
m_mm]);
771 writelogf(
" YY(%1d): % -14.7E % -14.7E % -10.5E\n", m_skip, yval, yy, resid[m_skip]);
777 doublereal xval, doublereal yval,
int loglevel)
781 size_t len = x.size();
784 doublereal atol = 1.e-10;
788 m_doResPerturb =
false;
789 for (
size_t n = 0; n < len; n++) {
791 double dx = std::max(atol, fabs(xsave) * 1.0E-7);
800 for (
size_t m = 0; m < x.size(); m++) {
801 jac(m, n) = (r1[m] - r0[m])*rdx;
805 m_doResPerturb =
false;
811 double pressureConst)
813 double n_t_calc = 0.0;
822 for (
size_t k = 0; k <
m_kk; k++) {
823 double tmp = - (
m_muSS_RT[k] + log(actCoeff[k]));
824 for (
size_t m = 0; m <
m_mm; m++) {
825 tmp +=
nAtoms(k,m) * x[m];
827 tmp = std::min(tmp, 100.0);
831 n_i_calc[k] = n_t * exp(tmp);
833 n_t_calc += n_i_calc[k];
835 for (
size_t m = 0; m <
m_mm; m++) {
837 for (
size_t k = 0; k <
m_kk; k++) {
838 eMolesCalc[m] +=
nAtoms(k,m) * n_i_calc[k];
851 bool modifiedMatrix =
false;
862 double pressureConst = s.
pressure();
875 double elMolesTotal = 0.0;
876 for (
size_t m = 0; m <
m_mm; m++) {
877 elMolesTotal += elMoles[m];
878 for (
size_t k = 0; k <
m_kk; k++) {
879 eMolesFix[m] +=
nAtoms(k,m) * n_i[k];
883 for (
size_t m = 0; m <
m_mm; m++) {
884 if (elMoles[m] > 1.0E-70) {
885 x[m] =
clip(x[m], -100.0, 50.0);
887 x[m] =
clip(x[m], -1000.0, 50.0);
892 double nAtomsMax = 1.0;
896 for (
size_t k = 0; k <
m_kk; k++) {
897 double tmp = - (
m_muSS_RT[k] + log(actCoeff[k]));
899 for (
size_t m = 0; m <
m_mm; m++) {
903 nAtomsMax = std::max(nAtomsMax, sum2);
913 writelog(
"estimateEP_Brinkley::\n\n");
916 writelog(
"Initial mole numbers and mu_SS:\n");
917 writelog(
" Name MoleNum mu_SS actCoeff\n");
918 for (
size_t k = 0; k <
m_kk; k++) {
922 writelogf(
"Initial n_t = %10.5g\n", n_t);
923 writelog(
"Comparison of Goal Element Abundance with Initial Guess:\n");
924 writelog(
" eName eCurrent eGoal\n");
925 for (
size_t m = 0; m <
m_mm; m++) {
930 for (
size_t m = 0; m <
m_mm; m++) {
939 for (
size_t m = 0; m <
m_mm; m++) {
945 writelogf(
"START ITERATION %d:\n", iter);
948 double n_t_calc =
calcEmoles(s, x, n_t, Xmol_i_calc, eMolesCalc, n_i_calc,
951 for (
size_t k = 0; k <
m_kk; k++) {
952 Xmol_i_calc[k] = n_i_calc[k]/n_t_calc;
956 writelog(
" Species: Calculated_Moles Calculated_Mole_Fraction\n");
957 for (
size_t k = 0; k <
m_kk; k++) {
961 writelogf(
"%15s: %10.5g\n",
"Total Molar Sum", n_t_calc);
962 writelogf(
"(iter %d) element moles bal: Goal Calculated\n", iter);
963 for (
size_t m = 0; m <
m_mm; m++) {
969 bool normalStep =
true;
972 for (
size_t m = 0; m <
m_mm; m++) {
973 if (elMoles[m] > 0.001 * elMolesTotal) {
974 if (eMolesCalc[m] > 1000. * elMoles[m]) {
978 if (1000 * eMolesCalc[m] < elMoles[m]) {
985 writelogf(
" NOTE: iter(%d) Doing an abnormal step due to row %d\n", iter, iM);
990 for (
size_t im = 0; im <
m_mm; im++) {
991 size_t m = m_orderVectorElements[im];
993 if (im <
m_nComponents && elMoles[m] > 0.001 * elMolesTotal) {
994 if (eMolesCalc[m] > 1000. * elMoles[m]) {
998 if (1000 * eMolesCalc[m] < elMoles[m]) {
1004 if (n_t < (elMolesTotal / nAtomsMax)) {
1005 if (resid[
m_mm] < 0.0) {
1008 }
else if (n_t > elMolesTotal) {
1009 resid[
m_mm] = std::min(resid[
m_mm], 0.0);
1034 for (
size_t m = 0; m <
m_mm; m++) {
1038 double nCutoff = 1.0E-9 * n_t_calc;
1040 writelog(
" Lump Sum Elements Calculation: \n");
1042 for (
size_t m = 0; m <
m_mm; m++) {
1044 size_t kMSp2 =
npos;
1045 int nSpeciesWithElem = 0;
1046 for (
size_t k = 0; k <
m_kk; k++) {
1047 if (n_i_calc[k] > nCutoff && fabs(
nAtoms(k,m)) > 0.001) {
1051 double factor = fabs(
nAtoms(kMSp,m) /
nAtoms(kMSp2,m));
1052 for (
size_t n = 0; n <
m_mm; n++) {
1053 if (fabs(factor *
nAtoms(kMSp2,n) -
nAtoms(kMSp,n)) > 1.0E-8) {
1070 for (
size_t im = 0; im <
m_mm; im++) {
1071 size_t m = m_orderVectorElements[im];
1073 for (
size_t n = 0; n <
m_mm; n++) {
1075 for (
size_t k = 0; k <
m_kk; k++) {
1079 a1(m,
m_mm) = eMolesCalc[m];
1080 a1(
m_mm, m) = eMolesCalc[m];
1082 for (
size_t n = 0; n <=
m_mm; n++) {
1093 for (
size_t im = 0; im <
m_mm; im++) {
1094 size_t m = m_orderVectorElements[im];
1096 resid[m] = elMoles[m] - eMolesCalc[m];
1115 for (
size_t m = 0; m <
m_mm; m++) {
1116 if (a1(m,m) < 1.0E-50) {
1118 writelogf(
" NOTE: Diagonalizing the analytical Jac row %d\n", m);
1120 for (
size_t n = 0; n <
m_mm; n++) {
1124 if (resid[m] > 0.0) {
1126 }
else if (resid[m] < 0.0) {
1134 resid[
m_mm] = n_t - n_t_calc;
1138 for (
size_t m = 0; m <=
m_mm; m++) {
1140 for (
size_t n = 0; n <=
m_mm; n++) {
1147 sum += pow(resid[
m_mm] /(n_t + 1.0E-15), 2);
1149 writelogf(
"(it %d) Convergence = %g\n", iter, sum);
1162 for (
size_t m = 0; m <=
m_mm; m++) {
1164 for (
size_t n = 0; n <=
m_mm; n++) {
1165 tmp += fabs(a1(m,n));
1167 if (m <
m_mm && tmp < 1.0E-30) {
1169 writelogf(
" NOTE: Diagonalizing row %d\n", m);
1171 for (
size_t n = 0; n <=
m_mm; n++) {
1179 for (
size_t n = 0; n <=
m_mm; n++) {
1187 for (
size_t m = 0; m <=
m_mm; m++) {
1189 for (
size_t n = 0; n <=
m_mm; n++) {
1215 modifiedMatrix =
false;
1216 for (
size_t m = 0; m <
m_mm; m++) {
1217 size_t sameAsRow =
npos;
1218 for (
size_t im = 0; im < m; im++) {
1219 bool theSame =
true;
1220 for (
size_t n = 0; n <
m_mm; n++) {
1221 if (fabs(a1(m,n) - a1(im,n)) > 1.0E-7) {
1230 if (sameAsRow !=
npos || lumpSum[m]) {
1233 writelogf(
"Lump summing row %d, due to rank deficiency analysis\n", m);
1234 }
else if (sameAsRow !=
npos) {
1235 writelogf(
"Identified that rows %d and %d are the same\n", m, sameAsRow);
1238 modifiedMatrix =
true;
1239 for (
size_t n = 0; n <
m_mm; n++) {
1241 a1(m,m) += fabs(a1(m,n));
1249 writelog(
"Row Summed, MODIFIED Matrix:\n");
1250 for (
size_t m = 0; m <=
m_mm; m++) {
1252 for (
size_t n = 0; n <=
m_mm; n++) {
1263 throw CanteraError(
"equilibrate:estimateEP_Brinkley()",
1264 "Jacobian is singular. \nTry adding more species, " 1265 "changing the elemental composition slightly, \nor removing " 1272 for (
size_t m = 0; m <
m_mm; m++) {
1273 if (resid[m] > 1.0) {
1274 beta = std::min(beta, 1.0 / resid[m]);
1276 if (resid[m] < -1.0) {
1277 beta = std::min(beta, -1.0 / resid[m]);
1281 writelogf(
"(it %d) Beta = %g\n", iter, beta);
1285 for (
size_t m = 0; m <
m_mm; m++) {
1286 x[m] += beta * resid[m];
1288 n_t *= exp(beta * resid[
m_mm]);
1291 writelogf(
"(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter);
1292 for (
size_t m = 0; m <
m_mm; m++) {
1293 writelogf(
" %5s %10.5g %10.5g %10.5g\n",
1304 writelogf(
" ChemEquil::estimateEP_Brinkley() SUCCESS: equilibrium found at T = %g, Pres = %g\n",
1307 writelogf(
" ChemEquil::estimateEP_Brinkley() FAILURE: equilibrium not found at T = %g, Pres = %g\n",
1320 if (fabs(elMolesGoal[
m_eloc]) > 1.0E-20) {
1324 size_t maxPosEloc =
npos;
1325 size_t maxNegEloc =
npos;
1326 double maxPosVal = -1.0;
1327 double maxNegVal = -1.0;
1329 for (
size_t k = 0; k <
m_kk; k++) {
1341 double sumPos = 0.0;
1342 double sumNeg = 0.0;
1343 for (
size_t k = 0; k <
m_kk; k++) {
1353 if (sumPos >= sumNeg) {
1354 if (sumPos <= 0.0) {
1357 double factor = (elMolesGoal[
m_eloc] + sumNeg) / sumPos;
1358 if (
m_loglevel > 0 && factor < 0.9999999999) {
1359 writelogf(
"adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n",
1363 for (
size_t k = 0; k <
m_kk; k++) {
1369 double factor = (-elMolesGoal[
m_eloc] + sumPos) / sumNeg;
1370 if (
m_loglevel > 0 && factor < 0.9999999999) {
1371 writelogf(
"adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n",
1375 for (
size_t k = 0; k <
m_kk; k++) {
EquilOpt options
Options controlling how the calculation is carried out.
size_t nElements() const
Number of elements.
void restoreState(const vector_fp &state)
Restore a state saved on a previous call to saveState.
vector_fp m_muSS_RT
Dimensionless values of the Gibbs free energy for the standard state of each species, at the temperature and pressure of the solution (the star standard state).
size_t BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, MultiPhase *mphase, std::vector< size_t > &orderVectorSpecies, std::vector< size_t > &orderVectorElements, vector_fp &formRxnMatrix)
Choose the optimum basis of species for the equilibrium calculations.
thermo_t * m_phase
Pointer to the ThermoPhase object used to initialize this object.
doublereal relTolerance
Relative tolerance.
const doublereal OneAtm
One atmosphere [Pa].
doublereal temperature() const
Temperature (K).
int setInitialMoles(thermo_t &s, vector_fp &elMoleGoal, int loglevel=0)
Estimate the initial mole numbers.
doublereal nAtoms(size_t k, size_t m) const
number of atoms of element m in species k.
void saveState(vector_fp &state) const
Save the current internal state of the phase.
doublereal moleFraction(size_t k) const
Return the mole fraction of a single species.
const size_t npos
index returned by functions to indicate "no position"
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
virtual void getGibbs_RT(doublereal *grt) const
Get the nondimensional Gibbs functions for the species in their standard states at the current T and ...
size_t m_kk
number of species in the phase
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
void addPhase(ThermoPhase *p, doublereal moles)
Add a phase to the mixture.
doublereal m_elementTotalSum
Current value of the sum of the element abundances given the current element potentials.
T clip(const T &value, const T &lower, const T &upper)
Clip value such that lower <= value <= upper.
size_t nSpecies() const
Returns the number of species in the phase.
virtual doublereal density() const
Density (kg/m^3).
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
void initialize(thermo_t &s)
virtual doublereal minTemp(size_t k=npos) const
Minimum temperature for which the thermodynamic data for the species or phase are valid...
doublereal enthalpy_mass() const
Specific enthalpy. Units: J/kg.
size_t m_nComponents
This is equal to the rank of the stoichiometric coefficient matrix when it is computed.
vector_fp m_elementmolefracs
Current value of the element mole fractions.
void ElemRearrange(size_t nComponents, const vector_fp &elementAbundances, MultiPhase *mphase, std::vector< size_t > &orderVectorSpecies, std::vector< size_t > &orderVectorElements)
Handles the potential rearrangement of the constraint equations represented by the Formula Matrix...
void equilResidual(thermo_t &s, const vector_fp &x, const vector_fp &elmtotal, vector_fp &resid, double xval, double yval, int loglevel=0)
Evaluates the residual vector F, of length m_mm.
void update(const thermo_t &s)
Update internally stored state information.
int iterations
Iteration counter.
doublereal RT() const
Return the Gas Constant multiplied by the current temperature.
Base class for a phase with thermodynamic properties.
virtual std::string getMessage() const
Method overridden by derived classes to format the error message.
std::vector< int > vector_int
Vector of ints.
doublereal atomicWeight(size_t m) const
Atomic weight of element m.
void setToEquilState(thermo_t &s, const vector_fp &x, doublereal t)
A class for multiphase mixtures.
size_t m_eloc
Index of the element id corresponding to the electric charge of each species.
std::string speciesName(size_t k) const
Name of the species with index k.
vector_fp m_comp
Storage of the element compositions. natom(k,m) = m_comp[k*m_mm+ m];.
doublereal absElemTol
Abs Tol in element number.
double m_elemFracCutoff
element fractional cutoff, below which the element will be zeroed.
doublereal entropy_mass() const
Specific entropy. Units: J/kg/K.
int estimateEP_Brinkley(thermo_t &s, vector_fp &lambda, vector_fp &elMoles)
Base class for exceptions thrown by Cantera classes.
virtual void setMoleFractions(const doublereal *const x)
Set the mole fractions to the specified values.
vector_fp & data()
Return a reference to the data vector.
int estimateElementPotentials(thermo_t &s, vector_fp &lambda, vector_fp &elMolesGoal, int loglevel=0)
Generate a starting estimate for the element potentials.
size_t m_mm
number of elements in the phase
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
virtual void getActivityCoefficients(doublereal *ac) const
Get the array of non-dimensional molar-based activity coefficients at the current solution temperatur...
doublereal dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
int equilibrate(thermo_t &s, const char *XY, int loglevel=0)
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
int m_loglevel
Verbosity of printed output.
virtual void setPressure(doublereal p)
Set the internally stored pressure (Pa) at constant temperature and composition.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
virtual void setTemperature(const doublereal temp)
Set the internally stored temperature of the phase (K).
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
vector_fp m_lambda
Current value of the dimensional element potentials.
virtual doublereal refPressure() const
Returns the reference pressure in Pa.
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
void init()
Process phases and build atomic composition array.
Contains declarations for string manipulation functions within Cantera.
vector_fp m_molefractions
Current value of the mole fractions in the single phase. length = m_kk.
double calcEmoles(thermo_t &s, vector_fp &x, const double &n_t, const vector_fp &Xmol_i_calc, vector_fp &eMolesCalc, vector_fp &n_i_calc, double pressureConst)
Given a vector of dimensionless element abundances, this routine calculates the moles of the elements...
int maxIterations
Maximum number of iterations.
int _equilflag(const char *xy)
map property strings to integers
virtual void setToEquilState(const doublereal *lambda_RT)
This method is used by the ChemEquil equilibrium solver.
Namespace for the Cantera kernel.
doublereal intEnergy_mass() const
Specific internal energy. Units: J/kg.
virtual doublereal maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
doublereal nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
std::string elementName(size_t m) const
Name of the element with index m.
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
int dampStep(thermo_t &s, vector_fp &oldx, double oldf, vector_fp &grad, vector_fp &step, vector_fp &x, double &f, vector_fp &elmols, double xval, double yval)
Find an acceptable step size and take it.