20 enum stages {MAIN, EQUILIB_CHECK, ELEM_ABUND_CHECK,
21 RECHECK_DELETED, RETURN_A, RETURN_B};
27 void VCS_SOLVE::checkDelta1(
double*
const dsLocal,
28 double*
const delTPhMoles,
size_t kspec)
31 for (
size_t k = 0; k < kspec; k++) {
33 size_t iph = m_phaseID[k];
34 dchange[iph] += dsLocal[k];
37 for (
size_t iphase = 0; iphase < m_numPhases; iphase++) {
38 double denom = max(m_totalMolNum, 1.0E-4);
39 if (!
vcs_doubleEqual(dchange[iphase]/denom, delTPhMoles[iphase]/denom)) {
40 throw CanteraError(
"VCS_SOLVE::checkDelta1",
41 "we have found a problem");
46 int VCS_SOLVE::vcs_solve_TP(
int print_lvl,
int printDetails,
int maxit)
49 bool allMinorZeroedSpecies =
false;
52 int rangeErrorFound = 0;
53 bool giveUpOnElemAbund =
false;
54 int finalElemAbundAttempts = 0;
55 bool uptodate_minors =
true;
56 int forceComponentCalc = 1;
60 m_debug_print_lvl = printDetails;
61 if (printDetails > 0 && print_lvl == 0) {
69 m_sm.assign(m_nelem * m_nelem, 0.0);
70 m_ss.assign(m_nelem, 0.0);
71 m_sa.assign(m_nelem, 0.0);
72 m_aw.assign(m_nsp, 0.0);
73 m_wx.assign(m_nelem, 0.0);
75 int solveFail =
false;
82 plogf(
"VCS CALCULATION METHOD\n\n ");
83 plogf(
"MultiPhase Object\n");
84 plogf(
"\n\n%5d SPECIES\n%5d ELEMENTS\n", m_nsp, m_nelem);
85 plogf(
"%5d COMPONENTS\n", m_numComponents);
86 plogf(
"%5d PHASES\n", m_numPhases);
87 plogf(
" PRESSURE%22.8g %3s\n", m_pressurePA,
"Pa ");
88 plogf(
" TEMPERATURE%19.3f K\n", m_temperature);
91 plogf(
" PHASE1 INERTS%17.3f\n", TPhInertMoles[0]);
93 if (m_numPhases > 1) {
94 plogf(
" PHASE2 INERTS%17.3f\n", TPhInertMoles[1]);
96 plogf(
"\n ELEMENTAL ABUNDANCES CORRECT");
97 plogf(
" FROM ESTIMATE Type\n\n");
98 for (
size_t i = 0; i < m_nelem; ++i) {
99 writeline(
' ', 26,
false);
100 plogf(
"%-2.2s", m_elementName[i]);
101 plogf(
"%20.12E%20.12E %3d\n", m_elemAbundancesGoal[i], m_elemAbundances[i],
104 if (m_doEstimateEquil < 0) {
105 plogf(
"\n MODIFIED LINEAR PROGRAMMING ESTIMATE OF EQUILIBRIUM - forced\n");
106 }
else if (m_doEstimateEquil > 0) {
107 plogf(
"\n MODIFIED LINEAR PROGRAMMING ESTIMATE OF EQUILIBRIUM - where necessary\n");
109 if (m_doEstimateEquil == 0) {
110 plogf(
"\n USER ESTIMATE OF EQUILIBRIUM\n");
112 plogf(
" Stan. Chem. Pot. in J/kmol\n");
113 plogf(
"\n SPECIES FORMULA VECTOR ");
114 writeline(
' ', 41,
false);
115 plogf(
" STAN_CHEM_POT EQUILIBRIUM_EST. Species_Type\n\n");
116 writeline(
' ', 20,
false);
117 for (
size_t i = 0; i < m_nelem; ++i) {
118 plogf(
"%-4.4s ", m_elementName[i]);
122 for (
size_t i = 0; i < m_nsp; ++i) {
123 plogf(
" %-18.18s", m_speciesName[i]);
124 for (
size_t j = 0; j < m_nelem; ++j) {
125 plogf(
"% -7.3g ", m_formulaMatrix(i,j));
127 plogf(
" %3d ", m_phaseID[i]);
128 writeline(
' ', std::max(55-
int(m_nelem)*8, 0),
false);
129 plogf(
"%12.5E %12.5E", RT * m_SSfeSpecies[i], m_molNumSpecies_old[i]);
140 for (
size_t i = 0; i < m_nsp; ++i) {
141 if (m_molNumSpecies_old[i] < 0.0) {
142 plogf(
"On Input species %-12s has a negative MF, setting it small\n",
144 size_t iph = m_phaseID[i];
147 m_molNumSpecies_old[i] = tmp;
163 if (forceComponentCalc) {
164 int retn = solve_tp_component_calc(allMinorZeroedSpecies);
169 forceComponentCalc = 0;
175 if (m_VCount->Its > maxit) {
178 solve_tp_inner(iti, it1, uptodate_minors, allMinorZeroedSpecies,
179 forceComponentCalc, stage, printDetails > 0, ANOTE);
181 }
else if (stage == EQUILIB_CHECK) {
183 solve_tp_equilib_check(allMinorZeroedSpecies, uptodate_minors,
184 giveUpOnElemAbund, solveFail, iti, it1,
186 }
else if (stage == ELEM_ABUND_CHECK) {
188 solve_tp_elem_abund_check(iti, stage, lec, giveUpOnElemAbund,
189 finalElemAbundAttempts, rangeErrorFound);
190 }
else if (stage == RECHECK_DELETED) {
200 size_t npb = vcs_recheck_deleted();
214 }
else if (stage == RETURN_A) {
216 size_t npb = vcs_recheck_deleted();
230 }
else if (stage == RETURN_B) {
233 size_t npb = vcs_add_all_deleted();
236 if (m_debug_print_lvl >= 1) {
237 plogf(
" --- add_all_deleted(): some rxns not converged. RETURNING TO LOOP!\n");
252 m_molNumSpecies_new.assign(m_molNumSpecies_new.size(), 0.0);
253 for (
size_t kspec = 0; kspec < m_nsp; ++kspec) {
254 if (m_SSPhase[kspec]) {
255 m_molNumSpecies_new[kspec] = 1.0;
257 size_t iph = m_phaseID[kspec];
258 if (m_tPhaseMoles_old[iph] != 0.0) {
259 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] / m_tPhaseMoles_old[iph];
265 size_t i = m_speciesLocalPhaseIndex[kspec];
266 m_molNumSpecies_new[kspec] = m_VolPhaseList[iph]->molefraction(i);
271 if (rangeErrorFound) {
277 m_VCount->Time_vcs_TP = tsecond;
278 m_VCount->T_Time_vcs_TP += m_VCount->Time_vcs_TP;
279 m_VCount->T_Calls_vcs_TP++;
280 m_VCount->T_Its += m_VCount->Its;
281 m_VCount->T_Basis_Opts += m_VCount->Basis_Opts;
282 m_VCount->T_Time_basopt += m_VCount->Time_basopt;
288 int VCS_SOLVE::solve_tp_component_calc(
bool& allMinorZeroedSpecies)
290 double test = -1.0e-10;
291 bool usedZeroedSpecies;
292 int retn = vcs_basopt(
false, &m_aw[0], &m_sa[0], &m_sm[0], &m_ss[0],
293 test, &usedZeroedSpecies);
303 allMinorZeroedSpecies = vcs_evaluate_speciesType();
306 if (! vcs_elabcheck(0)) {
307 debuglog(
" --- Element Abundance check failed\n", m_debug_print_lvl >= 2);
308 vcs_elcorr(&m_sm[0], &m_wx[0]);
315 debuglog(
" --- Element Abundance check passed\n", m_debug_print_lvl >= 2);
320 void VCS_SOLVE::solve_tp_inner(
size_t& iti,
size_t& it1,
321 bool& uptodate_minors,
322 bool& allMinorZeroedSpecies,
323 int& forceComponentCalc,
324 int& stage,
bool printDetails,
char* ANOTE)
333 if (!uptodate_minors) {
338 uptodate_minors =
true;
340 uptodate_minors =
false;
346 plogf(
" Iteration = %3d, Iterations since last evaluation of " 347 "optimal basis = %3d",
348 m_VCount->Its, it1 - 1);
350 plogf(
" (all species)\n");
352 plogf(
" (only major species)\n");
363 m_feSpecies_new = m_feSpecies_old;
364 m_actCoeffSpecies_new = m_actCoeffSpecies_old;
365 m_deltaGRxn_new = m_deltaGRxn_old;
366 m_deltaGRxn_Deficient = m_deltaGRxn_old;
367 m_tPhaseMoles_new = m_tPhaseMoles_old;
372 m_deltaMolNumSpecies.assign(m_deltaMolNumSpecies.size(), 0.0);
378 std::vector<size_t> phasePopPhaseIDs(0);
379 size_t iphasePop = vcs_popPhaseID(phasePopPhaseIDs);
380 if (iphasePop !=
npos) {
381 int soldel = vcs_popPhaseRxnStepSizes(iphasePop);
384 if (m_debug_print_lvl >= 2) {
385 plogf(
" --- vcs_popPhaseRxnStepSizes() was called but stoich " 386 "prevented phase %d popping\n");
394 size_t iphaseDelete =
npos;
396 if (iphasePop ==
npos) {
400 iphaseDelete = vcs_RxnStepSizes(forceComponentCalc, kspec);
401 }
else if (m_debug_print_lvl >= 2) {
402 plogf(
" --- vcs_RxnStepSizes not called because alternative" 403 "phase creation delta was used instead\n");
405 size_t doPhaseDeleteKspec =
npos;
406 size_t doPhaseDeleteIph =
npos;
409 m_deltaPhaseMoles.assign(m_deltaPhaseMoles.size(), 0.0);
425 if (iphaseDelete !=
npos) {
426 debuglog(
" --- Main Loop Treatment -> Circumvented due to Phase Deletion\n", m_debug_print_lvl >= 2);
427 for (
size_t k = 0; k < m_nsp; k++) {
428 m_molNumSpecies_new[k] = m_molNumSpecies_old[k] + m_deltaMolNumSpecies[k];
429 size_t iph = m_phaseID[k];
430 m_tPhaseMoles_new[iph] += m_deltaMolNumSpecies[k];
432 if (kspec >= m_numComponents) {
433 if (m_SSPhase[kspec] == 1) {
436 throw CanteraError(
"VCS_SOLVE::solve_tp_inner",
437 "we shouldn't be here!");
439 ++m_numRxnMinorZeroed;
440 allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc);
456 if (m_debug_print_lvl >= 2) {
457 plogf(
" --- Main Loop Treatment of each non-component species ");
459 plogf(
"- Full Calculation:\n");
461 plogf(
"- Major Components Calculation:\n");
463 plogf(
" --- Species IC ");
464 plogf(
" KMoles Tent_KMoles Rxn_Adj | Comment \n");
466 for (
size_t irxn = 0; irxn < m_numRxnRdc; irxn++) {
467 size_t kspec = m_indexRxnToSpecies[irxn];
468 double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
469 size_t iph = m_phaseID[kspec];
470 vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
474 if (iphasePop !=
npos) {
475 if (iph == iphasePop) {
476 dx = m_deltaMolNumSpecies[kspec];
477 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec];
478 sprintf(ANOTE,
"Phase pop");
481 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
486 dx = vcs_minor_alt_calc(kspec, irxn, &soldel_ret, ANOTE);
487 m_deltaMolNumSpecies[kspec] = dx;
490 bool resurrect = (m_deltaMolNumSpecies[kspec] > 0.0);
491 if (m_debug_print_lvl >= 3) {
492 plogf(
" --- %s currently zeroed (SpStatus=%-2d):",
493 m_speciesName[kspec], m_speciesStatus[kspec]);
494 plogf(
"%3d DG = %11.4E WT = %11.4E W = %11.4E DS = %11.4E\n",
495 irxn, m_deltaGRxn_new[irxn], m_molNumSpecies_new[kspec],
496 m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec]);
498 if (m_deltaGRxn_new[irxn] >= 0.0 || !resurrect) {
499 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
500 m_deltaMolNumSpecies[kspec] = 0.0;
502 sprintf(ANOTE,
"Species stays zeroed: DG = %11.4E", m_deltaGRxn_new[irxn]);
503 if (m_deltaGRxn_new[irxn] < 0.0) {
505 sprintf(ANOTE,
"Species stays zeroed even though dg neg due to " 506 "STOICH/PHASEPOP constraint: DG = %11.4E",
507 m_deltaGRxn_new[irxn]);
509 sprintf(ANOTE,
"Species stays zeroed even though dg neg: DG = %11.4E, ds zeroed",
510 m_deltaGRxn_new[irxn]);
514 for (
size_t j = 0; j < m_nelem; ++j) {
515 int elType = m_elType[j];
517 double atomComp = m_formulaMatrix(kspec,j);
518 if (atomComp > 0.0) {
519 double maxPermissible = m_elemAbundancesGoal[j] / atomComp;
521 sprintf(ANOTE,
"Species stays zeroed even though dG " 522 "neg, because of %s elemAbund",
523 m_elementName[j].c_str());
535 if (m_debug_print_lvl >= 2) {
536 plogf(
" --- Zeroed species changed to major: ");
537 plogf(
"%-12s\n", m_speciesName[kspec]);
540 allMinorZeroedSpecies =
false;
542 if (m_debug_print_lvl >= 2) {
543 plogf(
" --- Zeroed species changed to minor: ");
544 plogf(
"%-12s\n", m_speciesName[kspec]);
548 if (m_deltaMolNumSpecies[kspec] > 0.0) {
549 dx = m_deltaMolNumSpecies[kspec] * 0.01;
550 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
553 dx = m_molNumSpecies_new[kspec] - m_molNumSpecies_old[kspec];
555 m_deltaMolNumSpecies[kspec] = dx;
556 sprintf(ANOTE,
"Born:IC=-1 to IC=1:DG=%11.4E", m_deltaGRxn_new[irxn]);
558 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
559 m_deltaMolNumSpecies[kspec] = 0.0;
568 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
569 m_deltaMolNumSpecies[kspec] = 0.0;
571 sprintf(ANOTE,
"minor species not considered");
572 if (m_debug_print_lvl >= 2) {
573 plogf(
" --- %-12s%3d% 11.4E %11.4E %11.4E | %s\n",
574 m_speciesName[kspec], m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
575 m_deltaMolNumSpecies[kspec], ANOTE);
593 dx = vcs_minor_alt_calc(kspec, irxn, &soldel_ret, ANOTE);
594 m_deltaMolNumSpecies[kspec] = dx;
595 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
599 if (m_debug_print_lvl >= 2) {
600 plogf(
" --- Delete minor species in multispec phase: %-12s\n",
601 m_speciesName[kspec]);
603 m_deltaMolNumSpecies[kspec] = 0.0;
609 size_t lnospec = vcs_delete_species(kspec);
611 stage = RECHECK_DELETED;
625 sprintf(ANOTE,
"Normal Major Calc");
630 if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) {
631 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
632 m_deltaMolNumSpecies[kspec] = 0.0;
634 sprintf(ANOTE,
"major species is converged");
635 if (m_debug_print_lvl >= 2) {
636 plogf(
" --- %-12s %3d %11.4E %11.4E %11.4E | %s\n",
637 m_speciesName[kspec], m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
638 m_deltaMolNumSpecies[kspec], ANOTE);
650 if ((m_deltaGRxn_new[irxn] * m_deltaMolNumSpecies[kspec]) <= 0.0) {
651 dx = m_deltaMolNumSpecies[kspec];
654 m_deltaMolNumSpecies[kspec] = 0.0;
655 sprintf(ANOTE,
"dx set to 0, DG flipped sign due to " 656 "changed initial point");
660 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
666 if (m_molNumSpecies_new[kspec] <= 0.0) {
667 sprintf(ANOTE,
"initial nonpos kmoles= %11.3E",
668 m_molNumSpecies_new[kspec]);
676 if (!m_SSPhase[kspec]) {
681 dx = -0.9 * m_molNumSpecies_old[kspec];
682 m_deltaMolNumSpecies[kspec] = dx;
683 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
688 dx = -m_molNumSpecies_old[kspec];
694 for (
size_t j = 0; j < m_numComponents; ++j) {
695 if (sc_irxn[j] != 0.0) {
696 m_wx[j] = m_molNumSpecies_old[j] + sc_irxn[j] * dx;
697 if (m_wx[j] <= m_molNumSpecies_old[j] * 0.01 - 1.0E-150) {
698 dx = std::max(dx, m_molNumSpecies_old[j] * -0.99 / sc_irxn[j]);
701 m_wx[j] = m_molNumSpecies_old[j];
704 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
705 if (m_molNumSpecies_new[kspec] > 0.0) {
706 m_deltaMolNumSpecies[kspec] = dx;
708 "zeroing SS phase created a neg component species " 709 "-> reducing step size instead");
713 iph = m_phaseID[kspec];
714 Vphase = m_VolPhaseList[iph].get();
715 sprintf(ANOTE,
"zeroing out SS phase: ");
721 m_molNumSpecies_new[kspec] = 0.0;
722 doPhaseDeleteIph = iph;
724 doPhaseDeleteKspec = kspec;
725 if (m_debug_print_lvl >= 2 && m_speciesStatus[kspec] >= 0) {
726 plogf(
" --- SS species changed to zeroedss: ");
727 plogf(
"%-12s\n", m_speciesName[kspec]);
730 ++m_numRxnMinorZeroed;
731 allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc);
733 for (
size_t kk = 0; kk < m_nsp; kk++) {
734 m_deltaMolNumSpecies[kk] = 0.0;
735 m_molNumSpecies_new[kk] = m_molNumSpecies_old[kk];
737 m_deltaMolNumSpecies[kspec] = dx;
738 m_molNumSpecies_new[kspec] = 0.0;
740 for (
size_t k = 0; k < m_numComponents; ++k) {
741 m_deltaMolNumSpecies[k] = 0.0;
743 for (iph = 0; iph < m_numPhases; iph++) {
744 m_deltaPhaseMoles[iph] = 0.0;
752 if (dx != 0.0 && (m_speciesUnknownType[kspec] !=
759 1.0E-14*(fabs(m_deltaMolNumSpecies[kspec]) + fabs(dx) + 1.0E-32),
760 "VCS_SOLVE::solve_tp_inner",
761 "ds[kspec] = {}, dx = {}, kspec = {}\nwe have a problem!",
762 m_deltaMolNumSpecies[kspec], dx, kspec);
763 for (
size_t k = 0; k < m_numComponents; ++k) {
764 m_deltaMolNumSpecies[k] += sc_irxn[k] * dx;
769 for (iph = 0; iph < m_numPhases; iph++) {
770 m_deltaPhaseMoles[iph] += dx * m_deltaMolNumPhase(iph,irxn);
774 checkDelta1(&m_deltaMolNumSpecies[0],
775 &m_deltaPhaseMoles[0], kspec+1);
778 if (m_debug_print_lvl >= 2) {
779 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec];
780 plogf(
" --- %-12.12s%3d %11.4E %11.4E %11.4E | %s\n",
781 m_speciesName[kspec], m_speciesStatus[kspec],
782 m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
783 m_deltaMolNumSpecies[kspec], ANOTE);
786 if (doPhaseDeleteIph !=
npos) {
787 if (m_debug_print_lvl >= 2) {
788 plogf(
" --- %-12.12s Main Loop Special Case deleting phase with species:\n",
789 m_speciesName[doPhaseDeleteKspec]);
794 if (m_debug_print_lvl >= 2) {
795 for (
size_t k = 0; k < m_numComponents; k++) {
797 plogf(
"%-12.12s", m_speciesName[k]);
798 plogf(
" c %11.4E %11.4E %11.4E |\n",
799 m_molNumSpecies_old[k],
800 m_molNumSpecies_old[k]+m_deltaMolNumSpecies[k], m_deltaMolNumSpecies[k]);
804 plogf(
" --- Finished Main Loop\n");
813 for (
size_t k = 0; k < m_numComponents; ++k) {
814 if (m_molNumSpecies_old[k] > 0.0) {
815 double xx = -m_deltaMolNumSpecies[k] / m_molNumSpecies_old[k];
820 }
else if (m_deltaMolNumSpecies[k] < 0.0) {
823 size_t iph = m_phaseID[k];
824 m_deltaPhaseMoles[iph] -= m_deltaMolNumSpecies[k];
825 m_deltaMolNumSpecies[k] = 0.0;
829 if (par <= 1.01 && par > 0.0) {
832 if (m_debug_print_lvl >= 2) {
833 plogf(
" --- Reduction in step size due to component %s going negative = %11.3E\n",
834 m_speciesName[ll], par);
836 for (
size_t i = 0; i < m_nsp; ++i) {
837 m_deltaMolNumSpecies[i] *= par;
839 for (
size_t iph = 0; iph < m_numPhases; iph++) {
840 m_deltaPhaseMoles[iph] *= par;
845 checkDelta1(&m_deltaMolNumSpecies[0],
846 &m_deltaPhaseMoles[0], m_nsp);
853 for (
size_t kspec = 0; kspec < m_nsp; ++kspec) {
854 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec];
855 if (m_molNumSpecies_new[kspec] < 0.0 && (m_speciesUnknownType[kspec]
857 throw CanteraError(
"VCS_SOLVE::solve_tp_inner",
858 "vcs_solve_TP: ERROR on step change wt[{}:{}]: {} < 0.0",
859 kspec, m_speciesName[kspec], m_molNumSpecies_new[kspec]);
864 for (
size_t iph = 0; iph < m_numPhases; iph++) {
865 m_tPhaseMoles_new[iph] = m_tPhaseMoles_old[iph] + m_deltaPhaseMoles[iph];
883 plogf(
" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n",
884 vcs_Total_Gibbs(&m_molNumSpecies_old[0], &m_feSpecies_old[0],
885 &m_tPhaseMoles_old[0]));
886 plogf(
" --- Total tentative Dimensionless Gibbs Free Energy = %20.13E\n",
887 vcs_Total_Gibbs(&m_molNumSpecies_new[0], &m_feSpecies_new[0],
888 &m_tPhaseMoles_new[0]));
891 bool forced = vcs_globStepDamp();
894 if (printDetails && forced) {
895 plogf(
" -----------------------------------------------------\n");
896 plogf(
" --- FORCER SUBROUTINE changed the solution:\n");
897 plogf(
" --- SPECIES Status INIT MOLES TENT_MOLES");
898 plogf(
" FINAL KMOLES INIT_DEL_G/RT TENT_DEL_G/RT FINAL_DELTA_G/RT\n");
899 for (
size_t i = 0; i < m_numComponents; ++i) {
900 plogf(
" --- %-12.12s", m_speciesName[i]);
901 plogf(
" %14.6E %14.6E %14.6E\n", m_molNumSpecies_old[i],
902 m_molNumSpecies_old[i] + m_deltaMolNumSpecies[i], m_molNumSpecies_new[i]);
904 for (
size_t kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) {
905 size_t irxn = kspec - m_numComponents;
906 plogf(
" --- %-12.12s", m_speciesName[kspec]);
907 plogf(
" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", m_speciesStatus[kspec],
908 m_molNumSpecies_old[kspec],
909 m_molNumSpecies_old[kspec]+m_deltaMolNumSpecies[kspec],
910 m_molNumSpecies_new[kspec], m_deltaGRxn_old[irxn],
911 m_deltaGRxn_tmp[irxn], m_deltaGRxn_new[irxn]);
913 writeline(
' ', 26,
false);
914 plogf(
"Norms of Delta G():%14.6E%14.6E\n",
915 l2normdg(&m_deltaGRxn_old[0]),
916 l2normdg(&m_deltaGRxn_new[0]));
917 plogf(
" Total kmoles of gas = %15.7E\n", m_tPhaseMoles_old[0]);
918 if ((m_numPhases > 1) && (!m_VolPhaseList[1]->m_singleSpecies)) {
919 plogf(
" Total kmoles of liquid = %15.7E\n", m_tPhaseMoles_old[1]);
921 plogf(
" Total kmoles of liquid = %15.7E\n", 0.0);
923 plogf(
" Total New Dimensionless Gibbs Free Energy = %20.13E\n",
924 vcs_Total_Gibbs(&m_molNumSpecies_new[0], &m_feSpecies_new[0],
925 &m_tPhaseMoles_new[0]));
926 plogf(
" -----------------------------------------------------\n");
935 plogf(
" --- Summary of the Update ");
937 plogf(
" (all species):");
939 plogf(
" (only major species):");
942 plogf(
" --- Species Status Initial_KMoles Final_KMoles Initial_Mu/RT");
943 plogf(
" Mu/RT Init_Del_G/RT Delta_G/RT\n");
944 for (
size_t i = 0; i < m_numComponents; ++i) {
945 plogf(
" --- %-12.12s", m_speciesName[i]);
947 plogf(
"%14.6E%14.6E%14.6E%14.6E\n", m_molNumSpecies_old[i],
948 m_molNumSpecies_new[i], m_feSpecies_old[i], m_feSpecies_new[i]);
950 for (
size_t i = m_numComponents; i < m_numSpeciesRdc; ++i) {
951 size_t l1 = i - m_numComponents;
952 plogf(
" --- %-12.12s", m_speciesName[i]);
953 plogf(
" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n",
954 m_speciesStatus[i], m_molNumSpecies_old[i],
955 m_molNumSpecies_new[i], m_feSpecies_old[i], m_feSpecies_new[i],
956 m_deltaGRxn_old[l1], m_deltaGRxn_new[l1]);
958 for (
size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) {
959 size_t l1 = kspec - m_numComponents;
960 plogf(
" --- %-12.12s", m_speciesName[kspec]);
961 plogf(
" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n",
962 m_speciesStatus[kspec], m_molNumSpecies_old[kspec],
963 m_molNumSpecies_new[kspec], m_feSpecies_old[kspec], m_feSpecies_new[kspec],
964 m_deltaGRxn_old[l1], m_deltaGRxn_new[l1]);
967 writeline(
' ', 56,
false);
968 plogf(
"Norms of Delta G():%14.6E%14.6E\n",
969 l2normdg(&m_deltaGRxn_old[0]),
970 l2normdg(&m_deltaGRxn_new[0]));
972 plogf(
" --- Phase_Name KMoles(after update)\n");
975 for (
size_t iph = 0; iph < m_numPhases; iph++) {
976 vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
977 plogf(
" --- %18s = %15.7E\n", Vphase->PhaseName, m_tPhaseMoles_new[iph]);
981 plogf(
" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n",
982 vcs_Total_Gibbs(&m_molNumSpecies_old[0], &m_feSpecies_old[0],
983 &m_tPhaseMoles_old[0]));
984 plogf(
" --- Total New Dimensionless Gibbs Free Energy = %20.13E\n",
985 vcs_Total_Gibbs(&m_molNumSpecies_new[0], &m_feSpecies_new[0],
986 &m_tPhaseMoles_new[0]));
987 debuglog(
" --- Troublesome solve\n", m_VCount->Its > 550);
1000 m_tPhaseMoles_old = m_tPhaseMoles_new;
1001 m_molNumSpecies_old = m_molNumSpecies_new;
1002 m_actCoeffSpecies_old = m_actCoeffSpecies_new;
1003 m_deltaGRxn_old = m_deltaGRxn_new;
1004 m_feSpecies_old = m_feSpecies_new;
1011 if (m_debug_print_lvl >= 2) {
1012 plogf(
" --- Increment counter increased, step is accepted: %4d\n",
1022 bool justDeletedMultiPhase =
false;
1023 for (
size_t iph = 0; iph < m_numPhases; iph++) {
1024 if (!m_VolPhaseList[iph]->m_singleSpecies && m_tPhaseMoles_old[iph] != 0.0 &&
1026 if (m_debug_print_lvl >= 1) {
1027 plogf(
" --- Setting microscopic phase %d to zero\n", iph);
1029 justDeletedMultiPhase =
true;
1030 vcs_delete_multiphase(iph);
1037 if (justDeletedMultiPhase) {
1038 bool usedZeroedSpecies;
1039 double test = -1.0e-10;
1040 int retn = vcs_basopt(
false, &m_aw[0], &m_sa[0], &m_sm[0], &m_ss[0],
1041 test, &usedZeroedSpecies);
1043 throw CanteraError(
"VCS_SOLVE::solve_tp_inner",
1044 "BASOPT returned with an error condition");
1054 if (m_debug_print_lvl >= 2) {
1055 plogf(
" --- Normal element abundance check");
1058 if (! vcs_elabcheck(0)) {
1059 debuglog(
" - failed -> redoing element abundances.\n", m_debug_print_lvl >= 2);
1060 vcs_elcorr(&m_sm[0], &m_wx[0]);
1064 uptodate_minors =
true;
1066 debuglog(
" - passed\n", m_debug_print_lvl >= 2);
1070 for (
size_t i = 0; i < m_numRxnRdc; ++i) {
1071 size_t k = m_indexRxnToSpecies[i];
1075 for (
size_t j = 0; j < m_numComponents; ++j) {
1076 bool doSwap =
false;
1078 doSwap = (m_molNumSpecies_old[k] * m_spSize[k]) >
1079 (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
1080 if (!m_SSPhase[k] && doSwap) {
1081 doSwap = m_molNumSpecies_old[k] > (m_molNumSpecies_old[j] * 1.01);
1085 doSwap = (m_molNumSpecies_old[k] * m_spSize[k]) >
1086 (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
1088 doSwap = m_molNumSpecies_old[k] > (m_molNumSpecies_old[j] * 1.01);
1091 doSwap = (m_molNumSpecies_old[k] * m_spSize[k]) >
1092 (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
1095 if (doSwap && m_stoichCoeffRxnMatrix(j,i) != 0.0) {
1096 if (m_debug_print_lvl >= 2) {
1097 plogf(
" --- Get a new basis because %s", m_speciesName[k]);
1098 plogf(
" is better than comp %s", m_speciesName[j]);
1099 plogf(
" and share nonzero stoic: %-9.1f\n",
1100 m_stoichCoeffRxnMatrix(j,i));
1102 forceComponentCalc = 1;
1107 debuglog(
" --- Check for an optimum basis passed\n", m_debug_print_lvl >= 2);
1108 stage = EQUILIB_CHECK;
1117 if (m_debug_print_lvl >= 2) {
1118 plogf(
" --- Reevaluate major-minor status of noncomponents:\n");
1120 m_numRxnMinorZeroed = 0;
1121 for (
size_t irxn = 0; irxn < m_numRxnRdc; irxn++) {
1122 size_t kspec = m_indexRxnToSpecies[irxn];
1123 int speciesType = vcs_species_type(kspec);
1126 plogf(
" --- major/minor species is now zeroed out: %s\n",
1127 m_speciesName[kspec]);
1129 ++m_numRxnMinorZeroed;
1133 plogf(
" --- Noncomponent turned from major to minor: ");
1134 }
else if (kspec < m_numComponents) {
1135 plogf(
" --- Component turned into a minor species: ");
1137 plogf(
" --- Zeroed Species turned into a " 1140 plogf(
"%s\n", m_speciesName[kspec]);
1142 ++m_numRxnMinorZeroed;
1145 if (m_debug_print_lvl >= 2) {
1147 plogf(
" --- Noncomponent turned from minor to major: ");
1148 }
else if (kspec < m_numComponents) {
1149 plogf(
" --- Component turned into a major: ");
1151 plogf(
" --- Noncomponent turned from zeroed to major: ");
1153 plogf(
"%s\n", m_speciesName[kspec]);
1158 m_speciesStatus[kspec] = speciesType;
1163 allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc);
1166 void VCS_SOLVE::solve_tp_equilib_check(
bool& allMinorZeroedSpecies,
1167 bool& uptodate_minors,
1168 bool& giveUpOnElemAbund,
int& solveFail,
1169 size_t& iti,
size_t& it1,
int maxit,
1170 int& stage,
bool& lec)
1172 if (! allMinorZeroedSpecies) {
1173 if (m_debug_print_lvl >= 2) {
1174 plogf(
" --- Equilibrium check for major species: ");
1176 for (
size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
1177 size_t kspec = irxn + m_numComponents;
1178 if (m_speciesStatus[kspec] ==
VCS_SPECIES_MAJOR && (fabs(m_deltaGRxn_new[irxn]) > m_tolmaj)) {
1179 if (m_VCount->Its >= maxit) {
1186 if (m_debug_print_lvl >= 2) {
1187 plogf(
"%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]]);
1191 iti = ((it1/4) *4) - it1;
1197 debuglog(
" MAJOR SPECIES CONVERGENCE achieved", m_debug_print_lvl >= 2);
1199 debuglog(
" MAJOR SPECIES CONVERGENCE achieved " 1200 "(because there are no major species)\n", m_debug_print_lvl >= 2);
1205 if (m_numRxnMinorZeroed != 0) {
1212 uptodate_minors =
true;
1214 if (m_debug_print_lvl >= 2) {
1215 plogf(
" --- Equilibrium check for minor species: ");
1217 for (
size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
1218 size_t kspec = irxn + m_numComponents;
1219 if (m_speciesStatus[kspec] ==
VCS_SPECIES_MINOR && (fabs(m_deltaGRxn_new[irxn]) > m_tolmin)) {
1220 if (m_VCount->Its >= maxit) {
1227 if (m_debug_print_lvl >= 2) {
1228 plogf(
"%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]]);
1238 if (m_debug_print_lvl >= 2) {
1239 plogf(
" CONVERGENCE achieved\n");
1250 if (!giveUpOnElemAbund) {
1251 if (m_debug_print_lvl >= 2) {
1252 plogf(
" --- Check the Full Element Abundances: ");
1257 if (! vcs_elabcheck(1)) {
1258 if (m_debug_print_lvl >= 2) {
1259 if (! vcs_elabcheck(0)) {
1262 plogf(
" passed for NC but failed for NE: RANGE ERROR\n");
1266 stage = ELEM_ABUND_CHECK;
1269 if (m_debug_print_lvl >= 2) {
1276 if (m_numSpeciesRdc != m_nsp) {
1277 stage = RECHECK_DELETED;
1286 void VCS_SOLVE::solve_tp_elem_abund_check(
size_t& iti,
int& stage,
bool& lec,
1287 bool& giveUpOnElemAbund,
1288 int& finalElemAbundAttempts,
1289 int& rangeErrorFound)
1296 rangeErrorFound = 0;
1297 if (! vcs_elabcheck(1)) {
1298 bool ncBefore = vcs_elabcheck(0);
1299 vcs_elcorr(&m_sm[0], &m_wx[0]);
1300 bool ncAfter = vcs_elabcheck(0);
1301 bool neAfter = vcs_elabcheck(1);
1317 if (finalElemAbundAttempts >= 3) {
1318 giveUpOnElemAbund =
true;
1319 stage = EQUILIB_CHECK;
1321 finalElemAbundAttempts++;
1328 }
else if (ncAfter) {
1331 if (m_debug_print_lvl >= 2) {
1332 plogf(
" --- vcs_solve_tp: RANGE SPACE ERROR ENCOUNTERED\n");
1333 plogf(
" --- vcs_solve_tp: - Giving up on NE Element Abundance satisfaction \n");
1334 plogf(
" --- vcs_solve_tp: - However, NC Element Abundance criteria is satisfied \n");
1335 plogf(
" --- vcs_solve_tp: - Returning the calculated equilibrium condition\n");
1337 rangeErrorFound = 1;
1338 giveUpOnElemAbund =
true;
1343 stage = EQUILIB_CHECK;
1350 stage = EQUILIB_CHECK;
1353 double VCS_SOLVE::vcs_minor_alt_calc(
size_t kspec,
size_t irxn,
bool* do_delete,
1356 double w_kspec = m_molNumSpecies_old[kspec];
1357 double dg_irxn = m_deltaGRxn_old[irxn];
1358 size_t iph = m_phaseID[kspec];
1362 if (w_kspec <= 0.0) {
1365 dg_irxn = std::max(dg_irxn, -200.0);
1367 sprintf(ANOTE,
"minor species alternative calc");
1369 if (dg_irxn >= 23.0) {
1370 double molNum_kspec_new = w_kspec * 1.0e-10;
1377 return molNum_kspec_new - w_kspec;
1379 if (fabs(dg_irxn) <= m_tolmin2) {
1385 double s = m_np_dLnActCoeffdMolNum(kspec,kspec) / m_tPhaseMoles_old[iph];
1396 double a =
clip(w_kspec * s, -1.0+1e-8, 100.0);
1397 double tmp =
clip(-dg_irxn / (1.0 + a), -200.0, 200.0);
1398 double wTrial = w_kspec * exp(tmp);
1399 double molNum_kspec_new = wTrial;
1401 if (wTrial > 100. * w_kspec) {
1402 double molNumMax = 0.0001 * m_tPhaseMoles_old[iph];
1403 if (molNumMax < 100. * w_kspec) {
1404 molNumMax = 100. * w_kspec;
1406 if (wTrial > molNumMax) {
1407 molNum_kspec_new = molNumMax;
1409 molNum_kspec_new = wTrial;
1411 }
else if (1.0E10 * wTrial < w_kspec) {
1412 molNum_kspec_new= 1.0E-10 * w_kspec;
1414 molNum_kspec_new = wTrial;
1423 return molNum_kspec_new - w_kspec;
1428 double dx = m_deltaGRxn_old[irxn]/ m_Faraday_dim;
1430 sprintf(ANOTE,
"voltage species alternative calc");
1436 int VCS_SOLVE::delta_species(
const size_t kspec,
double*
const delta_ptr)
1438 size_t irxn = kspec - m_numComponents;
1440 AssertThrowMsg(kspec >= m_numComponents,
"VCS_SOLVE::delta_species",
1441 "delete_species() ERROR: called for a component {}", kspec);
1445 double dx = *delta_ptr;
1446 double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
1447 for (
size_t j = 0; j < m_numComponents; ++j) {
1448 if (m_molNumSpecies_old[j] > 0.0) {
1449 double tmp = sc_irxn[j] * dx;
1450 if (-tmp > m_molNumSpecies_old[j]) {
1452 dx = std::min(dx, - m_molNumSpecies_old[j] / sc_irxn[j]);
1458 if (m_molNumSpecies_old[j] <= 0.0 && sc_irxn[j] < 0.0) {
1466 m_molNumSpecies_old[kspec] += dx;
1467 size_t iph = m_phaseID[kspec];
1468 m_tPhaseMoles_old[iph] += dx;
1471 for (
size_t j = 0; j < m_numComponents; ++j) {
1472 double tmp = sc_irxn[j] * dx;
1475 m_molNumSpecies_old[j] += tmp;
1476 m_tPhaseMoles_old[iph] += tmp;
1478 m_molNumSpecies_old[j] = std::max(m_molNumSpecies_old[j], 0.0);
1485 int VCS_SOLVE::vcs_zero_species(
const size_t kspec)
1491 double dx = -m_molNumSpecies_old[kspec];
1493 retn = delta_species(kspec, &dx);
1494 if (!retn && m_debug_print_lvl >= 1) {
1495 plogf(
"vcs_zero_species: Couldn't zero the species %d, " 1496 "did delta of %g. orig conc of %g\n",
1497 kspec, dx, m_molNumSpecies_old[kspec] + dx);
1504 int VCS_SOLVE::vcs_delete_species(
const size_t kspec)
1506 const size_t klast = m_numSpeciesRdc - 1;
1507 const size_t iph = m_phaseID[kspec];
1508 vcs_VolPhase*
const Vphase = m_VolPhaseList[iph].get();
1509 const size_t irxn = kspec - m_numComponents;
1513 const int retn = vcs_zero_species(kspec);
1516 "Failed to delete a species!");
1522 --m_numRxnMinorZeroed;
1525 m_deltaGRxn_new[irxn] = 0.0;
1526 m_deltaGRxn_old[irxn] = 0.0;
1527 m_feSpecies_new[kspec] = 0.0;
1528 m_feSpecies_old[kspec] = 0.0;
1529 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
1532 if (kspec != klast) {
1533 vcs_switch_pos(
true, klast, kspec);
1538 &m_tPhaseMoles_old[0]);
1547 bool stillExists =
false;
1548 for (
size_t k = 0; k < m_numSpeciesRdc; k++) {
1550 m_phaseID[k] == iph && m_molNumSpecies_old[k] > 0.0) {
1556 vcs_delete_multiphase(iph);
1562 return (m_numRxnRdc == 0);
1565 void VCS_SOLVE::vcs_reinsert_deleted(
size_t kspec)
1567 size_t iph = m_phaseID[kspec];
1568 if (m_debug_print_lvl >= 2) {
1569 plogf(
" --- Add back a deleted species: %-12s\n", m_speciesName[kspec]);
1576 delta_species(kspec, &dx);
1579 if (m_SSPhase[kspec]) {
1581 --m_numRxnMinorZeroed;
1586 &m_molNumSpecies_old[0],
1587 &m_tPhaseMoles_old[0]);
1593 if (! m_SSPhase[kspec]) {
1596 for (
size_t k = 0; k < m_nsp; k++) {
1608 ++m_numRxnMinorZeroed;
1610 if (kspec != (m_numSpeciesRdc - 1)) {
1612 vcs_switch_pos(
true, (m_numSpeciesRdc - 1), kspec);
1616 bool VCS_SOLVE::vcs_delete_multiphase(
const size_t iph)
1619 bool successful =
true;
1623 if (m_debug_print_lvl >= 2) {
1624 plogf(
" --- delete_multiphase %d, %s\n", iph, Vphase->
PhaseName);
1628 for (
size_t kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) {
1629 if (m_phaseID[kspec] == iph) {
1632 double dx = - m_molNumSpecies_old[kspec];
1635 int retn = delta_species(kspec, &dxTent);
1638 if (m_debug_print_lvl >= 2) {
1639 plogf(
" --- delete_multiphase %d, %s ERROR problems deleting species %s\n",
1640 iph, Vphase->
PhaseName, m_speciesName[kspec]);
1641 plogf(
" --- delta attempted: %g achieved: %g " 1642 " Zeroing it manually\n", dx, dxTent);
1644 m_molNumSpecies_old[kspec] = 0.0;
1645 m_molNumSpecies_new[kspec] = 0.0;
1646 m_deltaMolNumSpecies[kspec] = 0.0;
1651 m_molNumSpecies_old[kspec] = 0.0;
1652 m_molNumSpecies_new[kspec] = 0.0;
1653 m_deltaMolNumSpecies[kspec] = 0.0;
1663 double dxPerm = 0.0, dxPerm2 = 0.0;
1664 for (
size_t kcomp = 0; kcomp < m_numComponents; ++kcomp) {
1665 if (m_phaseID[kcomp] == iph) {
1666 if (m_debug_print_lvl >= 2) {
1667 plogf(
" --- delete_multiphase One of the species is a component %d - %s with mole number %g\n",
1668 kcomp, m_speciesName[kcomp], m_molNumSpecies_old[kcomp]);
1670 if (m_molNumSpecies_old[kcomp] != 0.0) {
1671 for (
size_t kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) {
1672 size_t irxn = kspec - m_numComponents;
1673 if (m_phaseID[kspec] != iph) {
1674 if (m_stoichCoeffRxnMatrix(kcomp,irxn) != 0.0) {
1675 double dxWant = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix(kcomp,irxn);
1676 if (dxWant + m_molNumSpecies_old[kspec] < 0.0) {
1677 dxPerm = -m_molNumSpecies_old[kspec];
1679 for (
size_t jcomp = 0; kcomp < m_numComponents; ++kcomp) {
1680 if (jcomp != kcomp) {
1681 if (m_phaseID[jcomp] == iph) {
1684 double dj = dxWant * m_stoichCoeffRxnMatrix(jcomp,irxn);
1685 if (dj + m_molNumSpecies_old[kcomp] < 0.0) {
1686 dxPerm2 = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix(jcomp,irxn);
1688 if (fabs(dxPerm2) < fabs(dxPerm)) {
1695 if (dxPerm != 0.0) {
1696 delta_species(kspec, &dxPerm);
1702 if (m_molNumSpecies_old[kcomp] != 0.0) {
1703 if (m_debug_print_lvl >= 2) {
1704 plogf(
" --- delete_multiphase One of the species is a component %d - %s still with mole number %g\n",
1705 kcomp, m_speciesName[kcomp], m_molNumSpecies_old[kcomp]);
1706 plogf(
" --- zeroing it \n");
1708 m_molNumSpecies_old[kcomp] = 0.0;
1719 for (
size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) {
1720 if (m_phaseID[kspec] == iph) {
1721 m_molNumSpecies_old[kspec] = 0.0;
1722 m_molNumSpecies_new[kspec] = 0.0;
1723 m_deltaMolNumSpecies[kspec] = 0.0;
1728 if (m_debug_print_lvl >= 2) {
1729 plogf(
" --- Make %s", m_speciesName[kspec]);
1730 plogf(
" an active but zeroed species because its phase " 1733 if (kspec != (m_numSpeciesRdc - 1)) {
1735 vcs_switch_pos(
true, (m_numSpeciesRdc - 1), kspec);
1741 m_tPhaseMoles_old[iph] = 0.0;
1742 m_tPhaseMoles_new[iph] = 0.0;
1743 m_deltaPhaseMoles[iph] = 0.0;
1750 int VCS_SOLVE::vcs_recheck_deleted()
1753 if (m_debug_print_lvl >= 2) {
1754 plogf(
" --- Start rechecking deleted species in multispec phases\n");
1756 if (m_numSpeciesRdc == m_nsp) {
1763 for (
size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) {
1764 size_t iph = m_phaseID[kspec];
1765 m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec])
1766 - m_lnMnaughtSpecies[kspec]
1767 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]);
1774 for (
size_t iph = 0; iph < m_numPhases; iph++) {
1775 if (m_tPhaseMoles_old[iph] > 0.0) {
1778 xtcutoff[iph] = 0.0;
1804 for (
size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
1805 size_t kspec = m_indexRxnToSpecies[irxn];
1806 size_t iph = m_phaseID[kspec];
1807 if (m_tPhaseMoles_old[iph] == 0.0) {
1808 if (m_deltaGRxn_new[irxn] < 0.0) {
1809 vcs_reinsert_deleted(kspec);
1812 m_molNumSpecies_old[kspec] = 0.0;
1814 }
else if (m_tPhaseMoles_old[iph] > 0.0) {
1815 if (m_deltaGRxn_new[irxn] < xtcutoff[iph]) {
1816 vcs_reinsert_deleted(kspec);
1824 size_t VCS_SOLVE::vcs_add_all_deleted()
1826 if (m_numSpeciesRdc == m_nsp) {
1835 m_molNumSpecies_new = m_molNumSpecies_old;
1836 for (
int cits = 0; cits < 3; cits++) {
1837 for (
size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) {
1838 size_t iph = m_phaseID[kspec];
1840 if (m_molNumSpecies_new[kspec] == 0.0) {
1846 m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_new[kspec]) - m_lnMnaughtSpecies[kspec]
1847 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]);
1853 for (
size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
1854 size_t kspec = m_indexRxnToSpecies[irxn];
1855 size_t iph = m_phaseID[kspec];
1856 if (m_tPhaseMoles_old[iph] > 0.0) {
1857 double maxDG = std::min(m_deltaGRxn_new[irxn], 690.0);
1858 double dx = m_tPhaseMoles_old[iph] * exp(- maxDG);
1859 m_molNumSpecies_new[kspec] = dx;
1863 for (
size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
1864 size_t kspec = m_indexRxnToSpecies[irxn];
1865 size_t iph = m_phaseID[kspec];
1866 if (m_tPhaseMoles_old[iph] > 0.0) {
1867 double dx = m_molNumSpecies_new[kspec];
1868 size_t retn = delta_species(kspec, &dx);
1870 if (m_debug_print_lvl) {
1871 plogf(
" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
1872 m_speciesName[kspec], kspec, dx);
1876 retn = delta_species(kspec, &dx);
1877 if (retn == 0 && m_debug_print_lvl) {
1878 plogf(
" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
1879 m_speciesName[kspec], kspec, dx);
1883 if (m_debug_print_lvl >= 2) {
1885 plogf(
" --- add_deleted(): species %s added back in with mol number %g\n",
1886 m_speciesName[kspec], dx);
1888 plogf(
" --- add_deleted(): species %s failed to be added back in\n");
1898 for (
size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
1899 size_t kspec = m_indexRxnToSpecies[irxn];
1900 size_t iph = m_phaseID[kspec];
1901 if (m_tPhaseMoles_old[iph] > 0.0 && fabs(m_deltaGRxn_old[irxn]) > m_tolmin) {
1902 if (((m_molNumSpecies_old[kspec] * exp(-m_deltaGRxn_old[irxn])) >
1906 if (m_debug_print_lvl >= 2) {
1907 plogf(
" --- add_deleted(): species %s with mol number %g not converged: DG = %g\n",
1908 m_speciesName[kspec], m_molNumSpecies_old[kspec],
1909 m_deltaGRxn_old[irxn]);
1917 bool VCS_SOLVE::vcs_globStepDamp()
1919 double* dptr = &m_deltaGRxn_new[0];
1923 for (
size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
1924 size_t kspec = irxn + m_numComponents;
1926 s2 += dptr[irxn] * m_deltaMolNumSpecies[kspec];
1932 dptr = &m_deltaGRxn_old[0];
1933 for (
size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
1934 size_t kspec = irxn + m_numComponents;
1936 s1 += dptr[irxn] * m_deltaMolNumSpecies[kspec];
1940 if (m_debug_print_lvl >= 2) {
1941 plogf(
" --- subroutine FORCE: Beginning Slope = %g\n", s1);
1942 plogf(
" --- subroutine FORCE: End Slope = %g\n", s2);
1946 if (m_debug_print_lvl >= 2) {
1947 plogf(
" --- subroutine FORCE produced no adjustments,");
1949 plogf(
" s1 positive but really small\n");
1951 plogf(
" failed s1 test\n");
1958 debuglog(
" --- subroutine FORCE produced no adjustments, s2 < 0\n", m_debug_print_lvl >= 2);
1964 if (fabs(s1 -s2) > 1.0E-200) {
1965 al = s1 / (s1 - s2);
1967 if (al >= 0.95 || al < 0.0) {
1968 if (m_debug_print_lvl >= 2) {
1969 plogf(
" --- subroutine FORCE produced no adjustments (al = %g)\n", al);
1973 if (m_debug_print_lvl >= 2) {
1974 plogf(
" --- subroutine FORCE produced a damping factor = %g\n", al);
1978 if (m_debug_print_lvl >= 2) {
1979 m_deltaGRxn_tmp = m_deltaGRxn_new;
1982 dptr = &m_molNumSpecies_new[0];
1983 for (
size_t kspec = 0; kspec < m_numSpeciesRdc; ++kspec) {
1984 m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] +
1985 al * m_deltaMolNumSpecies[kspec];
1987 for (
size_t iph = 0; iph < m_numPhases; iph++) {
1988 m_tPhaseMoles_new[iph] = m_tPhaseMoles_old[iph] + al * m_deltaPhaseMoles[iph];
1992 if (m_debug_print_lvl >= 2) {
1993 plogf(
" --- subroutine FORCE adjusted the mole " 1994 "numbers, AL = %10.3f\n", al);
2007 dptr = &m_deltaGRxn_new[0];
2009 for (
size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
2010 size_t kspec = irxn + m_numComponents;
2012 s2 += dptr[irxn] * m_deltaMolNumSpecies[kspec];
2016 if (m_debug_print_lvl >= 2) {
2017 plogf(
" --- subroutine FORCE: Adj End Slope = %g\n", s2);
2022 int VCS_SOLVE::vcs_basopt(
const bool doJustComponents,
double aw[],
double sa[],
double sm[],
2023 double ss[],
double test,
bool*
const usedZeroedSpecies)
2027 size_t jlose =
npos;
2030 if (m_debug_print_lvl >= 2) {
2032 for (
size_t i=0; i<77; i++) {
2036 plogf(
" --- Subroutine BASOPT called to ");
2037 if (doJustComponents) {
2038 plogf(
"calculate the number of components\n");
2040 plogf(
"reevaluate the components\n");
2042 if (m_debug_print_lvl >= 2) {
2044 plogf(
" --- Formula Matrix used in BASOPT calculation\n");
2045 plogf(
" --- Active | ");
2046 for (
size_t j = 0; j < m_nelem; j++) {
2047 plogf(
" %1d ", m_elementActive[j]);
2050 plogf(
" --- Species | ");
2051 for (
size_t j = 0; j < m_nelem; j++) {
2052 writelog(
" {:>8.8s}", m_elementName[j]);
2055 for (k = 0; k < m_nsp; k++) {
2056 writelog(
" --- {:>11.11s} | ", m_speciesName[k]);
2057 for (
size_t j = 0; j < m_nelem; j++) {
2058 plogf(
" %8.2g", m_formulaMatrix(k,j));
2069 size_t ncTrial = std::min(m_nelem, m_nsp);
2070 m_numComponents = ncTrial;
2071 *usedZeroedSpecies =
false;
2075 std::copy(m_molNumSpecies_old.begin(),
2076 m_molNumSpecies_old.begin() + m_nsp, aw);
2079 for (k = 0; k < m_nsp; k++) {
2089 while (jr < ncTrial) {
2097 k = vcs_basisOptMax(aw, jr, m_nsp);
2137 *usedZeroedSpecies =
true;
2138 double maxConcPossKspec = 0.0;
2139 double maxConcPoss = 0.0;
2140 size_t kfound =
npos;
2141 int minNonZeroes = 100000;
2142 int nonZeroesKspec = 0;
2143 for (
size_t kspec = ncTrial; kspec < m_nsp; kspec++) {
2145 maxConcPossKspec = 1.0E10;
2147 for (
size_t j = 0; j < m_nelem; ++j) {
2149 double nu = m_formulaMatrix(kspec,j);
2152 maxConcPossKspec = std::min(m_elemAbundancesGoal[j] / nu, maxConcPossKspec);
2156 if ((maxConcPossKspec >= maxConcPoss) || (maxConcPossKspec > 1.0E-5)) {
2157 if (nonZeroesKspec <= minNonZeroes) {
2158 if (kfound ==
npos || nonZeroesKspec < minNonZeroes) {
2163 if (m_SSfeSpecies[kspec] <= m_SSfeSpecies[kfound]) {
2168 minNonZeroes = std::min(minNonZeroes, nonZeroesKspec);
2169 maxConcPoss = std::max(maxConcPoss, maxConcPossKspec);
2173 if (kfound ==
npos) {
2176 for (
size_t kspec = ncTrial; kspec < m_nsp; kspec++) {
2177 if (aw[kspec] >= 0.0) {
2178 size_t irxn = kspec - ncTrial;
2179 if (m_deltaGRxn_new[irxn] < gmin) {
2180 gmin = m_deltaGRxn_new[irxn];
2189 if (aw[k] == test) {
2190 m_numComponents = jr;
2191 ncTrial = m_numComponents;
2192 size_t numPreDeleted = m_numRxnTot - m_numRxnRdc;
2193 if (numPreDeleted != (m_nsp - m_numSpeciesRdc)) {
2194 throw CanteraError(
"VCS_SOLVE::vcs_basopt",
"we shouldn't be here");
2196 m_numRxnTot = m_nsp - ncTrial;
2197 m_numRxnRdc = m_numRxnTot - numPreDeleted;
2198 m_numSpeciesRdc = m_nsp - numPreDeleted;
2199 for (
size_t i = 0; i < m_nsp; ++i) {
2200 m_indexRxnToSpecies[i] = ncTrial + i;
2202 if (m_debug_print_lvl >= 2) {
2203 plogf(
" --- Total number of components found = %3d (ne = %d)\n ",
2218 for (
size_t j = 0; j < m_nelem; ++j) {
2219 sm[j + jr*m_nelem] = m_formulaMatrix(k,j);
2225 for (
size_t j = 0; j < jl; ++j) {
2227 for (
size_t i = 0; i < m_nelem; ++i) {
2228 ss[j] += sm[i + jr*m_nelem] * sm[i + j*m_nelem];
2234 for (
size_t j = 0; j < jl; ++j) {
2235 for (
size_t i = 0; i < m_nelem; ++i) {
2236 sm[i + jr*m_nelem] -= ss[j] * sm[i + j*m_nelem];
2244 for (
size_t ml = 0; ml < m_nelem; ++ml) {
2245 sa[jr] += pow(sm[ml + jr*m_nelem], 2);
2249 if (sa[jr] > 1.0e-6) {
2256 if (m_debug_print_lvl >= 2) {
2257 plogf(
" --- %-12.12s", m_speciesName[k]);
2259 plogf(
"(Volts = %9.2g)", m_molNumSpecies_old[k]);
2261 plogf(
"(%9.2g)", m_molNumSpecies_old[k]);
2263 plogf(
" replaces %-12.12s", m_speciesName[jr]);
2265 plogf(
"(Volts = %9.2g)", m_molNumSpecies_old[jr]);
2267 plogf(
"(%9.2g)", m_molNumSpecies_old[jr]);
2269 plogf(
" as component %3d\n", jr);
2271 vcs_switch_pos(
false, jr, k);
2272 std::swap(aw[jr], aw[k]);
2273 }
else if (m_debug_print_lvl >= 2) {
2274 plogf(
" --- %-12.12s", m_speciesName[k]);
2276 plogf(
"(Volts = %9.2g) remains ", m_molNumSpecies_old[k]);
2278 plogf(
"(%9.2g) remains ", m_molNumSpecies_old[k]);
2280 plogf(
" as component %3d\n", jr);
2289 if (doJustComponents) {
2319 C.
resize(ncTrial, ncTrial);
2320 for (
size_t j = 0; j < ncTrial; ++j) {
2321 for (
size_t i = 0; i < ncTrial; ++i) {
2322 C(i, j) = m_formulaMatrix(j,i);
2325 for (
size_t i = 0; i < m_numRxnTot; ++i) {
2326 k = m_indexRxnToSpecies[i];
2327 for (
size_t j = 0; j < ncTrial; ++j) {
2328 m_stoichCoeffRxnMatrix(j,i) = - m_formulaMatrix(k,j);
2333 solve(C, m_stoichCoeffRxnMatrix.
ptrColumn(0), m_numRxnTot, m_nelem);
2339 for (
size_t j = 0; j < m_nelem; j++) {
2340 if (!m_elementActive[j] && !strcmp(m_elementName[j].c_str(),
"E")) {
2344 for (
size_t j = 0; j < m_nelem; j++) {
2345 if (m_elementActive[j] && !strncmp((m_elementName[j]).c_str(),
"cn_", 3)) {
2349 for (k = 0; k < m_nsp; k++) {
2351 for (
size_t j = 0; j < ncTrial; ++j) {
2352 for (
size_t i = 0; i < ncTrial; ++i) {
2354 C(i, j) = m_formulaMatrix(j,juse);
2356 C(i, j) = m_formulaMatrix(j,i);
2360 for (
size_t i = 0; i < m_numRxnTot; ++i) {
2361 k = m_indexRxnToSpecies[i];
2362 for (
size_t j = 0; j < ncTrial; ++j) {
2364 aw[j] = - m_formulaMatrix(k,juse);
2366 aw[j] = - m_formulaMatrix(k,j);
2371 solve(C, aw, 1, m_nelem);
2372 size_t i = k - ncTrial;
2373 for (
size_t j = 0; j < ncTrial; j++) {
2374 m_stoichCoeffRxnMatrix(j,i) = aw[j];
2380 for (
size_t i = 0; i < m_numRxnTot; i++) {
2382 for (
size_t j = 0; j < ncTrial; j++) {
2383 szTmp += fabs(m_stoichCoeffRxnMatrix(j,i));
2385 m_scSize[i] = szTmp;
2388 if (m_debug_print_lvl >= 2) {
2389 plogf(
" --- Components:");
2390 for (
size_t j = 0; j < ncTrial; j++) {
2393 plogf(
"\n --- Components Moles:");
2394 for (
size_t j = 0; j < ncTrial; j++) {
2396 plogf(
" % -10.3E", 0.0);
2398 plogf(
" % -10.3E", m_molNumSpecies_old[j]);
2401 plogf(
"\n --- NonComponent| Moles |");
2402 for (
size_t j = 0; j < ncTrial; j++) {
2403 plogf(
" %10.10s", m_speciesName[j]);
2406 for (
size_t i = 0; i < m_numRxnTot; i++) {
2407 plogf(
" --- %3d ", m_indexRxnToSpecies[i]);
2408 plogf(
"%-10.10s", m_speciesName[m_indexRxnToSpecies[i]]);
2410 plogf(
"|% -10.3E|", 0.0);
2412 plogf(
"|% -10.3E|", m_molNumSpecies_old[m_indexRxnToSpecies[i]]);
2414 for (
size_t j = 0; j < ncTrial; j++) {
2415 plogf(
" %+7.3f", m_stoichCoeffRxnMatrix(j,i));
2422 double sumMax = -1.0;
2425 for (
size_t i = 0; i < m_numRxnTot; ++i) {
2426 k = m_indexRxnToSpecies[i];
2428 for (
size_t j = 0; j < ncTrial; ++j) {
2430 sum = m_formulaMatrix(k,juse);
2431 for (
size_t n = 0; n < ncTrial; n++) {
2432 double numElements = m_formulaMatrix(n,juse);
2433 double coeff = m_stoichCoeffRxnMatrix(n,i);
2434 sum += coeff * numElements;
2437 sum = m_formulaMatrix(k,j);
2438 for (
size_t n = 0; n < ncTrial; n++) {
2439 double numElements = m_formulaMatrix(n,j);
2440 double coeff = m_stoichCoeffRxnMatrix(n,i);
2441 sum += coeff * numElements;
2444 if (fabs(sum) > sumMax) {
2452 if (fabs(sum) > 1.0E-6) {
2453 throw CanteraError(
"VCS_SOLVE::vcs_basopt",
"we have a prob");
2457 plogf(
" --- largest error in Stoich coeff = %g at rxn = %d ", sumMax, iMax);
2458 plogf(
"%-10.10s", m_speciesName[m_indexRxnToSpecies[iMax]]);
2459 plogf(
" element = %d ", jMax);
2460 plogf(
"%-5.5s", m_elementName[jMax]);
2463 for (
size_t i=0; i<77; i++) {
2475 m_deltaMolNumPhase.zero();
2476 m_phaseParticipation.zero();
2481 for (
size_t irxn = 0; irxn < m_numRxnTot; ++irxn) {
2482 double* scrxn_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
2483 size_t kspec = m_indexRxnToSpecies[irxn];
2484 size_t iph = m_phaseID[kspec];
2485 m_deltaMolNumPhase(iph,irxn) = 1.0;
2486 m_phaseParticipation(iph,irxn)++;
2487 for (
size_t j = 0; j < ncTrial; ++j) {
2489 if (fabs(scrxn_ptr[j]) <= 1.0e-6) {
2492 m_deltaMolNumPhase(iph,irxn) += scrxn_ptr[j];
2493 m_phaseParticipation(iph,irxn)++;
2501 m_VCount->Time_basopt += tsecond;
2502 m_VCount->Basis_Opts++;
2506 size_t VCS_SOLVE::vcs_basisOptMax(
const double*
const molNum,
const size_t j,
2523 double big = molNum[j] * m_spSize[j] * 1.01;
2524 if (m_spSize[j] <= 0.0) {
2526 "spSize is nonpositive");
2528 for (
size_t i = j + 1; i < n; ++i) {
2529 if (m_spSize[i] <= 0.0) {
2531 "spSize is nonpositive");
2533 bool doSwap =
false;
2535 doSwap = (molNum[i] * m_spSize[i]) > big;
2536 if (!m_SSPhase[i] && doSwap) {
2537 doSwap = molNum[i] > (molNum[largest] * 1.001);
2541 doSwap = (molNum[i] * m_spSize[i]) > big;
2543 doSwap = molNum[i] > (molNum[largest] * 1.001);
2546 doSwap = (molNum[i] * m_spSize[i]) > big;
2551 big = molNum[i] * m_spSize[i] * 1.01;
2557 int VCS_SOLVE::vcs_species_type(
const size_t kspec)
const 2564 size_t iph = m_phaseID[kspec];
2565 int irxn = int(kspec) - int(m_numComponents);
2566 int phaseExist = m_VolPhaseList[iph]->exists();
2569 if (m_molNumSpecies_old[kspec] <= 0.0) {
2570 if (m_tPhaseMoles_old[iph] <= 0.0 && !m_SSPhase[kspec]) {
2576 for (
size_t j = 0; j < m_nelem; ++j) {
2578 double atomComp = m_formulaMatrix(kspec,j);
2579 if (atomComp > 0.0) {
2580 double maxPermissible = m_elemAbundancesGoal[j] / atomComp;
2582 if (m_debug_print_lvl >= 2) {
2583 plogf(
" --- %s can not be nonzero because" 2584 " needed element %s is zero\n",
2585 m_speciesName[kspec], m_elementName[j]);
2587 if (m_SSPhase[kspec]) {
2602 for (
size_t j = 0; j < m_numComponents; ++j) {
2603 double stoicC = m_stoichCoeffRxnMatrix(j,irxn);
2604 if (stoicC != 0.0) {
2605 double negChangeComp = - stoicC;
2606 if (negChangeComp > 0.0) {
2607 if (m_molNumSpecies_old[j] < 1.0E-60) {
2608 if (m_debug_print_lvl >= 2) {
2609 plogf(
" --- %s is prevented from popping into existence because" 2610 " a needed component to be consumed, %s, has a zero mole number\n",
2611 m_speciesName[kspec], m_speciesName[j]);
2613 if (m_SSPhase[kspec]) {
2619 }
else if (negChangeComp < 0.0) {
2620 if (m_VolPhaseList[m_phaseID[j]]->exists() <= 0) {
2621 if (m_debug_print_lvl >= 2) {
2622 plogf(
" --- %s is prevented from popping into existence because" 2623 " a needed component %s is in a zeroed-phase that would be " 2624 "popped into existence at the same time\n",
2625 m_speciesName[kspec], m_speciesName[j]);
2627 if (m_SSPhase[kspec]) {
2638 if (irxn >= 0 && m_deltaGRxn_old[irxn] >= 0.0) {
2640 if (m_SSPhase[kspec]) {
2654 if (m_tPhaseMoles_old[iph] > 0.0) {
2655 if (m_SSPhase[kspec]) {
2666 if (m_tPhaseMoles_old[iph] <= 0.0) {
2667 if (m_SSPhase[kspec]) {
2679 if (m_SSPhase[kspec]) {
2688 if (m_molNumSpecies_old[kspec] > (m_tPhaseMoles_old[iph] * 0.001)) {
2700 double szAdj = m_scSize[irxn] * std::sqrt((
double)m_numRxnTot);
2701 for (
size_t k = 0; k < m_numComponents; ++k) {
2702 if (!m_SSPhase[k] && m_stoichCoeffRxnMatrix(k,irxn) != 0.0 && m_molNumSpecies_old[kspec] * szAdj >= m_molNumSpecies_old[k] * 0.01) {
2710 void VCS_SOLVE::vcs_dfe(
const int stateCalc,
2711 const int ll,
const size_t lbot,
const size_t ltop)
2713 double* tPhMoles_ptr=0;
2714 double* actCoeff_ptr=0;
2715 double* feSpecies=0;
2718 feSpecies = &m_feSpecies_old[0];
2719 tPhMoles_ptr = &m_tPhaseMoles_old[0];
2720 actCoeff_ptr = &m_actCoeffSpecies_old[0];
2721 molNum = &m_molNumSpecies_old[0];
2723 feSpecies = &m_feSpecies_new[0];
2724 tPhMoles_ptr = &m_tPhaseMoles_new[0];
2725 actCoeff_ptr = &m_actCoeffSpecies_new[0];
2726 molNum = &m_molNumSpecies_new[0];
2729 "Subroutine vcs_dfe called with bad stateCalc value: {}", stateCalc);
2732 if (m_debug_print_lvl >= 2) {
2735 plogf(
" --- Subroutine vcs_dfe called for one species: ");
2736 plogf(
"%-12.12s", m_speciesName[lbot]);
2738 plogf(
" --- Subroutine vcs_dfe called for all species");
2740 }
else if (ll > 0) {
2741 plogf(
" --- Subroutine vcs_dfe called for components and minors");
2743 plogf(
" --- Subroutine vcs_dfe called for components and majors");
2746 plogf(
" using tentative solution");
2751 double* tlogMoles = &m_TmpPhase[0];
2755 double* tPhInertMoles = &TPhInertMoles[0];
2756 for (
size_t iph = 0; iph < m_numPhases; iph++) {
2757 tlogMoles[iph] = tPhInertMoles[iph];
2760 for (
size_t kspec = 0; kspec < m_nsp; kspec++) {
2762 size_t iph = m_phaseID[kspec];
2763 tlogMoles[iph] += molNum[kspec];
2766 for (
size_t iph = 0; iph < m_numPhases; iph++) {
2768 "VCS_SOLVE::vcs_dfe",
2769 "phase Moles may be off, iph = {}, {} {}",
2770 iph, tlogMoles[iph], tPhMoles_ptr[iph]);
2772 m_TmpPhase.assign(m_TmpPhase.size(), 0.0);
2773 for (
size_t iph = 0; iph < m_numPhases; iph++) {
2774 if (tPhMoles_ptr[iph] > 0.0) {
2775 tlogMoles[iph] = log(tPhMoles_ptr[iph]);
2782 l2 = m_numComponents;
2791 for (
size_t iphase = 0; iphase < m_numPhases; iphase++) {
2806 for (
size_t kspec = l1; kspec < l2; ++kspec) {
2807 size_t iphase = m_phaseID[kspec];
2809 AssertThrowMsg(molNum[kspec] == m_phasePhi[iphase],
"VCS_SOLVE::vcs_dfe",
2810 "We have an inconsistency!");
2811 AssertThrowMsg(m_chargeSpecies[kspec] == -1.0,
"VCS_SOLVE::vcs_dfe",
2812 "We have an unexpected situation!");
2813 feSpecies[kspec] = m_SSfeSpecies[kspec]
2814 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2816 if (m_SSPhase[kspec]) {
2817 feSpecies[kspec] = m_SSfeSpecies[kspec]
2818 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2821 feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
2822 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2825 size_t iph = m_phaseID[kspec];
2826 if (tPhMoles_ptr[iph] > 0.0) {
2827 feSpecies[kspec] = m_SSfeSpecies[kspec]
2829 - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
2830 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2832 feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
2833 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2836 feSpecies[kspec] = m_SSfeSpecies[kspec]
2837 + log(actCoeff_ptr[kspec] * molNum[kspec])
2838 - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
2839 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2847 for (
size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
2848 size_t kspec = m_indexRxnToSpecies[irxn];
2850 size_t iphase = m_phaseID[kspec];
2852 AssertThrowMsg(molNum[kspec] == m_phasePhi[iphase],
"VCS_SOLVE::vcs_dfe",
2853 "We have an inconsistency!");
2854 AssertThrowMsg(m_chargeSpecies[kspec] == -1.0,
"VCS_SOLVE::vcs_dfe",
2855 "We have an unexpected situation!");
2856 feSpecies[kspec] = m_SSfeSpecies[kspec]
2857 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2859 if (m_SSPhase[kspec]) {
2860 feSpecies[kspec] = m_SSfeSpecies[kspec]
2861 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2864 feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
2865 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2868 size_t iph = m_phaseID[kspec];
2869 if (tPhMoles_ptr[iph] > 0.0) {
2870 feSpecies[kspec] = m_SSfeSpecies[kspec]
2872 - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
2873 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2875 feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
2876 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2879 feSpecies[kspec] = m_SSfeSpecies[kspec]
2880 + log(actCoeff_ptr[kspec] * molNum[kspec])
2881 - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
2882 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2888 }
else if (ll > 0) {
2890 for (
size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
2891 size_t kspec = m_indexRxnToSpecies[irxn];
2893 size_t iphase = m_phaseID[kspec];
2895 AssertThrowMsg(molNum[kspec] == m_phasePhi[iphase],
"VCS_SOLVE::vcs_dfe",
2896 "We have an inconsistency!");
2897 AssertThrowMsg(m_chargeSpecies[kspec] == -1.0,
"VCS_SOLVE::vcs_dfe",
2898 "We have an unexpected situation!");
2899 feSpecies[kspec] = m_SSfeSpecies[kspec]
2900 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2902 if (m_SSPhase[kspec]) {
2903 feSpecies[kspec] = m_SSfeSpecies[kspec]
2904 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2907 feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
2908 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2911 size_t iph = m_phaseID[kspec];
2912 if (tPhMoles_ptr[iph] > 0.0) {
2913 feSpecies[kspec] = m_SSfeSpecies[kspec]
2915 - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
2916 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2918 feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
2919 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2922 feSpecies[kspec] = m_SSfeSpecies[kspec]
2923 + log(actCoeff_ptr[kspec] * molNum[kspec])
2924 - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
2925 + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
2934 double VCS_SOLVE::l2normdg(
double dgLocal[])
const 2936 if (m_numRxnRdc <= 0) {
2940 for (
size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
2941 size_t kspec = irxn + m_numComponents;
2943 dgLocal[irxn] < 0.0) {
2945 tmp += dgLocal[irxn] * dgLocal[irxn];
2949 return std::sqrt(tmp / m_numRxnRdc);
2952 double VCS_SOLVE::vcs_tmoles()
2954 for (
size_t i = 0; i < m_numPhases; i++) {
2955 m_tPhaseMoles_old[i] = TPhInertMoles[i];
2957 for (
size_t i = 0; i < m_nsp; i++) {
2959 m_tPhaseMoles_old[m_phaseID[i]] += m_molNumSpecies_old[i];
2963 for (
size_t i = 0; i < m_numPhases; i++) {
2964 sum += m_tPhaseMoles_old[i];
2966 if (m_tPhaseMoles_old[i] == 0.0) {
2972 m_totalMolNum = sum;
2973 return m_totalMolNum;
2976 void VCS_SOLVE::check_tmoles()
const 2979 for (
size_t i = 0; i < m_numPhases; i++) {
2980 double m_tPhaseMoles_old_a = TPhInertMoles[i];
2982 for (
size_t k = 0; k < m_nsp; k++) {
2984 m_tPhaseMoles_old_a += m_molNumSpecies_old[k];
2987 sum += m_tPhaseMoles_old_a;
2989 double denom = m_tPhaseMoles_old[i]+ m_tPhaseMoles_old_a + 1.0E-19;
2990 if (!
vcs_doubleEqual(m_tPhaseMoles_old[i]/denom, m_tPhaseMoles_old_a/denom)) {
2991 plogf(
"check_tmoles: we have found a problem with phase %d: %20.15g, %20.15g\n",
2992 i, m_tPhaseMoles_old[i], m_tPhaseMoles_old_a);
2997 void VCS_SOLVE::vcs_updateVP(
const int vcsState)
2999 for (
size_t i = 0; i < m_numPhases; i++) {
3003 &m_molNumSpecies_old[0],
3004 &m_tPhaseMoles_old[0]);
3007 &m_molNumSpecies_new[0],
3008 &m_tPhaseMoles_new[0]);
3011 "wrong stateCalc value: {}", vcsState);
3016 bool VCS_SOLVE::vcs_evaluate_speciesType()
3018 m_numRxnMinorZeroed = 0;
3019 if (m_debug_print_lvl >= 2) {
3020 plogf(
" --- Species Status decision is reevaluated: All species are minor except for:\n");
3021 }
else if (m_debug_print_lvl >= 5) {
3022 plogf(
" --- Species Status decision is reevaluated\n");
3024 for (
size_t kspec = 0; kspec < m_nsp; ++kspec) {
3025 m_speciesStatus[kspec] = vcs_species_type(kspec);
3026 if (m_debug_print_lvl >= 5) {
3027 plogf(
" --- %-16s: ", m_speciesName[kspec]);
3028 if (kspec < m_numComponents) {
3033 plogf(
" %10.3g ", m_molNumSpecies_old[kspec]);
3035 plogf(
"%s\n", sString);
3036 }
else if (m_debug_print_lvl >= 2) {
3038 switch (m_speciesStatus[kspec]) {
3042 plogf(
" --- Major Species : %-s\n", m_speciesName[kspec]);
3045 plogf(
" --- Purposely Zeroed-Phase Species (not in problem): %-s\n",
3046 m_speciesName[kspec]);
3049 plogf(
" --- Zeroed-MS Phase Species: %-s\n", m_speciesName[kspec]);
3052 plogf(
" --- Zeroed-SS Phase Species: %-s\n", m_speciesName[kspec]);
3055 plogf(
" --- Deleted-Small Species : %-s\n", m_speciesName[kspec]);
3058 plogf(
" --- Zeroed Species in an active MS phase (tmp): %-s\n",
3059 m_speciesName[kspec]);
3062 plogf(
" --- Zeroed Species in an active MS phase (Stoich Constraint): %-s\n",
3063 m_speciesName[kspec]);
3066 plogf(
" --- InterfaceVoltage Species: %-s\n", m_speciesName[kspec]);
3069 throw CanteraError(
"VCS_SOLVE::vcs_evaluate_speciesType",
3070 "Unknown type: {}", m_speciesStatus[kspec]);
3075 ++m_numRxnMinorZeroed;
3078 debuglog(
" ---\n", m_debug_print_lvl >= 2);
3079 return (m_numRxnMinorZeroed >= m_numRxnRdc);
3082 void VCS_SOLVE::vcs_deltag(
const int L,
const bool doDeleted,
3083 const int vcsState,
const bool alterZeroedPhases)
3085 size_t irxnl = m_numRxnRdc;
3087 irxnl = m_numRxnTot;
3092 double* molNumSpecies;
3093 double* actCoeffSpecies;
3095 deltaGRxn = &m_deltaGRxn_new[0];
3096 feSpecies = &m_feSpecies_new[0];
3097 molNumSpecies = &m_molNumSpecies_new[0];
3098 actCoeffSpecies = &m_actCoeffSpecies_new[0];
3100 deltaGRxn = &m_deltaGRxn_old[0];
3101 feSpecies = &m_feSpecies_old[0];
3102 molNumSpecies = &m_molNumSpecies_old[0];
3103 actCoeffSpecies = &m_actCoeffSpecies_old[0];
3105 throw CanteraError(
"VCS_SOLVE::vcs_deltag",
"bad vcsState");
3108 if (m_debug_print_lvl >= 2) {
3109 plogf(
" --- Subroutine vcs_deltag called for ");
3111 plogf(
"major noncomponents\n");
3112 }
else if (L == 0) {
3113 plogf(
"all noncomponents\n");
3115 plogf(
"minor noncomponents\n");
3121 for (
size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
3122 size_t kspec = irxn + m_numComponents;
3125 deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]];
3126 double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
3127 for (kspec = 0; kspec < m_numComponents; ++kspec) {
3128 deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec];
3134 deltaGRxn[irxn] = std::max(0.0, deltaGRxn[irxn]);
3138 }
else if (L == 0) {
3140 for (
size_t irxn = 0; irxn < irxnl; ++irxn) {
3142 deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]];
3143 double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
3144 for (
size_t kspec = 0; kspec < m_numComponents; ++kspec) {
3145 deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec];
3147 dtmp_ptr[kspec] < 0.0) {
3152 deltaGRxn[irxn] = std::max(0.0, deltaGRxn[irxn]);
3157 for (
size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
3158 size_t kspec = irxn + m_numComponents;
3161 deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]];
3162 double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
3163 for (kspec = 0; kspec < m_numComponents; ++kspec) {
3164 deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec];
3166 dtmp_ptr[kspec] < 0.0) {
3171 deltaGRxn[irxn] = std::max(0.0, deltaGRxn[irxn]);
3216 if (alterZeroedPhases &&
false) {
3217 for (
size_t iph = 0; iph < m_numPhases; iph++) {
3222 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
3225 sum += molNumSpecies[kspec];
3238 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
3241 if (kspec >= m_numComponents) {
3242 size_t irxn = kspec - m_numComponents;
3243 deltaGRxn[irxn] =
clip(deltaGRxn[irxn], -50.0, 50.0);
3244 poly += exp(-deltaGRxn[irxn])/actCoeffSpecies[kspec];
3252 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
3254 if (kspec >= m_numComponents) {
3255 size_t irxn = kspec - m_numComponents;
3256 deltaGRxn[irxn] = 1.0 - poly;
3264 void VCS_SOLVE::vcs_printDeltaG(
const int stateCalc)
3266 double* deltaGRxn = &m_deltaGRxn_old[0];
3267 double* feSpecies = &m_feSpecies_old[0];
3268 double* molNumSpecies = &m_molNumSpecies_old[0];
3269 const double* tPhMoles_ptr = &m_tPhaseMoles_old[0];
3270 const double* actCoeff_ptr = &m_actCoeffSpecies_old[0];
3272 deltaGRxn = &m_deltaGRxn_new[0];
3273 feSpecies = &m_feSpecies_new[0];
3274 molNumSpecies = &m_molNumSpecies_new[0];
3275 actCoeff_ptr = &m_actCoeffSpecies_new[0];
3276 tPhMoles_ptr = &m_tPhaseMoles_new[0];
3279 bool zeroedPhase =
false;
3280 if (m_debug_print_lvl >= 2) {
3281 plogf(
" --- DELTA_G TABLE Components:");
3282 for (
size_t j = 0; j < m_numComponents; j++) {
3285 plogf(
"\n --- Components Moles:");
3286 for (
size_t j = 0; j < m_numComponents; j++) {
3287 plogf(
"%10.3g", m_molNumSpecies_old[j]);
3289 plogf(
"\n --- NonComponent| Moles | ");
3290 for (
size_t j = 0; j < m_numComponents; j++) {
3291 plogf(
"%-10.10s", m_speciesName[j]);
3294 for (
size_t i = 0; i < m_numRxnTot; i++) {
3295 plogf(
" --- %3d ", m_indexRxnToSpecies[i]);
3296 plogf(
"%-10.10s", m_speciesName[m_indexRxnToSpecies[i]]);
3300 plogf(
"|%10.3g|", m_molNumSpecies_old[m_indexRxnToSpecies[i]]);
3302 for (
size_t j = 0; j < m_numComponents; j++) {
3303 plogf(
" %6.2f", m_stoichCoeffRxnMatrix(j,i));
3308 for (
int i=0; i<77; i++) {
3314 writelog(
" --- DeltaG Table (J/kmol) Name PhID MoleNum MolFR " 3315 " ElectrChemStar ElectrChem DeltaGStar DeltaG(Pred) Stability\n");
3317 writeline(
'-', 132);
3319 for (
size_t kspec = 0; kspec < m_nsp; kspec++) {
3321 if (kspec >= m_numComponents) {
3322 irxn = kspec - m_numComponents;
3324 double mfValue = 1.0;
3325 size_t iphase = m_phaseID[kspec];
3326 const vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get();
3332 zeroedPhase =
false;
3334 if (tPhMoles_ptr[iphase] > 0.0) {
3338 mfValue = molNumSpecies[kspec] / tPhMoles_ptr[iphase];
3341 size_t klocal = m_speciesLocalPhaseIndex[kspec];
3342 mfValue = Vphase->moleFraction(klocal);
3349 double feFull = feSpecies[kspec];
3352 feFull += log(actCoeff_ptr[kspec]) + log(mfValue);
3354 writelogf(
"%-24.24s", m_speciesName[kspec]);
3359 writelogf(
" % -12.4e", molNumSpecies[kspec]);
3362 writelogf(
" % -12.4e", feSpecies[kspec] * RT);
3365 writelogf(
" % -12.4e", deltaGRxn[irxn] * RT);
3366 writelogf(
" % -12.4e", (deltaGRxn[irxn] + feFull - feSpecies[kspec]) * RT);
3368 if (deltaGRxn[irxn] < 0.0) {
3369 if (molNumSpecies[kspec] > 0.0) {
3374 }
else if (deltaGRxn[irxn] > 0.0) {
3375 if (molNumSpecies[kspec] > 0.0) {
3387 writeline(
'-', 132);
3390 void VCS_SOLVE::vcs_switch_pos(
const bool ifunc,
const size_t k1,
const size_t k2)
3395 if (k1 >= m_nsp || k2 >= m_nsp) {
3396 plogf(
"vcs_switch_pos: ifunc = 0: inappropriate args: %d %d\n",
3401 vcs_VolPhase* pv1 = m_VolPhaseList[m_phaseID[k1]].get();
3402 vcs_VolPhase* pv2 = m_VolPhaseList[m_phaseID[k2]].get();
3403 size_t kp1 = m_speciesLocalPhaseIndex[k1];
3404 size_t kp2 = m_speciesLocalPhaseIndex[k2];
3411 std::swap(m_speciesName[k1], m_speciesName[k2]);
3412 std::swap(m_molNumSpecies_old[k1], m_molNumSpecies_old[k2]);
3413 std::swap(m_speciesUnknownType[k1], m_speciesUnknownType[k2]);
3414 std::swap(m_molNumSpecies_new[k1], m_molNumSpecies_new[k2]);
3415 std::swap(m_SSfeSpecies[k1], m_SSfeSpecies[k2]);
3416 std::swap(m_spSize[k1], m_spSize[k2]);
3417 std::swap(m_deltaMolNumSpecies[k1], m_deltaMolNumSpecies[k2]);
3418 std::swap(m_feSpecies_old[k1], m_feSpecies_old[k2]);
3419 std::swap(m_feSpecies_new[k1], m_feSpecies_new[k2]);
3420 std::swap(m_SSPhase[k1], m_SSPhase[k2]);
3421 std::swap(m_phaseID[k1], m_phaseID[k2]);
3422 std::swap(m_speciesMapIndex[k1], m_speciesMapIndex[k2]);
3423 std::swap(m_speciesLocalPhaseIndex[k1], m_speciesLocalPhaseIndex[k2]);
3424 std::swap(m_actConventionSpecies[k1], m_actConventionSpecies[k2]);
3425 std::swap(m_lnMnaughtSpecies[k1], m_lnMnaughtSpecies[k2]);
3426 std::swap(m_actCoeffSpecies_new[k1], m_actCoeffSpecies_new[k2]);
3427 std::swap(m_actCoeffSpecies_old[k1], m_actCoeffSpecies_old[k2]);
3428 std::swap(m_wtSpecies[k1], m_wtSpecies[k2]);
3429 std::swap(m_chargeSpecies[k1], m_chargeSpecies[k2]);
3430 std::swap(m_speciesThermoList[k1], m_speciesThermoList[k2]);
3431 std::swap(m_PMVolumeSpecies[k1], m_PMVolumeSpecies[k2]);
3433 for (
size_t j = 0; j < m_nelem; ++j) {
3434 std::swap(m_formulaMatrix(k1,j), m_formulaMatrix(k2,j));
3436 if (m_useActCoeffJac && k1 != k2) {
3437 for (
size_t i = 0; i < m_nsp; i++) {
3438 std::swap(m_np_dLnActCoeffdMolNum(k1,i), m_np_dLnActCoeffdMolNum(k2,i));
3440 for (
size_t i = 0; i < m_nsp; i++) {
3441 std::swap(m_np_dLnActCoeffdMolNum(i,k1), m_np_dLnActCoeffdMolNum(i,k2));
3444 std::swap(m_speciesStatus[k1], m_speciesStatus[k2]);
3449 size_t i1 = k1 - m_numComponents;
3450 size_t i2 = k2 - m_numComponents;
3451 if (i1 > m_numRxnTot || i2 >= m_numRxnTot) {
3452 plogf(
"switch_pos: ifunc = 1: inappropriate noncomp values: %d %d\n",
3455 for (
size_t j = 0; j < m_numComponents; ++j) {
3456 std::swap(m_stoichCoeffRxnMatrix(j,i1), m_stoichCoeffRxnMatrix(j,i2));
3458 std::swap(m_scSize[i1], m_scSize[i2]);
3459 for (
size_t iph = 0; iph < m_numPhases; iph++) {
3460 std::swap(m_deltaMolNumPhase(iph,i1), m_deltaMolNumPhase(iph,i2));
3461 std::swap(m_phaseParticipation(iph,i1),
3462 m_phaseParticipation(iph,i2));
3464 std::swap(m_deltaGRxn_new[i1], m_deltaGRxn_new[i2]);
3465 std::swap(m_deltaGRxn_old[i1], m_deltaGRxn_old[i2]);
3466 std::swap(m_deltaGRxn_tmp[i1], m_deltaGRxn_tmp[i2]);
3470 void VCS_SOLVE::vcs_setFlagsVolPhases(
const bool upToDate,
const int stateCalc)
3473 for (
size_t iph = 0; iph < m_numPhases; iph++) {
3474 m_VolPhaseList[iph]->setMolesOutOfDate(stateCalc);
3477 for (
size_t iph = 0; iph < m_numPhases; iph++) {
3478 m_VolPhaseList[iph]->setMolesCurrent(stateCalc);
3483 void VCS_SOLVE::vcs_setFlagsVolPhase(
const size_t iph,
const bool upToDate,
3484 const int stateCalc)
3487 m_VolPhaseList[iph]->setMolesOutOfDate(stateCalc);
3489 m_VolPhaseList[iph]->setMolesCurrent(stateCalc);
3493 void VCS_SOLVE::vcs_updateMolNumVolPhases(
const int stateCalc)
3495 for (
size_t iph = 0; iph < m_numPhases; iph++) {
3496 m_VolPhaseList[iph]->updateFromVCS_MoleNumbers(stateCalc);
#define VCS_PHASE_EXIST_NO
Phase doesn't currently exist in the mixture.
bool vcs_doubleEqual(double d1, double d2)
Simple routine to check whether two doubles are equal up to roundoff error.
bool m_singleSpecies
If true, this phase consists of a single species.
void updateFromVCS_MoleNumbers(const int stateCalc)
Update the moles within the phase, if necessary.
#define VCS_SPECIES_ZEROEDMS
Species lies in a multicomponent phase with concentration zero.
#define VCS_SPECIES_MINOR
Species is a major species.
Declarations for a simple class that implements an Ansi C wall clock timer (see Cantera::clockWC).
#define VCS_PHASE_EXIST_ALWAYS
Always exists because it contains inerts which can't exist in any other phase.
void setMolesFromVCSCheck(const int vcsStateStatus, const double *molesSpeciesVCS, const double *const TPhMoles)
Set the moles within the phase.
#define VCS_PHASE_EXIST_ZEROEDPHASE
Phase currently is zeroed due to a programmatic issue.
const size_t npos
index returned by functions to indicate "no position"
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
double electricPotential() const
Returns the electric field of the phase.
T clip(const T &value, const T &lower, const T &upper)
Clip value such that lower <= value <= upper.
size_t spGlobalIndexVCS(const size_t spIndex) const
Return the Global VCS index of the kth species in the phase.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
The class provides the wall clock timer in seconds.
#define VCS_SPECIES_INTERFACIALVOLTAGE
Species refers to an electron in the metal.
std::string PhaseName
String name for the phase.
doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
#define VCS_PHASE_EXIST_YES
Phase is a normal phase that currently exists.
void setExistence(const int existence)
Set the existence flag in the object.
std::vector< int > vector_int
Vector of ints.
Header file for the internal object that holds the vcs equilibrium problem (see Class VCS_SOLVE and E...
Header for the object representing each phase within vcs.
#define VCS_SPECIES_ZEROEDSS
Species is a SS phase, that is currently zeroed out.
#define VCS_SPECIES_TYPE_MOLNUM
Unknown refers to mole number of a single species.
#define VCS_ELEM_TYPE_ABSPOS
Normal element constraint consisting of positive coefficients for the formula matrix.
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
#define VCS_SPECIES_ZEROEDPHASE
Species lies in a multicomponent phase that is zeroed atm.
#define VCS_DELETE_MINORSPECIES_CUTOFF
Cutoff relative mole number value, below which species are deleted from the equilibrium problem...
double secondsWC()
Returns the wall clock time in seconds since the last reset.
#define VCS_DELETE_PHASE_CUTOFF
Cutoff relative moles below which a phase is deleted from the equilibrium problem.
#define VCS_STATECALC_NEW
State Calculation based on the new or tentative mole numbers.
size_t nSpecies() const
Return the number of species in the phase.
Base class for exceptions thrown by Cantera classes.
#define VCS_SPECIES_ACTIVEBUTZERO
Species lies in a multicomponent phase that is active, but species concentration is zero...
void sendToVCS_ActCoeff(const int stateCalc, double *const AC)
Fill in an activity coefficients vector within a VCS_SOLVE object.
#define VCS_SPECIES_MAJOR
Species is a major species.
#define VCS_RELDELETE_SPECIES_CUTOFF
Cutoff relative mole fraction value, below which species are deleted from the equilibrium problem...
#define VCS_SPECIES_STOICHZERO
Species lies in a multicomponent phase that is active, but species concentration is zero due to stoic...
void debuglog(const std::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.
Phase information and Phase calculations for vcs.
void writelogendl()
Write an end of line character to the screen and flush output.
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
#define VCS_STATECALC_OLD
State Calculation based on the old or base mole numbers.
int exists() const
int indicating whether the phase exists or not
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Headers for the DenseMatrix object, which deals with dense rectangular matrices and description of th...
void setSpGlobalIndexVCS(const size_t spIndex, const size_t spGlobalIndex)
set the Global VCS index of the kth species in the phase
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Contains declarations for string manipulation functions within Cantera.
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...
#define plogf
define this Cantera function to replace printf
#define VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
Unknown refers to the voltage level of a phase.
#define VCS_SPECIES_COMPONENT
Species is a component which can be nonzero.
Namespace for the Cantera kernel.
#define VCS_SPECIES_DELETED
Species has such a small mole fraction it is deleted even though its phase may possibly exist...
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
void setTotalMoles(const double totalMols)
Sets the total moles in the phase.