21enum stages {MAIN, EQUILIB_CHECK, ELEM_ABUND_CHECK,
22 RECHECK_DELETED, RETURN_A, RETURN_B};
24const int anote_size = 128;
30void VCS_SOLVE::checkDelta1(
double*
const dsLocal,
31 double*
const delTPhMoles,
size_t kspec)
34 for (
size_t k = 0; k < kspec; k++) {
37 dchange[iph] += dsLocal[k];
40 for (
size_t iphase = 0; iphase <
m_numPhases; iphase++) {
42 if (!
vcs_doubleEqual(dchange[iphase]/denom, delTPhMoles[iphase]/denom)) {
43 throw CanteraError(
"VCS_SOLVE::checkDelta1",
44 "we have found a problem");
52 bool allMinorZeroedSpecies =
false;
55 int rangeErrorFound = 0;
56 bool giveUpOnElemAbund =
false;
57 int finalElemAbundAttempts = 0;
58 bool uptodate_minors =
true;
59 int forceComponentCalc = 1;
61 char ANOTE[anote_size];
64 if (printDetails > 0 && print_lvl == 0) {
75 m_aw.assign(
m_nsp, 0.0);
78 int solveFail =
false;
85 plogf(
"VCS CALCULATION METHOD\n\n ");
86 plogf(
"MultiPhase Object\n");
99 plogf(
"\n ELEMENTAL ABUNDANCES CORRECT");
100 plogf(
" FROM ESTIMATE Type\n\n");
101 for (
size_t i = 0; i <
m_nelem; ++i) {
102 writeline(
' ', 26,
false);
108 plogf(
"\n MODIFIED LINEAR PROGRAMMING ESTIMATE OF EQUILIBRIUM - forced\n");
110 plogf(
"\n MODIFIED LINEAR PROGRAMMING ESTIMATE OF EQUILIBRIUM - where necessary\n");
113 plogf(
"\n USER ESTIMATE OF EQUILIBRIUM\n");
115 plogf(
" Stan. Chem. Pot. in J/kmol\n");
116 plogf(
"\n SPECIES FORMULA VECTOR ");
117 writeline(
' ', 41,
false);
118 plogf(
" STAN_CHEM_POT EQUILIBRIUM_EST. Species_Type\n\n");
119 writeline(
' ', 20,
false);
120 for (
size_t i = 0; i <
m_nelem; ++i) {
125 for (
size_t i = 0; i <
m_nsp; ++i) {
127 for (
size_t j = 0; j <
m_nelem; ++j) {
131 writeline(
' ', std::max(55-
int(
m_nelem)*8, 0),
false);
143 for (
size_t i = 0; i <
m_nsp; ++i) {
145 plogf(
"On Input species %-12s has a negative MF, setting it small\n",
166 if (forceComponentCalc) {
167 int retn = solve_tp_component_calc(allMinorZeroedSpecies);
172 forceComponentCalc = 0;
181 solve_tp_inner(iti, it1, uptodate_minors, allMinorZeroedSpecies,
182 forceComponentCalc, stage, printDetails > 0, ANOTE);
184 }
else if (stage == EQUILIB_CHECK) {
186 solve_tp_equilib_check(allMinorZeroedSpecies, uptodate_minors,
187 giveUpOnElemAbund, solveFail, iti, it1,
189 }
else if (stage == ELEM_ABUND_CHECK) {
191 solve_tp_elem_abund_check(iti, stage, lec, giveUpOnElemAbund,
192 finalElemAbundAttempts, rangeErrorFound);
193 }
else if (stage == RECHECK_DELETED) {
217 }
else if (stage == RETURN_A) {
233 }
else if (stage == RETURN_B) {
240 plogf(
" --- add_all_deleted(): some rxns not converged. RETURNING TO LOOP!\n");
256 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
274 if (rangeErrorFound) {
291int VCS_SOLVE::solve_tp_component_calc(
bool& allMinorZeroedSpecies)
293 double test = -1.0e-10;
294 bool usedZeroedSpecies;
295 int retn =
vcs_basopt(
false, &m_aw[0], &m_sa[0], &m_sm[0], &m_ss[0],
296 test, &usedZeroedSpecies);
323void VCS_SOLVE::solve_tp_inner(
size_t& iti,
size_t& it1,
324 bool& uptodate_minors,
325 bool& allMinorZeroedSpecies,
326 int& forceComponentCalc,
327 int& stage,
bool printDetails,
char* ANOTE)
336 if (!uptodate_minors) {
341 uptodate_minors =
true;
343 uptodate_minors =
false;
349 plogf(
" Iteration = %3d, Iterations since last evaluation of "
350 "optimal basis = %3d",
353 plogf(
" (all species)\n");
355 plogf(
" (only major species)\n");
381 vector<size_t> phasePopPhaseIDs(0);
383 if (iphasePop !=
npos) {
388 plogf(
" --- vcs_popPhaseRxnStepSizes() was called but stoich "
389 "prevented phase %d popping\n");
397 size_t iphaseDelete =
npos;
399 if (iphasePop ==
npos) {
405 plogf(
" --- vcs_RxnStepSizes not called because alternative"
406 "phase creation delta was used instead\n");
408 size_t doPhaseDeleteKspec =
npos;
409 size_t doPhaseDeleteIph =
npos;
428 if (iphaseDelete !=
npos) {
430 for (
size_t k = 0; k <
m_nsp; k++) {
439 throw CanteraError(
"VCS_SOLVE::solve_tp_inner",
440 "we shouldn't be here!");
460 plogf(
" --- Main Loop Treatment of each non-component species ");
462 plogf(
"- Full Calculation:\n");
464 plogf(
"- Major Components Calculation:\n");
466 plogf(
" --- Species IC ");
467 plogf(
" KMoles Tent_KMoles Rxn_Adj | Comment \n");
469 for (
size_t irxn = 0; irxn <
m_numRxnRdc; irxn++) {
477 if (iphasePop !=
npos) {
478 if (iph == iphasePop) {
481 snprintf(ANOTE, anote_size,
"Phase pop");
495 plogf(
" --- %s currently zeroed (SpStatus=%-2d):",
497 plogf(
"%3d DG = %11.4E WT = %11.4E W = %11.4E DS = %11.4E\n",
505 snprintf(ANOTE, anote_size,
"Species stays zeroed: DG = %11.4E",
509 snprintf(ANOTE, anote_size,
510 "Species stays zeroed even though dg neg due to "
511 "STOICH/PHASEPOP constraint: DG = %11.4E",
514 snprintf(ANOTE, anote_size,
"Species stays zeroed even"
515 " though dg neg: DG = %11.4E, ds zeroed",
520 for (
size_t j = 0; j <
m_nelem; ++j) {
524 if (atomComp > 0.0) {
527 snprintf(ANOTE, anote_size,
"Species stays zeroed"
528 " even though dG neg, because of %s elemAbund",
542 plogf(
" --- Zeroed species changed to major: ");
546 allMinorZeroedSpecies =
false;
549 plogf(
" --- Zeroed species changed to minor: ");
562 snprintf(ANOTE, anote_size,
"Born:IC=-1 to IC=1:DG=%11.4E",
578 snprintf(ANOTE, anote_size,
"minor species not considered");
580 plogf(
" --- %-12s%3d% 11.4E %11.4E %11.4E | %s\n",
607 plogf(
" --- Delete minor species in multispec phase: %-12s\n",
618 stage = RECHECK_DELETED;
632 snprintf(ANOTE, anote_size,
"Normal Major Calc");
641 snprintf(ANOTE, anote_size,
"major species is converged");
643 plogf(
" --- %-12s %3d %11.4E %11.4E %11.4E | %s\n",
662 snprintf(ANOTE, anote_size,
"dx set to 0, DG flipped sign due to "
663 "changed initial point");
674 snprintf(ANOTE, anote_size,
"initial nonpos kmoles= %11.3E",
702 if (sc_irxn[j] != 0.0) {
714 snprintf(ANOTE, anote_size,
715 "zeroing SS phase created a neg component species "
716 "-> reducing step size instead");
722 snprintf(ANOTE, anote_size,
"zeroing out SS phase: ");
729 doPhaseDeleteIph = iph;
731 doPhaseDeleteKspec = kspec;
733 plogf(
" --- SS species changed to zeroedss: ");
740 for (
size_t kk = 0; kk <
m_nsp; kk++) {
767 "VCS_SOLVE::solve_tp_inner",
768 "ds[kspec] = {}, dx = {}, kspec = {}\nwe have a problem!",
787 plogf(
" --- %-12.12s%3d %11.4E %11.4E %11.4E | %s\n",
793 if (doPhaseDeleteIph !=
npos) {
795 plogf(
" --- %-12.12s Main Loop Special Case deleting phase with species:\n",
805 plogf(
" c %11.4E %11.4E %11.4E |\n",
811 plogf(
" --- Finished Main Loop\n");
836 if (par <= 1.01 && par > 0.0) {
840 plogf(
" --- Reduction in step size due to component %s going negative = %11.3E\n",
843 for (
size_t i = 0; i <
m_nsp; ++i) {
860 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
864 throw CanteraError(
"VCS_SOLVE::solve_tp_inner",
865 "vcs_solve_TP: ERROR on step change wt[{}:{}]: {} < 0.0",
890 plogf(
" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n",
893 plogf(
" --- Total tentative Dimensionless Gibbs Free Energy = %20.13E\n",
901 if (printDetails && forced) {
902 plogf(
" -----------------------------------------------------\n");
903 plogf(
" --- FORCER SUBROUTINE changed the solution:\n");
904 plogf(
" --- SPECIES Status INIT MOLES TENT_MOLES");
905 plogf(
" FINAL KMOLES INIT_DEL_G/RT TENT_DEL_G/RT FINAL_DELTA_G/RT\n");
920 writeline(
' ', 26,
false);
921 plogf(
"Norms of Delta G():%14.6E%14.6E\n",
928 plogf(
" Total kmoles of liquid = %15.7E\n", 0.0);
930 plogf(
" Total New Dimensionless Gibbs Free Energy = %20.13E\n",
933 plogf(
" -----------------------------------------------------\n");
942 plogf(
" --- Summary of the Update ");
944 plogf(
" (all species):");
946 plogf(
" (only major species):");
949 plogf(
" --- Species Status Initial_KMoles Final_KMoles Initial_Mu/RT");
950 plogf(
" Mu/RT Init_Del_G/RT Delta_G/RT\n");
960 plogf(
" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n",
968 plogf(
" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n",
974 writeline(
' ', 56,
false);
975 plogf(
"Norms of Delta G():%14.6E%14.6E\n",
979 plogf(
" --- Phase_Name KMoles(after update)\n");
988 plogf(
" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n",
991 plogf(
" --- Total New Dimensionless Gibbs Free Energy = %20.13E\n",
1019 plogf(
" --- Increment counter increased, step is accepted: %4d\n",
1029 bool justDeletedMultiPhase =
false;
1034 plogf(
" --- Setting microscopic phase %d to zero\n", iph);
1036 justDeletedMultiPhase =
true;
1044 if (justDeletedMultiPhase) {
1045 bool usedZeroedSpecies;
1046 double test = -1.0e-10;
1047 int retn =
vcs_basopt(
false, &m_aw[0], &m_sa[0], &m_sm[0], &m_ss[0],
1048 test, &usedZeroedSpecies);
1050 throw CanteraError(
"VCS_SOLVE::solve_tp_inner",
1051 "BASOPT returned with an error condition");
1062 plogf(
" --- Normal element abundance check");
1071 uptodate_minors =
true;
1083 bool doSwap =
false;
1106 plogf(
" and share nonzero stoic: %-9.1f\n",
1109 forceComponentCalc = 1;
1115 stage = EQUILIB_CHECK;
1125 plogf(
" --- Reevaluate major-minor status of noncomponents:\n");
1128 for (
size_t irxn = 0; irxn <
m_numRxnRdc; irxn++) {
1133 plogf(
" --- major/minor species is now zeroed out: %s\n",
1140 plogf(
" --- Noncomponent turned from major to minor: ");
1142 plogf(
" --- Component turned into a minor species: ");
1144 plogf(
" --- Zeroed Species turned into a "
1154 plogf(
" --- Noncomponent turned from minor to major: ");
1156 plogf(
" --- Component turned into a major: ");
1158 plogf(
" --- Noncomponent turned from zeroed to major: ");
1173void VCS_SOLVE::solve_tp_equilib_check(
bool& allMinorZeroedSpecies,
1174 bool& uptodate_minors,
1175 bool& giveUpOnElemAbund,
int& solveFail,
1176 size_t& iti,
size_t& it1,
int maxit,
1177 int& stage,
bool& lec)
1179 if (! allMinorZeroedSpecies) {
1181 plogf(
" --- Equilibrium check for major species: ");
1183 for (
size_t irxn = 0; irxn <
m_numRxnRdc; ++irxn) {
1198 iti = ((it1/4) *4) - it1;
1206 debuglog(
" MAJOR SPECIES CONVERGENCE achieved "
1219 uptodate_minors =
true;
1222 plogf(
" --- Equilibrium check for minor species: ");
1224 for (
size_t irxn = 0; irxn <
m_numRxnRdc; ++irxn) {
1246 plogf(
" CONVERGENCE achieved\n");
1257 if (!giveUpOnElemAbund) {
1259 plogf(
" --- Check the Full Element Abundances: ");
1269 plogf(
" passed for NC but failed for NE: RANGE ERROR\n");
1273 stage = ELEM_ABUND_CHECK;
1284 stage = RECHECK_DELETED;
1293void VCS_SOLVE::solve_tp_elem_abund_check(
size_t& iti,
int& stage,
bool& lec,
1294 bool& giveUpOnElemAbund,
1295 int& finalElemAbundAttempts,
1296 int& rangeErrorFound)
1303 rangeErrorFound = 0;
1324 if (finalElemAbundAttempts >= 3) {
1325 giveUpOnElemAbund =
true;
1326 stage = EQUILIB_CHECK;
1328 finalElemAbundAttempts++;
1335 }
else if (ncAfter) {
1339 plogf(
" --- vcs_solve_tp: RANGE SPACE ERROR ENCOUNTERED\n");
1340 plogf(
" --- vcs_solve_tp: - Giving up on NE Element Abundance satisfaction \n");
1341 plogf(
" --- vcs_solve_tp: - However, NC Element Abundance criteria is satisfied \n");
1342 plogf(
" --- vcs_solve_tp: - Returning the calculated equilibrium condition\n");
1344 rangeErrorFound = 1;
1345 giveUpOnElemAbund =
true;
1350 stage = EQUILIB_CHECK;
1357 stage = EQUILIB_CHECK;
1369 if (w_kspec <= 0.0) {
1372 dg_irxn = std::max(dg_irxn, -200.0);
1374 snprintf(ANOTE, anote_size,
"minor species alternative calc");
1376 if (dg_irxn >= 23.0) {
1377 double molNum_kspec_new = w_kspec * 1.0e-10;
1384 return molNum_kspec_new - w_kspec;
1403 double a =
clip(w_kspec * s, -1.0+1e-8, 100.0);
1404 double tmp =
clip(-dg_irxn / (1.0 + a), -200.0, 200.0);
1405 double wTrial = w_kspec * exp(tmp);
1406 double molNum_kspec_new = wTrial;
1408 if (wTrial > 100. * w_kspec) {
1410 if (molNumMax < 100. * w_kspec) {
1411 molNumMax = 100. * w_kspec;
1413 if (wTrial > molNumMax) {
1414 molNum_kspec_new = molNumMax;
1416 molNum_kspec_new = wTrial;
1418 }
else if (1.0E10 * wTrial < w_kspec) {
1419 molNum_kspec_new= 1.0E-10 * w_kspec;
1421 molNum_kspec_new = wTrial;
1430 return molNum_kspec_new - w_kspec;
1437 snprintf(ANOTE, anote_size,
"voltage species alternative calc");
1448 "delete_species() ERROR: called for a component {}", kspec);
1452 double dx = *delta_ptr;
1456 double tmp = sc_irxn[j] * dx;
1479 double tmp = sc_irxn[j] * dx;
1502 plogf(
"vcs_zero_species: Couldn't zero the species %d, "
1503 "did delta of %g. orig conc of %g\n",
1523 "Failed to delete a species!");
1539 if (kspec != klast) {
1554 bool stillExists =
false;
1603 for (
size_t k = 0; k <
m_nsp; k++) {
1626 bool successful =
true;
1631 plogf(
" --- delete_multiphase %d, %s\n", iph, Vphase->
PhaseName);
1646 plogf(
" --- delete_multiphase %d, %s ERROR problems deleting species %s\n",
1648 plogf(
" --- delta attempted: %g achieved: %g "
1649 " Zeroing it manually\n", dx, dxTent);
1670 double dxPerm = 0.0, dxPerm2 = 0.0;
1674 plogf(
" --- delete_multiphase One of the species is a component %d - %s with mole number %g\n",
1687 if (jcomp != kcomp) {
1695 if (fabs(dxPerm2) < fabs(dxPerm)) {
1702 if (dxPerm != 0.0) {
1711 plogf(
" --- delete_multiphase One of the species is a component %d - %s still with mole number %g\n",
1713 plogf(
" --- zeroing it \n");
1737 plogf(
" an active but zeroed species because its phase "
1761 plogf(
" --- Start rechecking deleted species in multispec phases\n");
1785 xtcutoff[iph] = 0.0;
1843 for (
int cits = 0; cits < 3; cits++) {
1878 plogf(
" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
1885 plogf(
" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
1892 plogf(
" --- add_deleted(): species %s added back in with mol number %g\n",
1895 plogf(
" --- add_deleted(): species %s failed to be added back in\n");
1914 plogf(
" --- add_deleted(): species %s with mol number %g not converged: DG = %g\n",
1930 for (
size_t irxn = 0; irxn <
m_numRxnRdc; ++irxn) {
1940 for (
size_t irxn = 0; irxn <
m_numRxnRdc; ++irxn) {
1948 plogf(
" --- subroutine FORCE: Beginning Slope = %g\n", s1);
1949 plogf(
" --- subroutine FORCE: End Slope = %g\n", s2);
1952 if (s1 > 0.0 || s2 <= 0.0) {
1959 if (fabs(s1 -s2) > 1.0E-200) {
1960 al = s1 / (s1 - s2);
1962 if (al >= 0.95 || al < 0.0) {
1967 plogf(
" --- subroutine FORCE produced a damping factor = %g\n", al);
1986 plogf(
" --- subroutine FORCE adjusted the mole "
1987 "numbers, AL = %10.3f\n", al);
2002 for (
size_t irxn = 0; irxn <
m_numRxnRdc; ++irxn) {
2010 plogf(
" --- subroutine FORCE: Adj End Slope = %g\n", s2);
2016 double ss[],
double test,
bool*
const usedZeroedSpecies)
2020 size_t jlose =
npos;
2025 for (
size_t i=0; i<77; i++) {
2029 plogf(
" --- Subroutine BASOPT called to ");
2030 if (doJustComponents) {
2031 plogf(
"calculate the number of components\n");
2033 plogf(
"reevaluate the components\n");
2037 plogf(
" --- Formula Matrix used in BASOPT calculation\n");
2038 plogf(
" --- Active | ");
2039 for (
size_t j = 0; j <
m_nelem; j++) {
2043 plogf(
" --- Species | ");
2044 for (
size_t j = 0; j <
m_nelem; j++) {
2048 for (k = 0; k <
m_nsp; k++) {
2050 for (
size_t j = 0; j <
m_nelem; j++) {
2064 *usedZeroedSpecies =
false;
2065 vector<int> ipiv(ncTrial);
2072 for (k = 0; k <
m_nsp; k++) {
2082 while (jr < ncTrial) {
2130 *usedZeroedSpecies =
true;
2131 double maxConcPossKspec = 0.0;
2132 double maxConcPoss = 0.0;
2133 size_t kfound =
npos;
2134 int minNonZeroes = 100000;
2135 int nonZeroesKspec = 0;
2136 for (
size_t kspec = ncTrial; kspec <
m_nsp; kspec++) {
2138 maxConcPossKspec = 1.0E10;
2140 for (
size_t j = 0; j <
m_nelem; ++j) {
2149 if ((maxConcPossKspec >= maxConcPoss) || (maxConcPossKspec > 1.0E-5)) {
2150 if (nonZeroesKspec <= minNonZeroes) {
2151 if (kfound ==
npos || nonZeroesKspec < minNonZeroes) {
2161 minNonZeroes = std::min(minNonZeroes, nonZeroesKspec);
2162 maxConcPoss = std::max(maxConcPoss, maxConcPossKspec);
2166 if (kfound ==
npos) {
2169 for (
size_t kspec = ncTrial; kspec <
m_nsp; kspec++) {
2170 if (aw[kspec] >= 0.0) {
2171 size_t irxn = kspec - ncTrial;
2182 if (aw[k] == test) {
2187 throw CanteraError(
"VCS_SOLVE::vcs_basopt",
"we shouldn't be here");
2192 for (
size_t i = 0; i <
m_nsp; ++i) {
2196 plogf(
" --- Total number of components found = %3d (ne = %d)\n ",
2211 for (
size_t j = 0; j <
m_nelem; ++j) {
2218 for (
size_t j = 0; j < jl; ++j) {
2220 for (
size_t i = 0; i <
m_nelem; ++i) {
2227 for (
size_t j = 0; j < jl; ++j) {
2228 for (
size_t i = 0; i <
m_nelem; ++i) {
2237 for (
size_t ml = 0; ml <
m_nelem; ++ml) {
2238 sa[jr] += pow(sm[ml + jr*
m_nelem], 2);
2242 if (sa[jr] > 1.0e-6) {
2262 plogf(
" as component %3d\n", jr);
2265 std::swap(aw[jr], aw[k]);
2273 plogf(
" as component %3d\n", jr);
2282 if (doJustComponents) {
2312 C.
resize(ncTrial, ncTrial);
2313 for (
size_t j = 0; j < ncTrial; ++j) {
2314 for (
size_t i = 0; i < ncTrial; ++i) {
2320 for (
size_t j = 0; j < ncTrial; ++j) {
2332 for (
size_t j = 0; j <
m_nelem; j++) {
2337 for (
size_t j = 0; j <
m_nelem; j++) {
2342 for (k = 0; k <
m_nsp; k++) {
2344 for (
size_t j = 0; j < ncTrial; ++j) {
2345 for (
size_t i = 0; i < ncTrial; ++i) {
2355 for (
size_t j = 0; j < ncTrial; ++j) {
2365 size_t i = k - ncTrial;
2366 for (
size_t j = 0; j < ncTrial; j++) {
2375 for (
size_t j = 0; j < ncTrial; j++) {
2382 plogf(
" --- Components:");
2383 for (
size_t j = 0; j < ncTrial; j++) {
2386 plogf(
"\n --- Components Moles:");
2387 for (
size_t j = 0; j < ncTrial; j++) {
2389 plogf(
" % -10.3E", 0.0);
2394 plogf(
"\n --- NonComponent| Moles |");
2395 for (
size_t j = 0; j < ncTrial; j++) {
2403 plogf(
"|% -10.3E|", 0.0);
2407 for (
size_t j = 0; j < ncTrial; j++) {
2415 double sumMax = -1.0;
2421 for (
size_t j = 0; j < ncTrial; ++j) {
2424 for (
size_t n = 0; n < ncTrial; n++) {
2427 sum += coeff * numElements;
2431 for (
size_t n = 0; n < ncTrial; n++) {
2434 sum += coeff * numElements;
2437 if (fabs(sum) > sumMax) {
2445 if (fabs(sum) > 1.0E-6) {
2446 throw CanteraError(
"VCS_SOLVE::vcs_basopt",
"we have a prob");
2450 plogf(
" --- largest error in Stoich coeff = %g at rxn = %d ", sumMax, iMax);
2452 plogf(
" element = %d ", jMax);
2456 for (
size_t i=0; i<77; i++) {
2474 for (
size_t irxn = 0; irxn <
m_numRxnTot; ++irxn) {
2480 for (
size_t j = 0; j < ncTrial; ++j) {
2482 if (fabs(scrxn_ptr[j]) <= 1.0e-6) {
2516 double big = molNum[j] *
m_spSize[j] * 1.01;
2519 "spSize is nonpositive");
2521 for (
size_t i = j + 1; i < n; ++i) {
2524 "spSize is nonpositive");
2526 bool doSwap =
false;
2528 doSwap = (molNum[i] *
m_spSize[i]) > big;
2530 doSwap = molNum[i] > (molNum[largest] * 1.001);
2534 doSwap = (molNum[i] *
m_spSize[i]) > big;
2536 doSwap = molNum[i] > (molNum[largest] * 1.001);
2539 doSwap = (molNum[i] *
m_spSize[i]) > big;
2544 big = molNum[i] *
m_spSize[i] * 1.01;
2569 for (
size_t j = 0; j <
m_nelem; ++j) {
2572 if (atomComp > 0.0) {
2576 plogf(
" --- %s can not be nonzero because"
2577 " needed element %s is zero\n",
2597 if (stoicC != 0.0) {
2598 double negChangeComp = - stoicC;
2599 if (negChangeComp > 0.0) {
2602 plogf(
" --- %s is prevented from popping into existence because"
2603 " a needed component to be consumed, %s, has a zero mole number\n",
2612 }
else if (negChangeComp < 0.0) {
2615 plogf(
" --- %s is prevented from popping into existence because"
2616 " a needed component %s is in a zeroed-phase that would be "
2617 "popped into existence at the same time\n",
2704 const int ll,
const size_t lbot,
const size_t ltop)
2706 double* tPhMoles_ptr=0;
2707 double* actCoeff_ptr=0;
2708 double* feSpecies=0;
2722 "Subroutine vcs_dfe called with bad stateCalc value: {}", stateCalc);
2728 plogf(
" --- Subroutine vcs_dfe called for one species: ");
2731 plogf(
" --- Subroutine vcs_dfe called for all species");
2733 }
else if (ll > 0) {
2734 plogf(
" --- Subroutine vcs_dfe called for components and minors");
2736 plogf(
" --- Subroutine vcs_dfe called for components and majors");
2739 plogf(
" using tentative solution");
2750 tlogMoles[iph] = tPhInertMoles[iph];
2753 for (
size_t kspec = 0; kspec <
m_nsp; kspec++) {
2756 tlogMoles[iph] += molNum[kspec];
2761 "VCS_SOLVE::vcs_dfe",
2762 "phase Moles may be off, iph = {}, {} {}",
2763 iph, tlogMoles[iph], tPhMoles_ptr[iph]);
2767 if (tPhMoles_ptr[iph] > 0.0) {
2768 tlogMoles[iph] = log(tPhMoles_ptr[iph]);
2784 for (
size_t iphase = 0; iphase <
m_numPhases; iphase++) {
2799 for (
size_t kspec = l1; kspec < l2; ++kspec) {
2803 "We have an inconsistency!");
2805 "We have an unexpected situation!");
2819 if (tPhMoles_ptr[iph] > 0.0) {
2830 + log(actCoeff_ptr[kspec] * molNum[kspec])
2840 for (
size_t irxn = 0; irxn <
m_numRxnRdc; ++irxn) {
2846 "We have an inconsistency!");
2848 "We have an unexpected situation!");
2862 if (tPhMoles_ptr[iph] > 0.0) {
2873 + log(actCoeff_ptr[kspec] * molNum[kspec])
2881 }
else if (ll > 0) {
2883 for (
size_t irxn = 0; irxn <
m_numRxnRdc; ++irxn) {
2889 "We have an inconsistency!");
2891 "We have an unexpected situation!");
2905 if (tPhMoles_ptr[iph] > 0.0) {
2916 + log(actCoeff_ptr[kspec] * molNum[kspec])
2933 for (
size_t irxn = 0; irxn <
m_numRxnRdc; ++irxn) {
2936 dgLocal[irxn] < 0.0) {
2938 tmp += dgLocal[irxn] * dgLocal[irxn];
2950 for (
size_t i = 0; i <
m_nsp; i++) {
2969void VCS_SOLVE::check_tmoles()
const
2974 for (
size_t k = 0; k <
m_nsp; k++) {
2982 plogf(
"check_tmoles: we have found a problem with phase %d: %20.15g, %20.15g\n",
3002 "wrong stateCalc value: {}", vcsState);
3011 plogf(
" --- Species Status decision is reevaluated: All species are minor except for:\n");
3013 plogf(
" --- Species Status decision is reevaluated\n");
3015 for (
size_t kspec = 0; kspec <
m_nsp; ++kspec) {
3026 plogf(
"%s\n", sString);
3036 plogf(
" --- Purposely Zeroed-Phase Species (not in problem): %-s\n",
3049 plogf(
" --- Zeroed Species in an active MS phase (tmp): %-s\n",
3053 plogf(
" --- Zeroed Species in an active MS phase (Stoich Constraint): %-s\n",
3060 throw CanteraError(
"VCS_SOLVE::vcs_evaluate_speciesType",
3074 const int vcsState,
const bool alterZeroedPhases)
3083 double* molNumSpecies;
3084 double* actCoeffSpecies;
3096 throw CanteraError(
"VCS_SOLVE::vcs_deltag",
"bad vcsState");
3100 plogf(
" --- Subroutine vcs_deltag called for ");
3102 plogf(
"major noncomponents\n");
3103 }
else if (L == 0) {
3104 plogf(
"all noncomponents\n");
3106 plogf(
"minor noncomponents\n");
3112 for (
size_t irxn = 0; irxn <
m_numRxnRdc; ++irxn) {
3119 deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec];
3125 deltaGRxn[irxn] = std::max(0.0, deltaGRxn[irxn]);
3129 }
else if (L == 0) {
3131 for (
size_t irxn = 0; irxn < irxnl; ++irxn) {
3136 deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec];
3138 dtmp_ptr[kspec] < 0.0) {
3143 deltaGRxn[irxn] = std::max(0.0, deltaGRxn[irxn]);
3148 for (
size_t irxn = 0; irxn <
m_numRxnRdc; ++irxn) {
3155 deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec];
3157 dtmp_ptr[kspec] < 0.0) {
3162 deltaGRxn[irxn] = std::max(0.0, deltaGRxn[irxn]);
3207 if (alterZeroedPhases &&
false) {
3213 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
3216 sum += molNumSpecies[kspec];
3229 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
3234 deltaGRxn[irxn] =
clip(deltaGRxn[irxn], -50.0, 50.0);
3235 poly += exp(-deltaGRxn[irxn])/actCoeffSpecies[kspec];
3243 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
3247 deltaGRxn[irxn] = 1.0 - poly;
3255void VCS_SOLVE::vcs_printDeltaG(
const int stateCalc)
3270 bool zeroedPhase =
false;
3272 plogf(
" --- DELTA_G TABLE Components:");
3276 plogf(
"\n --- Components Moles:");
3280 plogf(
"\n --- NonComponent| Moles | ");
3299 for (
int i=0; i<77; i++) {
3305 writelog(
" --- DeltaG Table (J/kmol) Name PhID MoleNum MolFR "
3306 " ElectrChemStar ElectrChem DeltaGStar DeltaG(Pred) Stability\n");
3308 writeline(
'-', 132);
3310 for (
size_t kspec = 0; kspec <
m_nsp; kspec++) {
3315 double mfValue = 1.0;
3323 zeroedPhase =
false;
3325 if (tPhMoles_ptr[iphase] > 0.0) {
3329 mfValue = molNumSpecies[kspec] / tPhMoles_ptr[iphase];
3333 mfValue = Vphase->moleFraction(klocal);
3340 double feFull = feSpecies[kspec];
3343 feFull += log(actCoeff_ptr[kspec]) + log(mfValue);
3350 writelogf(
" % -12.4e", molNumSpecies[kspec]);
3353 writelogf(
" % -12.4e", feSpecies[kspec] * RT);
3356 writelogf(
" % -12.4e", deltaGRxn[irxn] * RT);
3357 writelogf(
" % -12.4e", (deltaGRxn[irxn] + feFull - feSpecies[kspec]) * RT);
3359 if (deltaGRxn[irxn] < 0.0) {
3360 if (molNumSpecies[kspec] > 0.0) {
3365 }
else if (deltaGRxn[irxn] > 0.0) {
3366 if (molNumSpecies[kspec] > 0.0) {
3378 writeline(
'-', 132);
3387 plogf(
"vcs_switch_pos: ifunc = 0: inappropriate args: %d %d\n",
3424 for (
size_t j = 0; j <
m_nelem; ++j) {
3428 for (
size_t i = 0; i <
m_nsp; i++) {
3431 for (
size_t i = 0; i <
m_nsp; i++) {
3443 plogf(
"switch_pos: ifunc = 1: inappropriate noncomp values: %d %d\n",
3461void VCS_SOLVE::vcs_setFlagsVolPhases(
const bool upToDate,
const int stateCalc)
3474void VCS_SOLVE::vcs_setFlagsVolPhase(
const size_t iph,
const bool upToDate,
3475 const int stateCalc)
Headers for the DenseMatrix object, which deals with dense rectangular matrices and description of th...
void zero()
Set all of the entries to zero.
double * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Base class for exceptions thrown by Cantera classes.
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
int Its
Current number of iterations in the main loop of vcs_TP() to solve for thermo equilibrium.
double T_Time_basopt
Total Time spent in basopt.
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
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.
vector< double > m_spSize
total size of the species
vector< double > m_scSize
Absolute size of the stoichiometric coefficients.
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.
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.
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...
void vcs_updateMolNumVolPhases(const int stateCalc)
Update all underlying vcs_VolPhase objects.
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.
void vcs_updateVP(const int stateCalc)
This routine uploads the state of the system into all of the vcs_VolumePhase objects in the current p...
size_t vcs_add_all_deleted()
Provide an estimate for the deleted species in phases that are not zeroed out.
bool vcs_globStepDamp()
This routine optimizes the minimization of the total Gibbs free energy by making sure the slope of th...
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.
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.
size_t vcs_basisOptMax(const double *const molNum, const size_t j, const size_t n)
Choose a species to test for the next component.
vector< size_t > m_speciesMapIndex
Index vector that keeps track of the species vector rearrangement.
Array2D m_phaseParticipation
This is 1 if the phase, iphase, participates in the formation reaction irxn, and zero otherwise.
vector< int > m_actConventionSpecies
specifies the activity convention of the phase containing the species
vector< double > m_elemAbundances
Element abundances vector.
vector< double > m_tPhaseMoles_new
total kmols of species in each phase in the tentative soln vector
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.
double l2normdg(double dg[]) const
Calculate the norm of a deltaGibbs free energy vector.
int delta_species(const size_t kspec, double *const delta_ptr)
Change the concentration of a species by delta moles.
double m_tolmin2
Below this, minor species aren't refined any more.
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.
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.
size_t m_numRxnMinorZeroed
Number of active species which are currently either treated as minor species.
size_t m_numSpeciesRdc
Current number of species in the problems.
int vcs_delete_species(const size_t kspec)
Change a single species from active to inactive status.
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.
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.
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.
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.
size_t m_numComponents
Number of components calculated for the problem.
int vcs_species_type(const size_t kspec) const
Evaluate the species category for the indicated species.
double m_pressurePA
Pressure.
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)
int vcs_recheck_deleted()
Recheck deleted species in multispecies phases.
double m_tolmin
Tolerance requirements for minor species.
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.
bool vcs_delete_multiphase(const size_t iph)
This routine handles the bookkeeping involved with the deletion of multiphase phases from the problem...
vector< double > m_phasePhi
electric potential of the iph phase
int vcs_zero_species(const size_t kspec)
Zero out the concentration of a species.
vector< double > m_elemAbundancesGoal
Element abundances vector Goals.
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_speciesStatus
Major -Minor status vector for the species in the problem.
vector< double > m_deltaGRxn_tmp
Temporary vector of Rxn DeltaG's.
double vcs_minor_alt_calc(size_t kspec, size_t irxn, bool *do_delete, char *ANOTE=0) const
Minor species alternative calculation.
double vcs_tmoles()
Calculates the total number of moles of species in all phases.
void vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k2)
Swaps the indices for all of the global data for two species, k1 and k2.
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.
double m_temperature
Temperature (Kelvin)
double m_tolmaj
Tolerance requirement for major species.
vector< double > m_deltaGRxn_old
Last deltag[irxn] from the previous step.
bool vcs_evaluate_speciesType()
This routine evaluates the species type for all species.
double m_totalMolNum
Total number of kmoles in all phases.
The class provides the wall clock timer in seconds.
double secondsWC()
Returns the wall clock time in seconds since the last reset.
Phase information and Phase calculations for vcs.
double electricPotential() const
Returns the electric field of the phase.
size_t nSpecies() const
Return the number of species in the phase.
string PhaseName
String name for the phase.
void setMolesFromVCSCheck(const int vcsStateStatus, const double *molesSpeciesVCS, const double *const TPhMoles)
Set the moles within the phase.
void updateFromVCS_MoleNumbers(const int stateCalc)
Update the moles within the phase, if necessary.
bool m_singleSpecies
If true, this phase consists of a single species.
int exists() const
int indicating whether the phase exists or not
size_t spGlobalIndexVCS(const size_t spIndex) const
Return the Global VCS index of the kth species in the phase.
void setTotalMoles(const double totalMols)
Sets the total moles in the phase.
void sendToVCS_ActCoeff(const int stateCalc, double *const AC)
Fill in an activity coefficients vector within a VCS_SOLVE object.
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.
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 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.
void writelogendl()
Write an end of line character to the screen and flush output.
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].
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
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.
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.
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_RELDELETE_SPECIES_CUTOFF
Cutoff relative mole fraction value, below which species are deleted from the equilibrium problem.
#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_ZEROEDMS
Species lies in a multicomponent phase with concentration zero.
#define VCS_SPECIES_ZEROEDPHASE
Species lies in a multicomponent phase that is zeroed atm.
#define VCS_SPECIES_ACTIVEBUTZERO
Species lies in a multicomponent phase that is active, but species concentration is zero.
#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_SPECIES_DELETED
Species has such a small mole fraction it is deleted even though its phase may possibly exist.
#define VCS_SPECIES_MINOR
Species is a major species.
#define VCS_PHASE_EXIST_YES
Phase is a normal phase that currently exists.
#define VCS_PHASE_EXIST_ZEROEDPHASE
Phase currently is zeroed due to a programmatic issue.
#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...