26void printProgress(
const vector<string>& spName,
const vector<double>& soln,
27 const vector<double>& ff)
30 plogf(
" --- Summary of current progress:\n");
31 plogf(
" --- Name Moles - SSGibbs \n");
32 plogf(
" -------------------------------------------------------------------------------------\n");
33 for (
size_t k = 0; k < soln.size(); k++) {
34 plogf(
" --- %20s %12.4g - %12.4g\n", spName[k], soln[k], ff[k]);
35 sum += soln[k] * ff[k];
37 plogf(
" --- Total sum to be minimized = %g\n", sum);
47 m_nsp(mphase->nSpecies()),
48 m_numSpeciesRdc(mphase->nSpecies()),
49 m_numPhases(mphase->nPhases()),
50 m_temperature(mphase->temperature()),
51 m_pressurePA(mphase->pressure()),
52 m_totalVol(mphase->volume()),
56 for (
size_t kspec = 0; kspec <
m_nsp; kspec++) {
60 string ser =
"VCS_SOLVE: ERROR:\n\t";
62 plogf(
"%s Number of species is nonpositive\n", ser);
64 " Number of species is nonpositive\n");
67 plogf(
"%s Number of phases is nonpositive\n", ser);
69 " Number of phases is nonpositive\n");
150 for (
size_t iphase = 0; iphase <
m_numPhases; iphase++) {
156 string eos = tPhase->
type();
157 bool gasPhase = (eos ==
"ideal-gas");
160 size_t nSpPhase = tPhase->
nSpecies();
162 string phaseName = tPhase->
name();
172 VolPhase->
resize(iphase, nSpPhase, nelem, phaseName.c_str(), 0.0);
188 if (eos ==
"ideal-gas") {
190 }
else if (eos ==
"fixed-stoichiometry") {
192 }
else if (eos ==
"ideal-condensed") {
194 }
else if (tPhase->
nDim() != 3) {
196 "Surface/edge phase not handled yet.");
199 writelog(
"Unknown Cantera EOS to VCSnonideal: '{}'\n", eos);
218 for (
size_t k = 0; k < nSpPhase; k++) {
258 for (
size_t e = 0; e <
m_nelem; e++) {
279 double minTemp, maxTemp, refPressure;
280 sp.
reportParams(k, spType, c, minTemp, maxTemp, refPressure);
295 plogf(
"vcs_Cantera_convert: Species Type %d not known \n",
317 for (
size_t k = 0; k < nSpPhase; k++) {
326 for (
size_t j = 0; j <
m_nelem; j++) {
327 for (
size_t kspec = 0; kspec <
m_nsp; kspec++) {
339 writeline(
'=', 80,
true,
true);
340 writeline(
'=', 16,
false);
341 plogf(
" Cantera_to_vprob: START OF PROBLEM STATEMENT ");
344 plogf(
" Phase IDs of species\n");
345 plogf(
" species phaseID phaseName ");
346 plogf(
" Initial_Estimated_kMols\n");
347 for (
size_t i = 0; i <
m_nsp; i++) {
361 writeline(
'-', 80,
true,
true);
362 plogf(
" Information about phases\n");
363 plogf(
" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
364 plogf(
" TMolesInert Tmoles(kmol)\n");
366 for (
size_t iphase = 0; iphase <
m_numPhases; iphase++) {
368 plogf(
"%16s %5d %5d %8d %16s %8d %16e ", VolPhase->
PhaseName.c_str(),
375 writeline(
'=', 80,
true,
true);
376 writeline(
'=', 16,
false);
377 plogf(
" Cantera_to_vprob: END OF PROBLEM STATEMENT ");
391 for (
size_t i = 0; i <
m_nsp; i++) {
396 for (
size_t i = 0; i <
m_nelem; i++) {
403 for (
size_t kspec = 0; kspec <
m_nsp; kspec++) {
407 "Species to Phase Mapping, PhaseID, has a bad value\n"
408 "\tm_phaseID[{}] = {}\n"
409 "Allowed values: 0 to {}", kspec, iph,
m_numPhases - 1);
417 if (numPhSp[iph] != Vphase->
nSpecies()) {
419 "Number of species in phase {}, {}, doesn't match ({} != {}) [vphase = {}]",
424 for (
size_t i = 0; i <
m_nelem; i++) {
429 "Charge neutrality condition {} is signicantly "
430 "nonzero, {}. Giving up",
434 plogf(
"Charge neutrality condition %s not zero, %g. Setting it zero\n",
444 for (
size_t i = 0; i <
m_nsp; i++) {
461 for (
size_t k = 1; k < Vphase->
nSpecies(); k++) {
471VCS_SOLVE::~VCS_SOLVE()
491 plogf(
"vcs_prep_oneTime returned a bad status, %d: bailing!\n",
515 if (ipr > 0 || ip1 > 0) {
522 }
else if (iconv == 1) {
523 plogf(
"WARNING: RANGE SPACE ERROR encountered\n");
532 "called for a phase that exists!");
538 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
541 "VCS_SOLVE::vcs_popPhasePossible",
542 "we shouldn't be here {}: {} > 0.0", kspec,
546 bool iPopPossible =
true;
556 double negChangeComp = - stoicC;
557 if (negChangeComp > 0.0) {
561 iPopPossible =
false;
576 for (
size_t jrxn = 0; jrxn <
m_numRxnRdc; jrxn++) {
577 bool foundJrxn =
false;
619 size_t iphasePop =
npos;
620 double FephaseMax = -1.0E30;
621 double Fephase = -1.0E30;
625 plogf(
" --- vcs_popPhaseID() called\n");
626 plogf(
" --- Phase Status F_e MoleNum\n");
627 plogf(
" --------------------------------------------------------------------------\n");
631 int existence = Vphase->
exists();
635 plogf(
" --- %18s %5d NA %11.3e\n",
645 "Index out of bounds due to logic error.");
648 Fephase = exp(-deltaGRxn) - 1.0;
650 strcpy(anote,
" (ready to be birthed)");
651 if (Fephase > FephaseMax) {
653 FephaseMax = Fephase;
654 strcpy(anote,
" (chosen to be birthed)");
658 strcpy(anote,
" (not stable)");
660 "VCS_SOLVE::vcs_popPhaseID",
"shouldn't be here");
664 plogf(
" --- %18s %5d %10.3g %10.3g %s\n",
673 if (Fephase > FephaseMax) {
675 FephaseMax = Fephase;
678 FephaseMax = std::max(FephaseMax, Fephase);
681 plogf(
" --- %18s %5d %11.3g %11.3g\n",
687 plogf(
" --- %18s %5d blocked %11.3g\n",
695 phasePopPhaseIDs.resize(0);
696 if (iphasePop !=
npos) {
697 phasePopPhaseIDs.push_back(iphasePop);
703 plogf(
" ---------------------------------------------------------------------\n");
715 vector<size_t> creationGlobalRxnNumbers;
723 "called for a phase that exists!");
725 plogf(
" --- vcs_popPhaseRxnStepSizes() called to pop phase %s %d into existence\n",
772 vector<double> fracDelta(Vphase->
nSpecies());
773 vector<double> X_est(Vphase->
nSpecies());
776 double sumFrac = 0.0;
777 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
778 sumFrac += fracDelta[k];
780 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
781 X_est[k] = fracDelta[k] / sumFrac;
784 double deltaMolNumPhase = tPhaseMoles;
789 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
791 double delmol = deltaMolNumPhase * X_est[k];
795 throw CanteraError(
"VCS_SOLVE::vcs_popPhaseRxnStepSizes",
796 "Index out of bounds due to logic error.");
801 molNumSpecies_tmp[j] += stoicC * delmol;
807 double ratioComp = 0.0;
810 if (molNumSpecies_tmp[j] < 0.0) {
814 damp = std::min(damp, delta0 / deltaJ * 0.9);
819 if ((jph != iphasePop) && (!
m_SSPhase[j])) {
820 double fdeltaJ = fabs(deltaJ);
833 if (ratioComp > 1.0E-30 && ratioComp < 0.001) {
834 damp = 0.001 / ratioComp;
836 if (damp <= 1.0E-6) {
840 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
846 if (X_est[k] > 1.0E-3) {
859 size_t iphDel =
npos;
864 for (
int j = 0; j < 82; j++) {
868 plogf(
" --- Subroutine vcs_RxnStepSizes called - Details:\n");
870 for (
int j = 0; j < 82; j++) {
874 plogf(
" --- Species KMoles Rxn_Adjustment DeltaG"
885 for (
size_t irxn = 0; irxn <
m_numRxnRdc; ++irxn) {
886 ANOTE =
"Normal Calc";
891 ANOTE =
"ZeroedPhase: Phase is artificially zeroed";
910 ANOTE = fmt::sprintf(
"MultSpec (%s): Species not born due to STOICH/PHASEPOP even though DG = %11.3E",
914 ANOTE = fmt::sprintf(
"MultSpec (%s): small species born again DG = %11.3E",
918 ANOTE = fmt::sprintf(
"MultSpec (%s):still dead, no phase pop, even though DG = %11.3E",
921 if (Vphase->
exists() > 0 && trphmoles > 0.0) {
923 ANOTE = fmt::sprintf(
"MultSpec (%s): birthed species because it was zero in a small existing phase with DG = %11.3E",
928 ANOTE = fmt::sprintf(
"MultSpec (%s): still dead DG = %11.3E",
939 ANOTE = fmt::sprintf(
"Skipped: superconverged DG = %11.3E",
m_deltaGRxn_new[irxn]);
942 plogf(
" %12.4E %12.4E %12.4E | %s\n",
952 ANOTE = fmt::sprintf(
"Skipped: IC = %3d and DG >0: %11.3E",
956 plogf(
" %12.4E %12.4E %12.4E | %s\n",
988 ANOTE = fmt::sprintf(
"Normal calc: diag adjusted from %g "
989 "to %g due to act coeff", s_old, s);
1000 ANOTE = fmt::sprintf(
"Delta damped from %g "
1005 ANOTE = fmt::sprintf(
"Delta damped from %g "
1016 ANOTE = fmt::sprintf(
"Delta damped from %g "
1064 if ((k == kspec) && (
m_SSPhase[kspec] != 1)) {
1071 ANOTE = fmt::sprintf(
"Delta damped from %g to %g due to delete %s",
m_deltaMolNumSpecies[kspec],
1076 plogf(
" %12.4E %12.4E %12.4E | %s\n",
1084 for (
size_t j = 0; j <
m_nsp; j++) {
1096 ANOTE = fmt::sprintf(
"Delete component SS phase %d named %s - SS phases only",
1099 ANOTE = fmt::sprintf(
"Delete this SS phase %d - SS components only", iphDel);
1103 plogf(
" %12.4E %12.4E %12.4E | %s\n",
1106 plogf(
" --- vcs_RxnStepSizes Special section to set up to delete %s\n",
1110 forceComponentCalc = 1;
1123 plogf(
" %12.4E %12.4E %12.4E | %s\n",
1140 const size_t nsp = Vphase->
nSpecies();
1141 int minNumberIterations = 3;
1143 minNumberIterations = 1;
1147 bool doSuccessiveSubstitution =
true;
1148 double funcPhaseStability;
1149 vector<double> X_est(nsp, 0.0);
1150 vector<double> delFrac(nsp, 0.0);
1151 vector<double> E_phi(nsp, 0.0);
1152 vector<double> fracDelta_old(nsp, 0.0);
1153 vector<double> fracDelta_raw(nsp, 0.0);
1154 vector<size_t> creationGlobalRxnNumbers(nsp,
npos);
1165 vector<size_t> componentList;
1166 for (
size_t k = 0; k < nsp; k++) {
1169 componentList.push_back(k);
1173 double normUpdate = 0.1 *
vcs_l2norm(fracDelta_new);
1174 double damp = 1.0E-2;
1176 if (doSuccessiveSubstitution) {
1179 plogf(
" --- vcs_phaseStabilityTest() called\n");
1180 plogf(
" --- Its X_old[%2d] FracDel_old[%2d] deltaF[%2d] FracDel_new[%2d]"
1181 " normUpdate damp FuncPhaseStability\n", KP, KP, KP, KP);
1182 plogf(
" --------------------------------------------------------------"
1183 "--------------------------------------------------------\n");
1185 plogf(
" --- vcs_phaseStabilityTest() called for phase %d\n", iph);
1188 for (
size_t k = 0; k < nsp; k++) {
1189 if (fracDelta_new[k] < 1.0E-13) {
1190 fracDelta_new[k] =1.0E-13;
1193 bool converged =
false;
1194 double dirProd = 0.0;
1195 for (
int its = 0; its < 200 && (!converged); its++) {
1196 double dampOld = damp;
1197 double normUpdateOld = normUpdate;
1198 fracDelta_old = fracDelta_new;
1199 double dirProdOld = dirProd;
1203 for (
size_t i = 0; i < componentList.size(); i++) {
1204 size_t kc = componentList[i];
1206 fracDelta_old[kc] = 0.0;
1207 for (
size_t k = 0; k < nsp; k++) {
1217 double sumFrac = 0.0;
1218 for (
size_t k = 0; k < nsp; k++) {
1219 sumFrac += fracDelta_old[k];
1223 if (sumFrac <= 0.0) {
1226 double sum_Xcomp = 0.0;
1227 for (
size_t k = 0; k < nsp; k++) {
1228 X_est[k] = fracDelta_old[k] / sumFrac;
1230 sum_Xcomp += X_est[k];
1243 for (
size_t i = 0; i < componentList.size(); i++) {
1244 size_t kc = componentList[i];
1255 for (
size_t i = 0; i < componentList.size(); i++) {
1257 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
1274 funcPhaseStability = sum_Xcomp - 1.0;
1275 for (
size_t k = 0; k < nsp; k++) {
1282 funcPhaseStability += E_phi[k];
1289 for (
size_t k = 0; k < nsp; k++) {
1291 double b = E_phi[k] / sum * (1.0 - sum_Xcomp);
1293 fracDelta_raw[k] = b;
1299 for (
size_t i = 0; i < componentList.size(); i++) {
1300 size_t kc = componentList[i];
1302 fracDelta_raw[kc] = 0.0;
1303 for (
size_t k = 0; k < nsp; k++) {
1313 for (
size_t k = 0; k < nsp; k++) {
1314 delFrac[k] = fracDelta_raw[k] - fracDelta_old[k];
1319 for (
size_t k = 0; k < nsp; k++) {
1320 dirProd += fracDelta_old[k] * delFrac[k];
1322 bool crossedSign =
false;
1323 if (dirProd * dirProdOld < 0.0) {
1328 if (dampOld < 0.25) {
1329 damp = 2.0 * dampOld;
1332 if (normUpdate *1.5 > normUpdateOld) {
1333 damp = 0.5 * dampOld;
1334 }
else if (normUpdate *2.0 > normUpdateOld) {
1335 damp = 0.8 * dampOld;
1338 if (normUpdate > normUpdateOld * 2.0) {
1339 damp = 0.6 * dampOld;
1340 }
else if (normUpdate > normUpdateOld * 1.2) {
1341 damp = 0.9 * dampOld;
1345 for (
size_t k = 0; k < nsp; k++) {
1346 if (fabs(damp * delFrac[k]) > 0.3*fabs(fracDelta_old[k])) {
1347 damp = std::max(0.3*fabs(fracDelta_old[k]) / fabs(delFrac[k]),
1348 1.0E-8/fabs(delFrac[k]));
1350 if (delFrac[k] < 0.0 && 2.0 * damp * (-delFrac[k]) > fracDelta_old[k]) {
1351 damp = fracDelta_old[k] / (2.0 * -delFrac[k]);
1353 if (delFrac[k] > 0.0 && 2.0 * damp * delFrac[k] > fracDelta_old[k]) {
1354 damp = fracDelta_old[k] / (2.0 * delFrac[k]);
1357 damp = std::max(damp, 0.000001);
1358 for (
size_t k = 0; k < nsp; k++) {
1359 fracDelta_new[k] = fracDelta_old[k] + damp * delFrac[k];
1363 plogf(
" --- %3d %12g %12g %12g %12g %12g %12g %12g\n", its, X_est[KP], fracDelta_old[KP],
1364 delFrac[KP], fracDelta_new[KP], normUpdate, damp, funcPhaseStability);
1367 if (normUpdate < 1.0E-5 * damp) {
1369 if (its < minNumberIterations) {
1384 throw CanteraError(
"VCS_SOLVE::vcs_phaseStabilityTest",
"not done yet");
1387 plogf(
" ------------------------------------------------------------"
1388 "-------------------------------------------------------------\n");
1390 if (funcPhaseStability > 0.0) {
1391 plogf(
" --- phase %d with func = %g is to be born\n", iph, funcPhaseStability);
1393 plogf(
" --- phase %d with func = %g stays dead\n", iph, funcPhaseStability);
1396 return funcPhaseStability;
1418 plogf(
"vcs_inest_TP returned a failure flag\n");
1438 for (
size_t k = 0; k <
m_nsp; k++) {
1447 for (
size_t i = 0; i <
m_nsp; ++i) {
1459 const double w[],
double volPM[])
1461 double VolTot = 0.0;
1462 for (
size_t iphase = 0; iphase <
m_numPhases; iphase++) {
1494 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
1501 for (
size_t e = 0; e < eSize; e++) {
1525 double test = -1.0e-10;
1526 bool modifiedSoln =
false;
1529 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
1534 if (fabs(sum) < 1.0E-6) {
1535 modifiedSoln =
true;
1538 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
1552 double* aw = &awSpace[0];
1554 plogf(
"vcs_prep_oneTime: failed to get memory: global bailout\n");
1555 return VCS_NOMEMORY;
1557 double* sa = aw +
m_nsp;
1561 retn =
vcs_basopt(
true, aw, sa, sm, ss, test, &conv);
1563 plogf(
"vcs_prep_oneTime:");
1564 plogf(
" Determination of number of components failed: %d\n",
1566 plogf(
" Global Bailout!\n");
1587 plogf(
"vcs_prep_oneTime:");
1588 plogf(
" Determination of element reordering failed: %d\n",
1590 plogf(
" Global Bailout!\n");
1597 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
1616 for (
size_t e = 0; e < m_mix->
nElements(); e++) {
1619 if (sum < 1.0E-20) {
1621 plogf(
"vcs has determined the problem is not well posed: Bailing\n");
1628 double*
const sm,
double*
const ss)
1634 plogf(
" --- Subroutine elem_rearrange() called to ");
1635 plogf(
"check stoich. coefficient matrix\n");
1636 plogf(
" --- and to rearrange the element ordering once\n");
1642 double test = -1.0E10;
1645 for (
size_t i = 0; i <
m_nelem; ++i) {
1648 if (test == aw[i]) {
1657 while (jr < ncomponents) {
1666 for (
size_t ielem = jr; ielem <
m_nelem; ielem++) {
1674 "Shouldn't be here. Algorithm misfired.");
1691 for (
size_t j = 0; j < ncomponents; ++j) {
1698 for (
size_t j = 0; j < jl; ++j) {
1700 for (
size_t i = 0; i < ncomponents; ++i) {
1701 ss[j] += sm[i + jr*ncomponents] * sm[i + j*ncomponents];
1708 for (
size_t j = 0; j < jl; ++j) {
1709 for (
size_t i = 0; i < ncomponents; ++i) {
1710 sm[i + jr*ncomponents] -= ss[j] * sm[i + j*ncomponents];
1718 for (
size_t ml = 0; ml < ncomponents; ++ml) {
1719 sa[jr] += pow(sm[ml + jr*ncomponents], 2);
1722 if (sa[jr] > 1.0e-6) {
1729 plogf(
" --- %-2.2s(%9.2g) replaces %-2.2s(%9.2g) as element %3d\n",
1734 std::swap(aw[jr], aw[k]);
1749 "vcs_switch_elem_pos",
1750 "inappropriate args: {} {}", ipos, jpos);
1770 for (
size_t j = 0; j <
m_nsp; ++j) {
1778 double diag = hessianDiag_Ideal;
1780 if (hessianDiag_Ideal <= 0.0) {
1782 "We shouldn't be here");
1784 if (hessActCoef >= 0.0) {
1785 diag += hessActCoef;
1786 }
else if (fabs(hessActCoef) < 0.6666 * hessianDiag_Ideal) {
1787 diag += hessActCoef;
1789 diag -= 0.6666 * hessianDiag_Ideal;
1827 for (
size_t iphase = 0; iphase <
m_numPhases; iphase++) {
1844 bool inertYes =
false;
1847 vector<pair<double, size_t>> x_order;
1848 for (
size_t i = 0; i <
m_nsp; i++) {
1861 plogf(
"\t\t VCS_TP REPORT\n");
1865 plogf(
" ERROR: CONVERGENCE CRITERION NOT SATISFIED.\n");
1866 }
else if (iconv == 1) {
1867 plogf(
" RANGE SPACE ERROR: Equilibrium Found but not all Element Abundances are Satisfied\n");
1882 plogf(
" Species Equilibrium kmoles ");
1883 plogf(
"Mole Fraction ChemPot/RT SpecUnkType\n");
1887 writeline(
' ', 13,
false);
1894 size_t j = x_order[i].second;
1896 writeline(
' ', 13,
false);
1906 throw CanteraError(
"VCS_SOLVE::vcs_report",
"we have a problem");
1914 plogf(
" Inert Gas Species ");
1916 plogf(
" Inert Species in phase %16s ",
1924 plogf(
"\n SPECIES WITH LESS THAN 1.0E-32 KMOLES:\n\n");
1928 plogf(
" %14.7E %14.7E %12.4E",
1946 plogf(
" |ComponentID|");
1951 plogf(
" | Components|");
1956 plogf(
" NonComponent | Moles |");
1960 plogf(
" | DG/RT Rxn |\n");
1962 for (
size_t irxn = 0; irxn <
m_numRxnTot; irxn++) {
1964 plogf(
" %3d ", kspec);
1977 vector<double> gaPhase(
m_nelem, 0.0);
1978 vector<double> gaTPhase(
m_nelem, 0.0);
1979 double totalMoles = 0.0;
1980 double gibbsPhase = 0.0;
1981 double gibbsTotal = 0.0;
1984 writeline(
'-',
m_nelem*10 + 58);
1985 plogf(
" | ElementID |");
1986 for (
size_t j = 0; j <
m_nelem; j++) {
1990 plogf(
" | Element |");
1991 for (
size_t j = 0; j <
m_nelem; j++) {
1995 plogf(
" PhaseName |KMolTarget |");
1996 for (
size_t j = 0; j <
m_nelem; j++) {
1999 plogf(
" | Gibbs Total |\n");
2000 writeline(
'-',
m_nelem*10 + 58);
2001 for (
size_t iphase = 0; iphase <
m_numPhases; iphase++) {
2002 plogf(
" %3d ", iphase);
2009 throw CanteraError(
"VCS_SOLVE::vcs_report",
"we have a problem");
2012 for (
size_t j = 0; j <
m_nelem; j++) {
2013 plogf(
" %10.3g", gaPhase[j]);
2014 gaTPhase[j] += gaPhase[j];
2018 gibbsTotal += gibbsPhase;
2019 plogf(
" | %18.11E |\n", gibbsPhase);
2021 writeline(
'-',
m_nelem*10 + 58);
2022 plogf(
" TOTAL |%10.3e |", totalMoles);
2023 for (
size_t j = 0; j <
m_nelem; j++) {
2024 plogf(
" %10.3g", gaTPhase[j]);
2026 plogf(
" | %18.11E |\n", gibbsTotal);
2028 writeline(
'-',
m_nelem*10 + 58);
2037 plogf(
"\n\tTotal Dimensionless Gibbs Free Energy = G/RT = %15.7E\n", g);
2039 plogf(
"\t\t(Inert species have standard free energy of zero)\n");
2042 plogf(
"\nElemental Abundances (kmol): ");
2043 plogf(
" Actual Target Type ElActive\n");
2044 for (
size_t i = 0; i <
m_nelem; ++i) {
2045 writeline(
' ', 26,
false);
2053 writeline(
'-', 93,
true,
true);
2054 plogf(
"Chemical Potentials of the Species: (dimensionless)\n");
2057 plogf(
" Name TKMoles StandStateChemPot "
2058 " ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)");
2059 plogf(
"| (MolNum ChemPot)|");
2060 writeline(
'-', 147,
true,
true);
2061 for (
size_t i = 0; i <
m_nsp; ++i) {
2062 size_t j = x_order[i].second;
2077 lx = log(tmp) - log(tpmoles);
2083 plogf(
"%14.7E |", lx);
2084 plogf(
"%14.7E | ", eContrib);
2089 "we have a problem - doesn't add up");
2101 for (
size_t i = 0; i < 125; i++) {
2104 plogf(
" %20.9E\n", g);
2105 writeline(
'-', 147);
2109 plogf(
"\nCounters: Iterations Time (seconds)\n");
2111 plogf(
" vcs_basopt: %5d %11.5E\n",
2113 plogf(
" vcs_TP: %5d %11.5E\n",
2116 plogf(
" vcs_basopt: %5d %11s\n",
2118 plogf(
" vcs_TP: %5d %11s\n",
2130 for (
size_t j = 0; j <
m_nelem; ++j) {
2132 for (
size_t i = 0; i <
m_nsp; ++i) {
2147 for (
size_t i = 0; i < top; ++i) {
2154 "Problem with charge neutrality condition");
2163 bool multisign =
false;
2164 for (
size_t kspec = 0; kspec <
m_nsp; kspec++) {
2198 for (
size_t j = 0; j <
m_nelem; ++j) {
2199 elemAbundPhase[j] = 0.0;
2200 for (
size_t i = 0; i <
m_nsp; ++i) {
2214 plogf(
" --- vcsc_elcorr: Element abundances correction routine");
2216 plogf(
" (m_numComponents != m_nelem)");
2221 for (
size_t i = 0; i <
m_nelem; ++i) {
2224 double l2before = 0.0;
2225 for (
size_t i = 0; i <
m_nelem; ++i) {
2226 l2before += x[i] * x[i];
2228 l2before = sqrt(l2before/
m_nelem);
2233 bool changed =
false;
2234 for (
size_t i = 0; i <
m_nelem; ++i) {
2236 bool multisign =
false;
2237 for (
size_t kspec = 0; kspec <
m_nsp; kspec++) {
2249 if (numNonZero < 2) {
2250 for (
size_t kspec = 0; kspec <
m_nsp; kspec++) {
2260 int numCompNonZero = 0;
2261 size_t compID =
npos;
2271 if (numCompNonZero == 1) {
2297 for (
size_t i = 0; i <
m_nelem; ++i) {
2300 for (
size_t kspec = 0; kspec <
m_nsp; kspec++) {
2303 if (atomComp > 0.0) {
2307 plogf(
" --- vcs_elcorr: Reduced species %s from %g to %g "
2308 "due to %s max bounds constraint\n",
2322 plogf(
" --- vcs_elcorr: Zeroed species %s and changed "
2323 "status to %d due to max bounds constraint\n",
2344 if (fabs(x[i]) > 1.0E-13) {
2361 par = std::min(par, 100.0);
2363 if (par < 1.0 && par > 0.0) {
2402 for (
size_t kspec = 0; kspec <
m_nsp; kspec++) {
2406 double saveDir = 0.0;
2407 bool goodSpec =
true;
2410 if (fabs(dir) > 1.0E-10) {
2412 if (saveDir < 0.0) {
2417 if (saveDir > 0.0) {
2462 for (
size_t i = 0; i <
m_nelem; ++i) {
2488 for (
size_t i = 0; i <
m_nelem; ++i) {
2491 bool useZeroed =
true;
2531 double l2after = 0.0;
2532 for (
size_t i = 0; i <
m_nelem; ++i) {
2535 l2after = sqrt(l2after/
m_nelem);
2537 plogf(
" --- Elem_Abund: Correct Initial "
2539 for (
size_t i = 0; i <
m_nelem; ++i) {
2544 plogf(
" --- Diff_Norm: %20.12E %20.12E\n",
2552 const char* pprefix =
" --- vcs_inest: ";
2560 plogf(
"%s Initial guess passed element abundances on input\n", pprefix);
2561 plogf(
"%s m_doEstimateEquil = 1 so will use the input mole "
2562 "numbers as estimates\n", pprefix);
2566 plogf(
"%s Initial guess failed element abundances on input\n", pprefix);
2567 plogf(
"%s m_doEstimateEquil = 1 so will discard input "
2568 "mole numbers and find our own estimate\n", pprefix);
2574 vector<double> ss(
m_nelem, 0.0);
2575 vector<double> sa(
m_nelem, 0.0);
2580 plogf(
"%sGo find an initial estimate for the equilibrium problem\n",
2583 double test = -1.0E20;
2584 vcs_inest(&aw[0], &sa[0], &sm[0], &ss[0], test);
2599 plogf(
"%sInitial guess failed element abundances\n", pprefix);
2600 plogf(
"%sCall vcs_elcorr to attempt fix\n", pprefix);
2605 plogf(
"%sInitial guess still fails element abundance equations\n",
2607 plogf(
"%s - Inability to ever satisfy element abundance "
2608 "constraints is probable\n", pprefix);
2613 plogf(
"%sInitial guess now satisfies element abundances\n", pprefix);
2615 plogf(
"%sElement Abundances RANGE ERROR\n", pprefix);
2616 plogf(
"%s - Initial guess satisfies NC=%d element abundances, "
2617 "BUT not NE=%d element abundances\n", pprefix,
2625 plogf(
"%sInitial guess satisfies element abundances\n", pprefix);
2627 plogf(
"%sElement Abundances RANGE ERROR\n", pprefix);
2628 plogf(
"%s - Initial guess satisfies NC=%d element abundances, "
2629 "BUT not NE=%d element abundances\n", pprefix,
2636 plogf(
"%sTotal Dimensionless Gibbs Free Energy = %15.7E\n", pprefix,
2649 double test = -1.0E-10;
2652 plogf(
" --- call setInitialMoles\n");
2655 double dxi_min = 1.0e10;
2658 bool abundancesOK =
true;
2659 bool usedZeroedSpecies;
2661 vector<double> ss(
m_nelem, 0.0);
2662 vector<double> sa(
m_nelem, 0.0);
2663 vector<double> wx(
m_nelem, 0.0);
2664 vector<double> aw(
m_nsp, 0.0);
2666 for (
size_t ik = 0; ik <
m_nsp; ik++) {
2680 plogf(
" --- seMolesLinProg Mole numbers failing element abundances\n");
2681 plogf(
" --- seMolesLinProg Call vcs_elcorr to attempt fix\n");
2685 abundancesOK =
false;
2687 abundancesOK =
true;
2690 abundancesOK =
true;
2696 retn =
vcs_basopt(
false, &aw[0], &sa[0], &sm[0], &ss[0],
2697 test, &usedZeroedSpecies);
2703 plogf(
"iteration %d\n", iter);
2712 for (
size_t irxn = 0; irxn <
m_numRxnTot; irxn++) {
2718 for (
size_t jcomp = 0; jcomp <
m_nelem; jcomp++) {
2724 int idir = (dg_rt < 0.0 ? 1 : -1);
2730 double nu = sc_irxn[jcomp];
2743 dxi_min = std::min(dxi_min, delta_xi);
2750 double dsLocal = idir*dxi_min;
2773 plogf(
" --- setInitialMoles end\n");
2776 if (!abundancesOK) {
2778 }
else if (iter > 15) {
2802 g += molesSp[kspec] * chemPot[kspec];
2810 const double*
const fe)
2813 double phaseMols = 0.0;
2816 g += w[kspec] * fe[kspec];
2817 phaseMols += w[kspec];
2855 vector<size_t> invSpecies(
m_nsp);
2856 for (
size_t k = 0; k <
m_nsp; k++) {
2860 for (
size_t iphase = 0; iphase <
m_numPhases; iphase++) {
2867 size_t nSpPhase = tPhase->
nSpecies();
2868 if ((nSpPhase == 1) && (volPhase->
phiVarIndex() == 0)) {
2879 writeline(
'=', 80,
true,
true);
2880 writeline(
'=', 20,
false);
2881 plogf(
" Cantera_to_vprob: START OF PROBLEM STATEMENT ");
2885 plogf(
" Phase IDs of species\n");
2886 plogf(
" species phaseID phaseName ");
2887 plogf(
" Initial_Estimated_kMols\n");
2888 for (
size_t i = 0; i <
m_nsp; i++) {
2901 writeline(
'-', 80,
true,
true);
2902 plogf(
" Information about phases\n");
2903 plogf(
" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
2904 plogf(
" TMolesInert Tmoles(kmol)\n");
2906 for (
size_t iphase = 0; iphase <
m_numPhases; iphase++) {
2908 plogf(
"%16s %5d %5d %8d %16s %8d %16e ", VolPhase->
PhaseName.c_str(),
2915 writeline(
'=', 80,
true,
true);
2916 writeline(
'=', 20,
false);
2917 plogf(
" Cantera_to_vprob: END OF PROBLEM STATEMENT ");
2936 double*
const ss,
double test)
2938 const char* pprefix =
" --- vcs_inest: ";
2946 plogf(
"%s Mole Numbers returned from linear programming (vcs_inest initial guess):\n",
2948 plogf(
"%s SPECIES MOLE_NUMBER -SS_ChemPotential\n", pprefix);
2949 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
2950 plogf(
"%s ", pprefix);
2954 plogf(
"%s Element Abundance Agreement returned from linear "
2955 "programming (vcs_inest initial guess):\n", pprefix);
2956 plogf(
"%s Element Goal Actual\n", pprefix);
2957 for (
size_t j = 0; j <
m_nelem; j++) {
2960 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
2963 plogf(
"%s ", pprefix);
2974 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
2990 vcs_basopt(
false, aw, sa, sm, ss, test, &conv);
3006 double TMolesMultiphase = 0.0;
3032 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
3033 plogf(
"%s", pprefix);
3039 plogf(
"fe* = %15.5g ff = %15.5g dg* = %15.5g\n",
3053 for (
size_t irxn = 0; irxn < nrxn; ++irxn) {
3087 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
3089 plogf(
"%sdirection (", pprefix);
3096 plogf(
" (ssPhase doesn't exist -> stability not checked)");
3113 if (par <= 1.0 && par > 0.0) {
3149 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
3155 if (s < 0.0 && ikl == 0) {
3168 double xl = (1.0 - s / (s1 - s)) * 0.5;
3178 par = par * 2.0 * xl;
3185 plogf(
"%s Final Mole Numbers produced by inest:\n",
3187 plogf(
"%s SPECIES MOLE_NUMBER\n", pprefix);
3188 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
3189 plogf(
"%s %-12.12s %g\n",
3198 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
3219 for (
size_t kspec = 0; kspec <
m_nsp; kspec++) {
3261 plogf(
"\nTCounters: Num_Calls Total_Its Total_Time (seconds)\n");
3262 if (timing_print_lvl > 0) {
3263 plogf(
" vcs_basopt: %5d %5d %11.5E\n",
3266 plogf(
" vcs_TP: %5d %5d %11.5E\n",
3269 plogf(
" vcs_inest: %5d %11.5E\n",
3271 plogf(
" vcs_TotalTime: %11.5E\n",
3274 plogf(
" vcs_basopt: %5d %5d %11s\n",
3276 plogf(
" vcs_TP: %5d %5d %11s\n",
3278 plogf(
" vcs_inest: %5d %11s\n",
3280 plogf(
" vcs_TotalTime: %11s\n",
3291 writeline(
'=', 80,
true,
true);
3292 writeline(
'=', 20,
false);
3293 plogf(
" VCS_PROB: PROBLEM STATEMENT ");
3297 plogf(
"\tSolve a constant T, P problem:\n");
3301 plogf(
"\t\tPres = %g atm\n", pres_atm);
3303 plogf(
" Phase IDs of species\n");
3304 plogf(
" species phaseID phaseName ");
3305 plogf(
" Initial_Estimated_Moles Species_Type\n");
3306 for (
size_t i = 0; i <
m_nsp; i++) {
3326 writeline(
'-', 80,
true,
true);
3327 plogf(
" Information about phases\n");
3328 plogf(
" PhaseName PhaseNum SingSpec GasPhase "
3329 " EqnState NumSpec");
3330 plogf(
" TMolesInert TKmoles\n");
3332 for (
size_t iphase = 0; iphase <
m_numPhases; iphase++) {
3345 plogf(
"\nElemental Abundances: ");
3346 plogf(
" Target_kmol ElemType ElActive\n");
3347 for (
size_t i = 0; i <
m_nelem; ++i) {
3348 writeline(
' ', 26,
false);
3354 plogf(
"\nChemical Potentials: (J/kmol)\n");
3355 plogf(
" Species (phase) "
3356 " SS0ChemPot StarChemPot\n");
3357 for (
size_t iphase = 0; iphase <
m_numPhases; iphase++) {
3360 for (
size_t kindex = 0; kindex < Vphase->
nSpecies(); kindex++) {
3373 writeline(
'=', 80,
true,
true);
3374 writeline(
'=', 20,
false);
3375 plogf(
" VCS_PROB: END OF PROBLEM STATEMENT ");
3386 for (
size_t eVP = 0; eVP < neVP; eVP++) {
3387 size_t foundPos =
npos;
3392 for (
size_t e = 0; e <
m_nelem; e++) {
3394 if (!strcmp(enVP.c_str(), en.c_str())) {
3399 if (foundPos ==
npos) {
3401 int elactive = volPhase->elementActive(eVP);
3402 size_t e =
addElement(enVP.c_str(), elType, elactive);
3412 throw CanteraError(
"VCS_SOLVE::addOnePhaseSpecies",
"Shouldn't be here");
3418 "element not found");
3432 "error: element must have a name");
Headers for the MultiPhase object that is used to set up multiphase equilibrium problems (see Chemica...
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
A class for 2D arrays stored in column-major (Fortran-compatible) form.
void zero()
Set all of the entries to zero.
size_t nColumns() const
Number of columns.
double * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
virtual void resize(size_t n, size_t m, double v=0.0)
Resize the array, and fill the new entries with 'v'.
Base class for exceptions thrown by Cantera classes.
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
A class for multiphase mixtures.
double speciesMoles(size_t kGlob) const
Returns the moles of global species k. units = kmol.
double pressure() const
Pressure [Pa].
double temperature() const
Temperature [K].
string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
double volume() const
The total mixture volume [m^3].
void uploadMoleFractionsFromPhases()
Update the locally-stored composition within this object to match the current compositions of the pha...
size_t nElements() const
Number of elements.
ThermoPhase & phase(size_t n)
Return a reference to phase n.
void setPhaseMoles(const size_t n, const double moles)
Set the number of moles of phase with index n.
double elementMoles(size_t m) const
Total moles of global element m, summed over all phases.
A species thermodynamic property manager for a phase.
virtual int reportType(size_t index) const
This utility function reports the type of parameterization used for the species with index number ind...
virtual void reportParams(size_t index, int &type, double *const c, double &minTemp, double &maxTemp, double &refPressure) const
This utility function reports back the type of parameterization and all of the parameters for the spe...
size_t nSpecies() const
Returns the number of species in the phase.
size_t nDim() const
Returns the number of spatial dimensions (1, 2, or 3)
size_t nElements() const
Number of elements.
double molecularWeight(size_t k) const
Molecular weight of species k.
double charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
string name() const
Return the name of the phase.
Base class for a phase with thermodynamic properties.
double electricPotential() const
Returns the electric potential of this phase (V).
string type() const override
String indicating the thermodynamic model implemented.
virtual int activityConvention() const
This method returns the convention used in specification of the activities, of which there are curren...
virtual MultiSpeciesThermo & speciesThermo(int k=-1)
Return a changeable reference to the calculation manager for species reference-state thermodynamic pr...
Class to keep track of time and iterations.
int Its
Current number of iterations in the main loop of vcs_TP() to solve for thermo equilibrium.
int T_Calls_Inest
Current number of times the initial thermo equilibrium estimator has been called.
double T_Time_basopt
Total Time spent in basopt.
double T_Time_inest
Time spent in initial estimator.
int T_Basis_Opts
Total number of optimizations of the components basis set done.
double T_Time_vcs_TP
Current time spent in vcs_TP.
int Basis_Opts
number of optimizations of the components basis set done
double T_Time_vcs
Time spent in the vcs suite of programs.
int T_Its
Total number of iterations in the main loop of vcs_TP() to solve for thermo equilibrium.
int T_Calls_vcs_TP
Current number of calls to vcs_TP.
double Time_basopt
Current Time spent in basopt.
double Time_vcs_TP
Current time spent in vcs_TP.
size_t vcs_popPhaseID(vector< size_t > &phasePopPhaseIDs)
Decision as to whether a phase pops back into existence.
vector< double > m_deltaPhaseMoles
Change in the total moles in each phase.
size_t m_numPhases
Number of Phases in the problem.
int m_doEstimateEquil
Setting for whether to do an initial estimate.
vector< double > m_deltaMolNumSpecies
Reaction Adjustments for each species during the current step.
VCS_SOLVE(MultiPhase *mphase, int printLvl=0)
Initialize the sizes within the VCS_SOLVE object.
vector< double > m_spSize
total size of the species
vector< double > m_scSize
Absolute size of the stoichiometric coefficients.
int vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
Evaluate the standard state free energies at the current temperature and pressure.
vector< double > m_wtSpecies
Molecular weight of each species.
vector< size_t > m_indexRxnToSpecies
Mapping between the species index for noncomponent species and the full species index.
int vcs_inest_TP()
Create an initial estimate of the solution to the thermodynamic equilibrium problem.
vector< size_t > m_phaseID
Mapping from the species number to the phase number.
vector< string > m_speciesName
Species string name for the kth species.
void vcs_deltag(const int L, const bool doDeleted, const int vcsState, const bool alterZeroedPhases=true)
This subroutine calculates reaction free energy changes for all noncomponent formation reactions.
double m_tolmaj2
Below this, major species aren't refined any more.
double vcs_VolTotal(const double tkelvin, const double pres, const double w[], double volPM[])
Calculation of the total volume and the partial molar volumes.
int vcs_elem_rearrange(double *const aw, double *const sa, double *const sm, double *const ss)
Rearrange the constraint equations represented by the Formula Matrix so that the operational ones are...
vector< double > m_deltaGRxn_new
Delta G(irxn) for the noncomponent species in the mechanism.
Array2D m_np_dLnActCoeffdMolNum
Change in the log of the activity coefficient with respect to the mole number multiplied by the phase...
size_t m_nelem
Number of element constraints in the problem.
vector< double > m_PMVolumeSpecies
Partial molar volumes of the species.
vector< int > m_elementActive
Specifies whether an element constraint is active.
int m_useActCoeffJac
Choice of Hessians.
vector< size_t > m_speciesLocalPhaseIndex
Index that keeps track of the index of the species within the local phase.
double m_totalVol
Total volume of all phases. Units are m^3.
int vcs_basopt(const bool doJustComponents, double aw[], double sa[], double sm[], double ss[], double test, bool *const usedZeroedSpecies)
Choose the optimum species basis for the calculations.
int vcs_prep(int printLvl)
This routine is mostly concerned with changing the private data to be consistent with what's needed f...
vector< string > m_elementName
Vector of strings containing the element names.
vector< int > m_speciesUnknownType
Specifies the species unknown type.
int vcs_elcorr(double aa[], double x[])
This subroutine corrects for element abundances.
vector< double > m_TmpPhase
Temporary vector of length NPhase.
vector< size_t > m_speciesMapIndex
Index vector that keeps track of the species vector rearrangement.
void vcs_TCounters_report(int timing_print_lvl=1)
Create a report on the plog file containing timing and its information.
Array2D m_phaseParticipation
This is 1 if the phase, iphase, participates in the formation reaction irxn, and zero otherwise.
static void disableTiming()
Disable printing of timing information.
vector< int > m_actConventionSpecies
specifies the activity convention of the phase containing the species
double vcs_phaseStabilityTest(const size_t iph)
Main program to test whether a deleted phase should be brought back into existence.
vector< double > m_elemAbundances
Element abundances vector.
vector< double > m_tPhaseMoles_new
total kmols of species in each phase in the tentative soln vector
void vcs_inest(double *const aw, double *const sa, double *const sm, double *const ss, double test)
Estimate equilibrium compositions.
vector< double > m_actCoeffSpecies_old
Molar-based Activity Coefficients for Species based on old mole numbers.
vector< double > m_SSfeSpecies
Standard state chemical potentials for species K at the current temperature and pressure.
bool vcs_popPhasePossible(const size_t iphasePop) const
Utility function that evaluates whether a phase can be popped into existence.
int vcs_TP(int ipr, int ip1, int maxit, double T, double pres)
Solve an equilibrium problem at a particular fixed temperature and pressure.
vector< double > m_lnMnaughtSpecies
specifies the ln(Mnaught) used to calculate the chemical potentials
vector< int > m_elType
Type of the element constraint.
void vcs_elab()
Computes the current elemental abundances vector.
void prob_report(int print_lvl)
Print out the problem specification in all generality as it currently exists in the VCS_SOLVE object.
int m_timing_print_lvl
printing level of timing information
int vcs_setMolesLinProg()
Estimate the initial mole numbers by constrained linear programming.
vector< double > m_molNumSpecies_old
Total moles of the species.
vector< double > TPhInertMoles
Total kmoles of inert to add to each phase.
double vcs_Total_Gibbs(double *w, double *fe, double *tPhMoles)
Calculate the total dimensionless Gibbs free energy.
void vcs_CalcLnActCoeffJac(const double *const moleSpeciesVCS)
Recalculate all of the activity coefficients in all of the phases based on input mole numbers.
vector< size_t > m_elementMapIndex
Index vector that keeps track of the rearrangement of the elements.
void vcs_elabPhase(size_t iphase, double *const elemAbundPhase)
Computes the elemental abundances vector for a single phase, elemAbundPhase[], and returns it through...
size_t m_numSpeciesRdc
Current number of species in the problems.
vector< double > m_feSpecies_old
Free energy vector from the start of the current iteration.
int m_debug_print_lvl
Debug printing lvl.
void vcs_counters_init(int ifunc)
Initialize the internal counters.
void vcs_reinsert_deleted(size_t kspec)
We make decisions on the initial mole number, and major-minor status here.
vector< double > m_chargeSpecies
Charge of each species. Length = number of species.
Array2D m_formulaMatrix
Formula matrix for the problem.
void addPhaseElements(vcs_VolPhase *volPhase)
Add elements to the local element list.
int vcs_solve_TP(int print_lvl, int printDetails, int maxit)
Main routine that solves for equilibrium at constant T and P using a variant of the VCS method.
size_t addOnePhaseSpecies(vcs_VolPhase *volPhase, size_t k, size_t kT)
This routines adds entries for the formula matrix for one species.
Array2D m_stoichCoeffRxnMatrix
Stoichiometric coefficient matrix for the reaction mechanism expressed in Reduced Canonical Form.
size_t m_numRxnRdc
Current number of non-component species in the problem.
VCS_COUNTERS * m_VCount
Timing and iteration counters for the vcs object.
int vcs(int ipr, int ip1, int maxit)
Solve an equilibrium problem.
vector< double > m_deltaGRxn_Deficient
Last deltag[irxn] from the previous step with additions for possible births of zeroed phases.
vector< unique_ptr< vcs_VolPhase > > m_VolPhaseList
Array of Phase Structures. Length = number of phases.
double vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal)
Calculates the diagonal contribution to the Hessian due to the dependence of the activity coefficient...
double vcs_GibbsPhase(size_t iphase, const double *const w, const double *const fe)
Calculate the total dimensionless Gibbs free energy of a single phase.
size_t m_numComponents
Number of components calculated for the problem.
void vcs_prob_specifyFully()
Fully specify the problem to be solved.
double m_pressurePA
Pressure.
void vcs_SSPhase()
Calculate the status of single species phases.
int vcs_report(int iconv)
Print out a report on the state of the equilibrium problem to standard output.
vector< double > m_feSpecies_new
Dimensionless new free energy for all the species in the mechanism at the new tentative T,...
double m_Faraday_dim
dimensionless value of Faraday's constant, F / RT (1/volt)
size_t m_nsp
Total number of species in the problems.
vector< double > m_actCoeffSpecies_new
Molar-based Activity Coefficients for Species.
bool vcs_elabcheck(int ibound)
Checks to see if the element abundances are in compliance.
vector< unique_ptr< VCS_SPECIES_THERMO > > m_speciesThermoList
Vector of pointers to thermo structures which identify the model and parameters for evaluating the th...
vector< char > m_SSPhase
Boolean indicating whether a species belongs to a single-species phase.
void vcs_delete_memory()
Delete memory that isn't just resizable STL containers.
vector< double > m_phasePhi
electric potential of the iph phase
double vcs_Hessian_actCoeff_diag(size_t irxn)
Calculates the diagonal contribution to the Hessian due to the dependence of the activity coefficient...
vector< double > m_elemAbundancesGoal
Element abundances vector Goals.
void vcs_prob_update()
Transfer the results of the equilibrium calculation back from VCS_SOLVE.
size_t vcs_RxnStepSizes(int &forceComponentCalc, size_t &kSpecial)
Calculates formation reaction step sizes.
Array2D m_deltaMolNumPhase
Change in the number of moles of phase, iphase, due to the noncomponent formation reaction,...
vector< int > m_phaseActConvention
specifies the activity convention of the phase.
vector< int > m_speciesStatus
Major -Minor status vector for the species in the problem.
vector< double > m_deltaGRxn_tmp
Temporary vector of Rxn DeltaG's.
vector< double > m_TmpPhase2
Temporary vector of length NPhase.
void vcs_switch_elem_pos(size_t ipos, size_t jpos)
Swaps the indices for all of the global data for two elements, ipos and jpos.
double vcs_tmoles()
Calculates the total number of moles of species in all phases.
int m_printLvl
Print level for print routines.
void vcs_dfe(const int stateCalc, const int ll, const size_t lbot, const size_t ltop)
Calculate the dimensionless chemical potentials of all species or of certain groups of species,...
vector< double > m_tPhaseMoles_old
Total kmols of species in each phase.
int vcs_popPhaseRxnStepSizes(const size_t iphasePop)
Calculates the deltas of the reactions due to phases popping into existence.
size_t m_numRxnTot
Total number of non-component species in the problem.
vector< double > m_molNumSpecies_new
Tentative value of the mole number vector.
void vcs_fePrep_TP()
Initialize the chemical potential of single species phases.
double m_temperature
Temperature (Kelvin)
vector< double > m_deltaGRxn_old
Last deltag[irxn] from the previous step.
size_t addElement(const char *elNameNew, int elType, int elactive)
This routine resizes the number of elements in the VCS_SOLVE object by adding a new element to the en...
double m_totalMolNum
Total number of kmoles in all phases.
Identifies the thermo model for the species.
size_t IndexSpeciesPhase
Index of this species in the current phase.
int SSStar_Vol_Model
Models for the standard state volume of each species.
double SS0_Cp0
Base heat capacity used in the VCS_SS0_CONSTANT_CP model.
double SSStar_Vol0
parameter that is used in the VCS_SSVOL_CONSTANT model.
double SS0_H0
Base enthalpy used in the VCS_SS0_CONSTANT_CP model.
double SS0_T0
Base temperature used in the VCS_SS0_CONSTANT_CP model.
double SS0_TSave
Internal storage of the last temperature used in the calculation of the reference naught Gibbs free e...
int SSStar_Model
Integer value representing the star state model.
double SS0_S0
Base entropy used in the VCS_SS0_CONSTANT_CP model.
int SS0_Model
Integer representing the models for the species standard state Naught temperature dependence.
double SS0_feSave
Internal storage of the last calculation of the reference naught Gibbs free energy at SS0_TSave.
size_t IndexPhase
Index of the phase that this species belongs to.
vcs_VolPhase * OwningPhase
Pointer to the owning phase object.
The class provides the wall clock timer in seconds.
double secondsWC()
Returns the wall clock time in seconds since the last reset.
Properties of a single species.
double VolPM
Partial molar volume of the species.
int SurfaceSpecies
True if this species belongs to a surface phase.
string SpName
Name of the species.
double WtSpecies
Molecular Weight of the species (gm/mol)
VCS_SPECIES_THERMO * SpeciesThermo
Pointer to the thermo structure for this species.
vector< double > FormulaMatrixCol
Column of the formula matrix, comprising the element composition of the species.
double Charge
Charge state of the species -> This may be duplication of what's in the FormulaMatrixCol entries.
Phase information and Phase calculations for vcs.
void setCreationMoleNumbers(const double *const n_k, const vector< size_t > &creationGlobalRxnNumbers)
Sets the creationMoleNum's within the phase object.
void setElectricPotential(const double phi)
set the electric potential of the phase
int speciesUnknownType(const size_t k) const
Returns the type of the species unknown.
size_t elemGlobalIndex(const size_t e) const
Returns the global index of the local element index for the phase.
size_t nSpecies() const
Return the number of species in the phase.
string PhaseName
String name for the phase.
void setMoleFractionsState(const double molNum, const double *const moleFracVec, const int vcsStateStatus)
Set the moles and/or mole fractions within the phase.
size_t VP_ID_
Original ID of the phase in the problem.
const Array2D & getFormulaMatrix() const
Get a constant form of the Species Formula Matrix.
vcs_SpeciesProperties * speciesProperty(const size_t kindex)
Retrieve the kth Species structure for the species belonging to this phase.
size_t phiVarIndex() const
Return the index of the species that represents the the voltage of the phase.
double totalMolesInert() const
Returns the value of the total kmol of inert in the phase.
bool m_singleSpecies
If true, this phase consists of a single species.
void setMolesFromVCS(const int stateCalc, const double *molesSpeciesVCS=0)
Set the moles within the phase.
void setElemGlobalIndex(const size_t eLocal, const size_t eGlobal)
sets a local phase element to a global index value
void resize(const size_t phaseNum, const size_t numSpecies, const size_t numElem, const char *const phaseName, const double molesInert=0.0)
The resize() function fills in all of the initial information if it is not given in the constructor.
int exists() const
int indicating whether the phase exists or not
size_t nElemConstraints() const
Returns the number of element constraints.
size_t spGlobalIndexVCS(const size_t spIndex) const
Return the Global VCS index of the kth species in the phase.
size_t transferElementsFM(const ThermoPhase *const tPhase)
Transfer all of the element information from the ThermoPhase object to the vcs_VolPhase object.
double GStar_calc_one(size_t kspec) const
Gibbs free energy calculation for standard state of one species.
bool m_gasPhase
If true, this phase is a gas-phase like phase.
void setTotalMoles(const double totalMols)
Sets the total moles in the phase.
int p_activityConvention
Convention for the activity formulation.
void sendToVCS_GStar(double *const gstar) const
Fill in the standard state Gibbs free energy vector for VCS.
void sendToVCS_ActCoeff(const int stateCalc, double *const AC)
Fill in an activity coefficients vector within a VCS_SOLVE object.
double sendToVCS_VolPM(double *const VolPM) const
Fill in the partial molar volume vector for VCS.
void sendToVCS_LnActCoeffJac(Array2D &LnACJac_VCS)
Downloads the ln ActCoeff Jacobian into the VCS version of the ln ActCoeff Jacobian.
int m_eqnState
Type of the equation of state.
double G0_calc_one(size_t kspec) const
Gibbs free energy calculation at a temperature for the reference state of a species,...
void setState_TP(const double temperature_Kelvin, const double pressure_PA)
Sets the temperature and pressure in this object and underlying ThermoPhase objects.
string eos_name() const
Return the name corresponding to the equation of state.
int elementType(const size_t e) const
Type of the element constraint with index e.
void setSpGlobalIndexVCS(const size_t spIndex, const size_t spGlobalIndex)
set the Global VCS index of the kth species in the phase
void setExistence(const int existence)
Set the existence flag in the object.
const vector< double > & creationMoleNumbers(vector< size_t > &creationGlobalRxnNumbers) const
Return a const reference to the creationMoleNumbers stored in the object.
string elementName(const size_t e) const
Name of the element constraint with index e.
bool isIdealSoln() const
Returns whether the phase is an ideal solution phase.
double totalMoles() const
Return the total moles in the phase.
void setPtrThermoPhase(ThermoPhase *tp_ptr)
Set the pointer for Cantera's ThermoPhase parameter.
Declarations for a simple class that implements an Ansi C wall clock timer (see Cantera::clockWC).
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
void debuglog(const string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
void writelogendl()
Write an end of line character to the screen and flush output.
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 Faraday
Faraday constant [C/kmol].
const double GasConstant
Universal Gas Constant [J/kmol/K].
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
int vcs_timing_print_lvl
Global hook for turning on and off time printing.
const char * vcs_speciesType_string(int speciesStatus, int length)
Returns a const char string representing the type of the species given by the first argument.
bool vcs_doubleEqual(double d1, double d2)
Simple routine to check whether two doubles are equal up to roundoff error.
double vcs_l2norm(const vector< double > &vec)
determine the l2 norm of a vector of doubles
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
Contains const definitions for types of species reference-state thermodynamics managers (see Species ...
#define SIMPLE
Constant Cp thermo.
Contains declarations for string manipulation functions within Cantera.
Header for the object representing each phase within vcs.
#define VCS_SPECIES_INTERFACIALVOLTAGE
Species refers to an electron in the metal.
#define VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
Unknown refers to the voltage level of a phase.
#define VCS_STATECALC_OLD
State Calculation based on the old or base mole numbers.
#define VCS_DELETE_MINORSPECIES_CUTOFF
Cutoff relative mole number value, below which species are deleted from the equilibrium problem.
#define VCS_SSVOL_IDEALGAS
Models for the standard state volume of each species.
#define VCS_SPECIES_COMPONENT
Species is a component which can be nonzero.
#define VCS_SPECIES_MAJOR
Species is a major species.
#define VCS_PHASE_EXIST_NO
Phase doesn't currently exist in the mixture.
#define VCS_PHASE_EXIST_ALWAYS
#define VCS_SPECIES_STOICHZERO
Species lies in a multicomponent phase that is active, but species concentration is zero due to stoic...
#define VCS_ELEM_TYPE_ABSPOS
Normal element constraint consisting of positive coefficients for the formula matrix.
#define VCS_STATECALC_NEW
State Calculation based on the new or tentative mole numbers.
#define VCS_SPECIES_TYPE_MOLNUM
Unknown refers to mole number of a single species.
#define VCS_SPECIES_ZEROEDPHASE
Species lies in a multicomponent phase that is zeroed atm.
#define VCS_ELEM_TYPE_LATTICERATIO
Constraint associated with maintaining a fixed lattice stoichiometry in the solids.
#define VCS_SPECIES_ACTIVEBUTZERO
Species lies in a multicomponent phase that is active, but species concentration is zero.
#define VCS_ELEM_TYPE_ELECTRONCHARGE
This refers to conservation of electrons.
#define VCS_SPECIES_ZEROEDSS
Species is a SS phase, that is currently zeroed out.
#define VCS_DELETE_PHASE_CUTOFF
Cutoff relative moles below which a phase is deleted from the equilibrium problem.
#define VCS_SMALL_MULTIPHASE_SPECIES
Relative value of multiphase species mole number for a multiphase species which is small.
#define VCS_ELEM_TYPE_CHARGENEUTRALITY
This refers to a charge neutrality of a single phase.
#define VCS_SPECIES_MINOR
Species is a major species.
#define VCS_DELETE_ELEMENTABS_CUTOFF
Cutoff moles below which a phase or species which comprises the bulk of an element's total concentrat...
#define VCS_PHASE_EXIST_YES
Phase is a normal phase that currently exists.
#define VCS_STATECALC_PHASESTABILITY
State Calculation based on tentative mole numbers for a phase which is currently zeroed,...
#define plogf
define this Cantera function to replace printf
Header file for the internal object that holds the vcs equilibrium problem (see Class VCS_SOLVE and C...