diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp index f4b0a3b7c..aa7eac6e6 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp @@ -1559,18 +1559,14 @@ namespace VCSnonideal { * states. */ m_mix->uploadMoleFractionsFromPhases(); - int kGlob = 0; - for (int ip = 0; ip < m_vprob->NPhase; ip++) { - double phaseMole = 0.0; - Cantera::ThermoPhase &tref = m_mix->phase(ip); - int nspPhase = tref.nSpecies(); - for (int k = 0; k < nspPhase; k++, kGlob++) { - phaseMole += m_vprob->w[kGlob]; - } - //phaseMole *= 1.0E-3; - m_mix->setPhaseMoles(ip, phaseMole); - } - + for (int i = 0; i < m_vprob->nspecies; i++) { + plogf("%d %15.3e\n", m_vprob->m_gibbsSpecies[i]); + } + m_mix->getChemPotentials(DATA_PTR(m_vprob->m_gibbsSpecies)); + for (int i = 0; i < m_vprob->nspecies; i++) { + plogf("%d %15.3e\n", m_vprob->m_gibbsSpecies[i]); + } + double te = tickTock.secondsWC(); if (printLvl > 0) { plogf("\n Results from vcs_PS:\n"); @@ -1615,11 +1611,11 @@ namespace VCSnonideal { if (m_vprob->w[i] <= 0.0) { int iph = m_vprob->PhaseID[i]; vcs_VolPhase *VPhase = m_vprob->VPhaseList[iph]; - if (VPhase->nSpecies() > 1) { - plogf(" -1.000e+300\n"); - } else { + //if (VPhase->nSpecies() > 1) { + // plogf(" -1.000e+300\n"); + //} else { plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]); - } + //} } else { plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]); } diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index dc790bb28..3f52c273c 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -896,7 +896,7 @@ namespace VCSnonideal { retn = VCS_PUB_BAD; } - if (vPhase->PhaseName == pub_phase_ptr->PhaseName) { + if (vPhase->PhaseName != pub_phase_ptr->PhaseName) { plogf("%sPhaseName value have changed:%s %s\n", yo.c_str(), vPhase->PhaseName.c_str(), pub_phase_ptr->PhaseName.c_str()); diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index cc47b2061..954bfebda 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -1257,16 +1257,16 @@ private: * the mole numbers of the component species. Therefore the following * approximation is valid for a small component of an ideal phase: * - * 0 = m_deltaGrxn_old(I) + log(molNum_new(I)/molNum_old(I)) + * 0 = m_deltaGRxn_old(I) + log(molNum_new(I)/molNum_old(I)) * - * m_deltaGrxn_old contains the contribution from + * m_deltaGRxn_old contains the contribution from * * m_feSpecies_old(I) = * m_SSfeSpecies(I) + * log(ActCoeff[i] * molNum_old(I) / m_tPhaseMoles_old(iph)) * Thus, * - * molNum_new(I)= molNum_old(I) * EXP(-m_deltaGrxn_old(I)) + * molNum_new(I)= molNum_old(I) * EXP(-m_deltaGRxn_old(I)) * * Most of this section is mainly restricting the update to reasonable * values. diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 2772202f7..1f4eaf3e3 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -2015,16 +2015,16 @@ namespace VCSnonideal { * the mole numbers of the component species. Therefore the following * approximation is valid for a small component of an ideal phase: * - * 0 = m_deltaGrxn_old(I) + log(molNum_new(I)/molNum_old(I)) + * 0 = m_deltaGRxn_old(I) + log(molNum_new(I)/molNum_old(I)) * - * m_deltaGrxn_old contains the contribution from + * m_deltaGRxn_old contains the contribution from * * m_feSpecies_old(I) = * m_SSfeSpecies(I) + * log(ActCoeff[i] * molNum_old(I) / m_tPhaseMoles_old(iph)) * Thus, * - * molNum_new(I)= molNum_old(I) * EXP(-m_deltaGrxn_old(I)) + * molNum_new(I)= molNum_old(I) * EXP(-m_deltaGRxn_old(I)) * * Most of this section is mainly restricting the update to reasonable * values. diff --git a/Cantera/src/equil/vcs_solve_phaseStability.cpp b/Cantera/src/equil/vcs_solve_phaseStability.cpp index 5244414c8..f3e6b35de 100644 --- a/Cantera/src/equil/vcs_solve_phaseStability.cpp +++ b/Cantera/src/equil/vcs_solve_phaseStability.cpp @@ -41,7 +41,16 @@ namespace VCSnonideal { int ifunc = 0; int iStab = 0; - + /* + * This function is called to create the private data + * using the public data. + */ + int nspecies0 = vprob->nspecies + 10; + int nelements0 = vprob->ne; + int nphase0 = vprob->NPhase; + + vcs_initSizes(nspecies0, nelements0, nphase0); + if (ifunc < 0 || ifunc > 2) { plogf("vcs: Unrecognized value of ifunc, %d: bailing!\n", @@ -131,8 +140,7 @@ namespace VCSnonideal { * (all information concerning Temperature and Pressure has already * been derived. The free energies are now in dimensionless form.) */ - double funcVal; - iStab = vcs_solve_phaseStability(iphase, ifunc, funcVal, printLvl); + iStab = vcs_solve_phaseStability(iphase, ifunc, feStable, printLvl); /* @@ -140,6 +148,23 @@ namespace VCSnonideal { * the reverse of vcs_nondim to add back units. */ vcs_redim_TP(); + + /* + vcs_VolPhase *Vphase = m_VolPhaseList[iphase]; + + std::vector mfPop = Vphase->moleFractions(); + int nsp = Vphase->nSpecies(); + + vcs_VolPhase *VPphase = vprob->VPhaseList[iphase]; + int kstart = Vphase->spGlobalIndexVCS(0); + for (int k = 0; k < nsp; k++) { + vprob->mf[kstart + k] = mfPop[k]; + } + VPphase->setMoleFractionsState(Vphase->totalMoles(), + VCS_DATA_PTR(Vphase->moleFractions()), + VCS_STATECALC_TMP); + */ + vcs_prob_update(vprob); /* * Return the convergence success flag. */ @@ -197,6 +222,11 @@ namespace VCSnonideal { retn = vcs_basopt(FALSE, VCS_DATA_PTR(aw), VCS_DATA_PTR(sa), VCS_DATA_PTR(sm), VCS_DATA_PTR(ss), test, &usedZeroedSpecies); + + + vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc); + vcs_deltag(0, true, VCS_STATECALC_OLD); + phasePopPhaseIDs.clear(); iphasePop = vcs_popPhaseID(phasePopPhaseIDs); funcVal = vcs_phaseStabilityTest(iph);