8 #include "vcs_species_thermo.h"
40 bool VCS_SOLVE::vcs_popPhasePossible(
const size_t iphasePop)
const
46 int existence = Vphase->
exists();
48 printf(
"ERROR vcs_popPhasePossible called for a phase that exists!");
58 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
61 if (m_molNumSpecies_old[kspec] > 0.0) {
62 printf(
"ERROR vcs_popPhasePossible we shouldn't be here %lu %g > 0.0",
63 kspec, m_molNumSpecies_old[kspec]);
67 size_t irxn = kspec - m_numComponents;
68 if (kspec >= m_numComponents) {
69 bool iPopPossible =
true;
70 for (
size_t j = 0; j < m_numComponents; ++j) {
72 double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
74 double negChangeComp = - stoicC * 1.0;
75 if (negChangeComp > 0.0) {
94 for (
size_t jrxn = 0; jrxn < m_numRxnRdc; jrxn++) {
95 bool foundJrxn =
false;
97 if (m_stoichCoeffRxnMatrix[jrxn][kspec] > 0.0) {
99 for (
size_t kcomp = 0; kcomp < m_numComponents; kcomp++) {
100 if (m_stoichCoeffRxnMatrix[jrxn][kcomp] < 0.0) {
112 else if (m_stoichCoeffRxnMatrix[jrxn][kspec] < 0.0) {
114 size_t jspec = jrxn + m_numComponents;
119 for (
size_t kcomp = 0; kcomp < m_numComponents; kcomp++) {
120 if (m_stoichCoeffRxnMatrix[jrxn][kcomp] > 0.0) {
147 int VCS_SOLVE::vcs_phasePopDeterminePossibleList()
153 std::vector<int> linkedPhases;
154 phasePopProblemLists_.clear();
166 std::vector< std::vector<size_t> > zeroedComponentLinkedPhasePops(m_numComponents);
170 for (
size_t j = 0; j < m_numComponents; j++) {
172 molComp = m_molNumSpecies_old[j];
173 if (molComp <= 0.0) {
174 std::vector<size_t> &jList = zeroedComponentLinkedPhasePops[j];
175 size_t iph = m_phaseID[j];
176 jList.push_back(iph);
177 for (
size_t irxn = 0; irxn < m_numRxnTot; irxn++) {
178 size_t kspec = irxn + m_numComponents;
179 iph = m_phaseID[kspec];
180 Vphase = m_VolPhaseList[iph];
181 int existence = Vphase->
exists();
183 stoicC = m_stoichCoeffRxnMatrix[irxn][j];
185 if (std::find(jList.begin(), jList.end(), iph) != jList.end()) {
186 jList.push_back(iph);
200 std::vector< std::vector<size_t> > zeroedPhaseLinkedZeroComponents(m_numPhases);
204 for (
size_t iph = 0; iph < m_numPhases; iph++) {
205 std::vector<size_t> &iphList = zeroedPhaseLinkedZeroComponents[iph];
207 Vphase = m_VolPhaseList[iph];
208 int existence = Vphase->
exists();
211 linkedPhases.clear();
213 for (
size_t k = 0; k < nsp; k++) {
215 size_t irxn = kspec - m_numComponents;
217 for (
size_t j = 0; j < m_numComponents; j++) {
219 molComp = m_molNumSpecies_old[j];
220 if (molComp <= 0.0) {
221 stoicC = m_stoichCoeffRxnMatrix[irxn][j];
223 bool foundPos =
false;
224 for (
size_t kk = 0; kk < nsp; kk++) {
226 if (kkspec >= m_numComponents) {
227 size_t iirxn = kkspec - m_numComponents;
228 if (m_stoichCoeffRxnMatrix[iirxn][j] > 0.0) {
234 if (std::find(iphList.begin(), iphList.end(), j) != iphList.end()) {
235 iphList.push_back(j);
250 for (
size_t iph = 0; iph < m_numPhases; iph++) {
251 Vphase = m_VolPhaseList[iph];
252 int existence = Vphase->
exists();
254 std::vector<size_t> &iphList = zeroedPhaseLinkedZeroComponents[iph];
255 std::vector<size_t> popProblem(0);
256 popProblem.push_back(iph);
257 for (
size_t i = 0; i < iphList.size(); i++) {
258 size_t j = iphList[i];
259 std::vector<size_t> &jList = zeroedComponentLinkedPhasePops[j];
260 for (
size_t jjl = 0; jjl < jList.size(); jjl++) {
261 size_t jph = jList[jjl];
262 if (std::find(popProblem.begin(), popProblem.end(), jph) != popProblem.end()) {
263 popProblem.push_back(jph);
267 phasePopProblemLists_.push_back(popProblem);
281 size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
283 size_t iphasePop =
npos;
285 doublereal FephaseMax = -1.0E30;
286 doublereal Fephase = -1.0E30;
292 if (m_debug_print_lvl >= 2) {
293 plogf(
" --- vcs_popPhaseID() called\n");
294 plogf(
" --- Phase Status F_e MoleNum\n");
295 plogf(
" --------------------------------------------------------------------------\n");
298 for (
size_t iph = 0; iph < m_numPhases; iph++) {
299 Vphase = m_VolPhaseList[iph];
300 int existence = Vphase->
exists();
307 if (m_debug_print_lvl >= 2) {
308 plogf(
" --- %18s %5d NA %11.3e\n",
311 m_tPhaseMoles_old[iph]);
322 irxn = kspec - m_numComponents;
323 doublereal deltaGRxn = m_deltaGRxn_old[irxn];
324 Fephase = exp(-deltaGRxn) - 1.0;
327 strcpy(anote,
" (ready to be birthed)");
329 if (Fephase > FephaseMax) {
331 FephaseMax = Fephase;
333 strcpy(anote,
" (chosen to be birthed)");
339 strcpy(anote,
" (not stable)");
340 if (m_tPhaseMoles_old[iph] > 0.0) {
341 printf(
"shouldn't be here\n");
348 if (m_debug_print_lvl >= 2) {
349 plogf(
" --- %18s %5d %10.3g %10.3g %s\n",
352 m_tPhaseMoles_old[iph], anote);
362 if (vcs_popPhasePossible(iph)) {
363 Fephase = vcs_phaseStabilityTest(iph);
365 if (Fephase > FephaseMax) {
367 FephaseMax = Fephase;
370 if (Fephase > FephaseMax) {
371 FephaseMax = Fephase;
375 if (m_debug_print_lvl >= 2) {
376 plogf(
" --- %18s %5d %11.3g %11.3g\n",
379 m_tPhaseMoles_old[iph]);
384 if (m_debug_print_lvl >= 2) {
385 plogf(
" --- %18s %5d blocked %11.3g\n",
387 existence, m_tPhaseMoles_old[iph]);
394 phasePopPhaseIDs.resize(0);
395 if (iphasePop !=
npos) {
396 phasePopPhaseIDs.push_back(iphasePop);
405 if (m_debug_print_lvl >= 2) {
406 plogf(
" ---------------------------------------------------------------------\n");
434 int VCS_SOLVE::vcs_popPhaseRxnStepSizes(
const size_t iphasePop)
440 size_t irxn = kspec - m_numComponents;
441 std::vector<size_t> creationGlobalRxnNumbers;
451 int existence = Vphase->
exists();
453 printf(
"ERROR vcs_popPhaseRxnStepSizes called for a phase that exists!");
457 if (m_debug_print_lvl >= 2) {
458 plogf(
" --- vcs_popPhaseRxnStepSizes() called to pop phase %s %d into existence\n",
465 double* dnPhase_irxn = m_deltaMolNumPhase[irxn];
466 for (
size_t j = 0; j < m_numComponents; ++j) {
468 if (m_molNumSpecies_old[j] > 0.0) {
469 s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j];
473 for (
size_t j = 0; j < m_numPhases; j++) {
474 Vphase = m_VolPhaseList[j];
476 if (m_tPhaseMoles_old[j] > 0.0) {
477 s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j];
483 s = vcs_Hessian_diag_adj(irxn, s_old);
486 sprintf(anote,
"Normal calc: diag adjusted from %g "
487 "to %g due to act coeff", s_old, s);
490 m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
494 m_deltaMolNumSpecies[kspec] = tPhaseMoles;
500 for (
size_t j = 0; j < m_numComponents; ++j) {
501 double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
504 double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec];
505 if (negChangeComp > m_molNumSpecies_old[j]) {
506 if (m_molNumSpecies_old[j] > 0.0) {
508 sprintf(anote,
"Delta damped from %g "
509 "to %g due to component %lu (%10s) going neg", m_deltaMolNumSpecies[kspec],
510 -m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str());
512 m_deltaMolNumSpecies[kspec] = - 0.5 * m_molNumSpecies_old[j] / stoicC;
515 sprintf(anote,
"Delta damped from %g "
516 "to %g due to component %lu (%10s) zero", m_deltaMolNumSpecies[kspec],
517 -m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str());
519 m_deltaMolNumSpecies[kspec] = 0.0;
526 if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) {
528 sprintf(anote,
"Delta damped from %g "
529 "to %g due to %s going negative", m_deltaMolNumSpecies[kspec],
530 -m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str());
532 m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec];
537 vector<doublereal> fracDelta(Vphase->
nSpecies());
538 vector<doublereal> X_est(Vphase->
nSpecies());
541 double sumFrac = 0.0;
542 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
543 sumFrac += fracDelta[k];
545 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
546 X_est[k] = fracDelta[k] / sumFrac;
549 doublereal deltaMolNumPhase = tPhaseMoles;
550 doublereal damp = 1.0;
551 m_deltaGRxn_tmp = m_molNumSpecies_old;
552 double* molNumSpecies_tmp =
DATA_PTR(m_deltaGRxn_tmp);
555 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
557 double delmol = deltaMolNumPhase * X_est[k];
558 if (kspec >= m_numComponents) {
559 irxn = kspec - m_numComponents;
560 for (
size_t j = 0; j < m_numComponents; ++j) {
561 double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
564 molNumSpecies_tmp[j] += stoicC * delmol;
571 doublereal ratioComp = 0.0;
572 for (
size_t j = 0; j < m_numComponents; ++j) {
573 double deltaJ = m_molNumSpecies_old[j] - molNumSpecies_tmp[j];
574 if (molNumSpecies_tmp[j] < 0.0) {
577 double delta0 = m_molNumSpecies_old[j];
578 double dampj = delta0 / deltaJ * 0.9;
585 size_t jph = m_phaseID[j];
586 if ((jph != iphasePop) && (!m_SSPhase[j])) {
587 double fdeltaJ = fabs(deltaJ);
588 if (m_molNumSpecies_old[j] > 0.0) {
589 ratioComp =
std::max(ratioComp, fdeltaJ/ m_molNumSpecies_old[j]);
600 if (ratioComp > 1.0E-30) {
601 if (ratioComp < 0.001) {
602 damp = 0.001 / ratioComp;
607 if (damp <= 1.0E-6) {
611 for (
size_t k = 0; k < Vphase->
nSpecies(); k++) {
613 if (kspec < m_numComponents) {
616 m_deltaMolNumSpecies[kspec] = deltaMolNumPhase * X_est[k] * damp;
617 if (X_est[k] > 1.0E-3) {
631 double VCS_SOLVE::vcs_phaseStabilityTest(
const size_t iph)
637 size_t kspec, irxn, k, i, kc, kc_spec;
639 doublereal deltaGRxn;
642 bool doSuccessiveSubstitution =
true;
643 double funcPhaseStability;
644 vector<doublereal> X_est(Vphase->
nSpecies(), 0.0);
645 vector<doublereal> delFrac(Vphase->
nSpecies(), 0.0);
646 vector<doublereal> E_phi(Vphase->
nSpecies(), 0.0);
647 vector<doublereal> fracDelta_new(Vphase->
nSpecies(), 0.0);
648 vector<doublereal> fracDelta_old(Vphase->
nSpecies(), 0.0);
649 vector<doublereal> fracDelta_raw(Vphase->
nSpecies(), 0.0);
650 vector<size_t> creationGlobalRxnNumbers(Vphase->
nSpecies(),
npos);
653 vector<doublereal> m_feSpecies_Deficient(m_numComponents, 0.0);
654 doublereal damp = 1.0;
655 doublereal dampOld = 1.0;
656 doublereal normUpdate = 1.0;
657 doublereal normUpdateOld = 1.0;
658 doublereal sum = 0.0;
659 doublereal dirProd = 0.0;
660 doublereal dirProdOld = 0.0;
670 std::vector<size_t> componentList;
672 for (k = 0; k < Vphase->
nSpecies(); k++) {
674 if (kspec < m_numComponents) {
675 componentList.push_back(k);
679 for (k = 0; k < m_numComponents; k++) {
680 m_feSpecies_Deficient[k] = m_feSpecies_old[k];
682 normUpdate = 0.1 * vcs_l2norm(fracDelta_new);
685 if (doSuccessiveSubstitution) {
689 if (m_debug_print_lvl >= 2) {
690 plogf(
" --- vcs_phaseStabilityTest() called\n");
691 plogf(
" --- Its X_old[%2d] FracDel_old[%2d] deltaF[%2d] FracDel_new[%2d]"
692 " normUpdate damp FuncPhaseStability\n", KP, KP, KP, KP);
693 plogf(
" --------------------------------------------------------------"
694 "--------------------------------------------------------\n");
695 }
else if (m_debug_print_lvl == 1) {
696 plogf(
" --- vcs_phaseStabilityTest() called for phase %d\n", iph);
700 for (k = 0; k < Vphase->
nSpecies(); k++) {
701 if (fracDelta_new[k] < 1.0E-13) {
702 fracDelta_new[k] = 1.0E-13;
705 bool converged =
false;
706 for (
int its = 0; its < 200 && (!converged); its++) {
709 normUpdateOld = normUpdate;
710 fracDelta_old = fracDelta_new;
711 dirProdOld = dirProd;
715 for (i = 0; i < componentList.size(); i++) {
716 kc = componentList[i];
718 fracDelta_old[kc] = 0.0;
719 for (k = 0; k < Vphase->
nSpecies(); k++) {
721 if (kspec >= m_numComponents) {
722 irxn = kspec - m_numComponents;
723 fracDelta_old[kc] += m_stoichCoeffRxnMatrix[irxn][kc_spec] * fracDelta_old[k];
729 double sumFrac = 0.0;
730 for (k = 0; k < Vphase->
nSpecies(); k++) {
731 sumFrac += fracDelta_old[k];
734 if (sumFrac <= 0.0) {
737 double sum_Xcomp = 0.0;
738 for (k = 0; k < Vphase->
nSpecies(); k++) {
739 X_est[k] = fracDelta_old[k] / sumFrac;
741 if (kc_spec < m_numComponents) {
742 sum_Xcomp += X_est[k];
762 for (i = 0; i < componentList.size(); i++) {
763 kc = componentList[i];
766 m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
767 + log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]);
769 m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
774 for (i = 0; i < componentList.size(); i++) {
775 kc = componentList[i];
778 for (k = 0; k < Vphase->
nSpecies(); k++) {
780 if (kspec >= m_numComponents) {
781 irxn = kspec - m_numComponents;
783 m_deltaGRxn_Deficient[irxn] = m_deltaGRxn_old[irxn];
785 double* dtmp_ptr = m_stoichCoeffRxnMatrix[irxn];
786 if (dtmp_ptr[kc_spec] != 0.0) {
787 m_deltaGRxn_Deficient[irxn] +=
788 dtmp_ptr[kc_spec] * (m_feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]);
799 funcPhaseStability = sum_Xcomp - 1.0;
800 for (k = 0; k < Vphase->
nSpecies(); k++) {
802 if (kspec >= m_numComponents) {
803 irxn = kspec - m_numComponents;
804 deltaGRxn = m_deltaGRxn_Deficient[irxn];
805 if (deltaGRxn > 50.0) {
808 if (deltaGRxn < -50.0) {
811 E_phi[k] = std::exp(-deltaGRxn) / m_actCoeffSpecies_new[kspec];
813 funcPhaseStability += E_phi[k];
822 for (k = 0; k < Vphase->
nSpecies(); k++) {
824 double b = E_phi[k] / sum * (1.0 - sum_Xcomp);
825 if (kspec >= m_numComponents) {
826 irxn = kspec - m_numComponents;
827 fracDelta_raw[k] = b;
834 for (i = 0; i < componentList.size(); i++) {
835 kc = componentList[i];
837 fracDelta_raw[kc] = 0.0;
838 for (k = 0; k < Vphase->
nSpecies(); k++) {
840 if (kspec >= m_numComponents) {
841 irxn = kspec - m_numComponents;
842 fracDelta_raw[kc] += m_stoichCoeffRxnMatrix[irxn][kc_spec] * fracDelta_raw[k];
852 doublereal sumADel = 0.0;
853 for (k = 0; k < Vphase->
nSpecies(); k++) {
854 delFrac[k] = fracDelta_raw[k] - fracDelta_old[k];
855 sumADel += fabs(delFrac[k]);
857 normUpdate = vcs_l2norm(delFrac);
860 for (k = 0; k < Vphase->
nSpecies(); k++) {
861 dirProd += fracDelta_old[k] * delFrac[k];
863 bool crossedSign =
false;
864 if (dirProd * dirProdOld < 0.0) {
870 if (dampOld < 0.25) {
871 damp = 2.0 * dampOld;
874 if (normUpdate *1.5 > normUpdateOld) {
875 damp = 0.5 * dampOld;
876 }
else if (normUpdate *2.0 > normUpdateOld) {
877 damp = 0.8 * dampOld;
880 if (normUpdate > normUpdateOld * 2.0) {
881 damp = 0.6 * dampOld;
882 }
else if (normUpdate > normUpdateOld * 1.2) {
883 damp = 0.9 * dampOld;
887 for (k = 0; k < Vphase->
nSpecies(); k++) {
888 if (fabs(damp * delFrac[k]) > 0.3*fabs(fracDelta_old[k])) {
889 damp =
std::max(0.3*fabs(fracDelta_old[k]) / fabs(delFrac[k]),
890 1.0E-8/fabs(delFrac[k]));
892 if (delFrac[k] < 0.0) {
893 if (2.0 * damp * (-delFrac[k]) > fracDelta_old[k]) {
894 damp = fracDelta_old[k] / (2.0 * (-delFrac[k]));
897 if (delFrac[k] > 0.0) {
898 if (2.0 * damp * delFrac[k] > fracDelta_old[k]) {
899 damp = fracDelta_old[k] / (2.0 * delFrac[k]);
903 if (damp < 0.000001) {
907 for (k = 0; k < Vphase->
nSpecies(); k++) {
908 fracDelta_new[k] = fracDelta_old[k] + damp * (delFrac[k]);
912 if (m_debug_print_lvl >= 2) {
913 plogf(
" --- %3d %12g %12g %12g %12g %12g %12g %12g\n", its, X_est[KP], fracDelta_old[KP],
914 delFrac[KP], fracDelta_new[KP], normUpdate, damp, funcPhaseStability);
918 if (normUpdate < 1.0E-5) {
932 printf(
"not done yet\n");
936 if (m_debug_print_lvl >= 2) {
937 plogf(
" ------------------------------------------------------------"
938 "-------------------------------------------------------------\n");
939 }
else if (m_debug_print_lvl == 1) {
940 if (funcPhaseStability > 0.0) {
941 plogf(
" --- phase %d with func = %g is to be born\n", iph, funcPhaseStability);
943 plogf(
" --- phase %d with func = %g stays dead\n", iph, funcPhaseStability);
947 return funcPhaseStability;