diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index 51d4515a0..664544779 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -507,6 +507,10 @@ namespace VCSnonideal { const std::vector & vcs_VolPhase::moleFractions() const { return Xmol_; } + + double vcs_VolPhase::moleFraction(int k) const { + return Xmol_[k]; + } /***************************************************************************/ // Set the moles and/or mole fractions within the phase diff --git a/Cantera/src/equil/vcs_VolPhase.h b/Cantera/src/equil/vcs_VolPhase.h index 03715db1e..15cd25d3c 100644 --- a/Cantera/src/equil/vcs_VolPhase.h +++ b/Cantera/src/equil/vcs_VolPhase.h @@ -419,6 +419,8 @@ namespace VCSnonideal { //! object. const std::vector & moleFractions() const; + double moleFraction(int klocal) const; + //! Sets the creationMoleNum's within the phase object /*! * @param F_k Pointer to a vector of n_k's diff --git a/Cantera/src/equil/vcs_phaseStability.cpp b/Cantera/src/equil/vcs_phaseStability.cpp index 997bbf830..6d323da8c 100644 --- a/Cantera/src/equil/vcs_phaseStability.cpp +++ b/Cantera/src/equil/vcs_phaseStability.cpp @@ -671,7 +671,7 @@ namespace VCSnonideal { vector fracDelta_old(Vphase->nSpecies(), 0.0); vector fracDelta_raw(Vphase->nSpecies(), 0.0); vector creationGlobalRxnNumbers(Vphase->nSpecies(), -1); - + vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_Deficient), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc); vector m_feSpecies_Deficient(m_numComponents, 0.0); doublereal damp = 1.0; @@ -779,18 +779,16 @@ namespace VCSnonideal { Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, VCS_DATA_PTR(m_actCoeffSpecies_new)); /* - * first Calculate altered chemical potentials for component species + * First calculate altered chemical potentials for component species * belonging to this phase. */ for (i = 0; i < (int) componentList.size(); i++) { kc = componentList[i]; kc_spec = Vphase->spGlobalIndexVCS(kc); if ( X_est[kc] > VCS_DELETE_MINORSPECIES_CUTOFF) { - m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec] - + log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]); + m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec] + log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]); } else { - m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec] - + log(m_actCoeffSpecies_new[kc_spec] * VCS_DELETE_MINORSPECIES_CUTOFF); + m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec] + log(m_actCoeffSpecies_new[kc_spec] * VCS_DELETE_MINORSPECIES_CUTOFF); } } @@ -807,8 +805,7 @@ namespace VCSnonideal { } double *dtmp_ptr = m_stoichCoeffRxnMatrix[irxn]; if (dtmp_ptr[kc_spec] != 0.0) { - m_deltaGRxn_Deficient[irxn] += - dtmp_ptr[kc_spec] * (m_feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]); + m_deltaGRxn_Deficient[irxn] += dtmp_ptr[kc_spec] * (m_feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]); } } diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 954bfebda..daf952d26 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -495,13 +495,20 @@ public: */ void vcs_dfe(const int stateCalc, const int ll, const int lbot, const int ltop); + //! Print out a table of chemical potentials + /*! + * @param vcsState Determines where to get the mole numbers from. + * - VCS_STATECALC_OLD -> from m_molNumSpecies_old + * - VCS_STATECALC_NEW -> from m_molNumSpecies_new + */ + void vcs_printSpeciesChemPot(const int stateCalc) const; + //! This routine uploads the state of the system into all of the //! vcs_VolumePhase objects in the current problem. /*! * @param vcsState Determines where to get the mole numbers from. * - VCS_STATECALC_OLD -> from m_molNumSpecies_old * - VCS_STATECALC_NEW -> from m_molNumSpecies_new - * */ void vcs_updateVP(const int stateCalc); @@ -621,6 +628,8 @@ public: void vcs_deltag(const int l, const bool doDeleted, const int vcsState, const bool alterZeroedPhases = true); + void vcs_printDeltaG(const int stateCalc); + //! Calculate deltag of formation for all species in a single phase. /*! * Calculate deltag of formation for all species in a single @@ -1618,8 +1627,10 @@ public: //! Last deltag[irxn] from the previous step std::vector m_deltaGRxn_old; - //! Last deltag[irxn] from the previous step with additions for - //! possible births of zeroed phases. + //! Last deltag[irxn] from the previous step with additions for possible births of zeroed phases for component species + /*! + * + */ std::vector m_deltaGRxn_Deficient; //! Temporary vector of Rxn DeltaG's @@ -1684,10 +1695,10 @@ public: std::vector m_tPhaseMoles_new; //! Temporary vector of length NPhase - std::vector m_TmpPhase; + mutable std::vector m_TmpPhase; //! Temporary vector of length NPhase - std::vector m_TmpPhase2; + mutable std::vector m_TmpPhase2; //! Change in the total moles in each phase /*! diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 1f4eaf3e3..7dfc0d821 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -4371,10 +4371,10 @@ namespace VCSnonideal { feSpecies[kspec] = m_SSfeSpecies[kspec] + log(actCoeff_ptr[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF) - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec] - + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; ; + + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; } else { feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec] - + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; ; + + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; } } else { feSpecies[kspec] = m_SSfeSpecies[kspec] @@ -4438,7 +4438,93 @@ namespace VCSnonideal { } } } + + } + //==================================================================================================================== + // Print out a table of chemical potentials + /* + * @param vcsState Determines where to get the mole numbers from. + * - VCS_STATECALC_OLD -> from m_molNumSpecies_old + * - VCS_STATECALC_NEW -> from m_molNumSpecies_new + */ + void VCS_SOLVE::vcs_printSpeciesChemPot(const int stateCalc) const { + double mfValue = 1.0; + bool zeroedPhase = false; + int kspec; + + const double * molNum = VCS_DATA_PTR(m_molNumSpecies_old); + const double * tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_old); + const double * actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_old); + if (stateCalc == VCS_STATECALC_NEW) { + tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new); + actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new); + molNum = VCS_DATA_PTR(m_molNumSpecies_new); + } + + double * tMoles = VCS_DATA_PTR(m_TmpPhase); + const double *tPhInertMoles = VCS_DATA_PTR(TPhInertMoles); + for (int iph = 0; iph < m_numPhases; iph++) { + tMoles[iph] = tPhInertMoles[iph]; + } + for (kspec = 0; kspec < m_numSpeciesTot; kspec++) { + if(m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + int iph = m_phaseID[kspec]; + tMoles[iph] += molNum[kspec]; + } + } + + double RT = m_temperature * Cantera::GasConstant; + printf(" --- CHEMICAL POT TABLE (J/kmol) Name PhID MolFR ChemoSS " + " logMF Gamma Elect extra ElectrChem\n"); + printf(" "); + vcs_print_line("-", 132); + + for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) { + mfValue = 1.0; + int iphase = m_phaseID[kspec]; + const vcs_VolPhase * Vphase = m_VolPhaseList[iphase]; + if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) || + (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) || + (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDSS) ) { + zeroedPhase = true; + } else { + zeroedPhase = false; + } + if (tMoles[iphase] > 0.0) { + if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) { + mfValue = VCS_DELETE_MINORSPECIES_CUTOFF / tMoles[iphase]; + } else { + mfValue = molNum[kspec]/tMoles[iphase]; + } + } else { + int klocal = m_speciesLocalPhaseIndex[kspec]; + mfValue = Vphase->moleFraction(klocal); + } + double volts = Vphase->electricPotential(); + double elect = m_chargeSpecies[kspec] * m_Faraday_dim * volts; + double comb = - m_lnMnaughtSpecies[kspec]; + double total = (m_SSfeSpecies[kspec] + log(mfValue) + elect + log(actCoeff_ptr[kspec]) + comb); + + if (zeroedPhase) { + printf(" --- ** zp *** "); + } else { + printf(" --- "); + } + printf("%-24.24s", m_speciesName[kspec].c_str()); + printf(" %-3d", iphase); + printf(" % -12.4e", mfValue); + printf(" % -12.4e", m_SSfeSpecies[kspec] * RT); + printf(" % -12.4e", log(mfValue) * RT); + printf(" % -12.4e", log(actCoeff_ptr[kspec]) * RT); + printf(" % -12.4e", elect * RT); + printf(" % -12.4e", comb * RT); + printf(" % -12.4e\n", total *RT); + } + printf(" "); + vcs_print_line("-", 132); + } + /*****************************************************************************/ #ifdef DEBUG_MODE @@ -4990,8 +5076,125 @@ namespace VCSnonideal { } #endif } - /*****************************************************************************/ + //==================================================================================================================== + void VCS_SOLVE::vcs_printDeltaG( const int stateCalc) { + int j; + double * deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_old); + double * feSpecies = VCS_DATA_PTR(m_feSpecies_old); + double * molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_old); + const double * tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_old); + const double * actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_old); + if (stateCalc == VCS_STATECALC_NEW) { + deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_new); + feSpecies = VCS_DATA_PTR(m_feSpecies_new); + molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_new); + actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new); + tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new); + } + double RT = m_temperature * Cantera::GasConstant; + bool zeroedPhase = false; + if (m_debug_print_lvl >= 2) { + plogf(" --- DELTA_G TABLE Components:"); + for (j = 0; j < m_numComponents; j++) { + plogf(" %3d ", j); + } + plogf("\n --- Components Moles:"); + for (j = 0; j < m_numComponents; j++) { + plogf("%10.3g", m_molNumSpecies_old[j]); + } + plogf("\n --- NonComponent| Moles | "); + for (j = 0; j < m_numComponents; j++) { + plogf("%-10.10s", m_speciesName[j].c_str()); + } + //plogf("| m_scSize"); + plogf("\n"); + for (int i = 0; i < m_numRxnTot; i++) { + plogf(" --- %3d ", m_indexRxnToSpecies[i]); + plogf("%-10.10s", m_speciesName[m_indexRxnToSpecies[i]].c_str()); + plogf("|%10.3g|", m_molNumSpecies_old[m_indexRxnToSpecies[i]]); + for (j = 0; j < m_numComponents; j++) { + plogf(" %6.2f", m_stoichCoeffRxnMatrix[i][j]); + } + //plogf(" | %6.2f", m_scSize[i]); + plogf("\n"); + } + plogf(" "); for(int i=0; i<77; i++) plogf("-"); plogf("\n"); + } + printf(" --- DeltaG Table (J/kmol) Name PhID MoleNum MolFR " + " ElectrChemStar ElectrChem DeltaGStar DeltaG(Pred) Stability\n"); + printf(" "); + vcs_print_line("-", 132); + + for (int kspec = 0; kspec < m_numSpeciesTot; kspec++) { + + int irxn = kspec - m_numComponents; + + double mfValue = 1.0; + int iphase = m_phaseID[kspec]; + const vcs_VolPhase * Vphase = m_VolPhaseList[iphase]; + if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) || + (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) || + (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDSS) ) { + zeroedPhase = true; + } else { + zeroedPhase = false; + } + if (tPhMoles_ptr[iphase] > 0.0) { + if (molNumSpecies[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) { + mfValue = VCS_DELETE_MINORSPECIES_CUTOFF / tPhMoles_ptr[iphase]; + } else { + mfValue = molNumSpecies[kspec] / tPhMoles_ptr[iphase]; + } + } else { + int klocal = m_speciesLocalPhaseIndex[kspec]; + mfValue = Vphase->moleFraction(klocal); + } + if (zeroedPhase) { + printf(" --- ** zp *** "); + } else { + printf(" --- "); + } + double feFull = feSpecies[kspec]; + if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) || + (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) ) { + feFull += log(actCoeff_ptr[kspec]) + log(mfValue); + } + printf("%-24.24s", m_speciesName[kspec].c_str()); + printf(" %-3d", iphase); + printf(" % -12.4e", molNumSpecies[kspec]); + printf(" % -12.4e", mfValue); + printf(" % -12.4e", feSpecies[kspec] * RT); + printf(" % -12.4e", feFull * RT); + if (irxn >= 0) { + printf(" % -12.4e", deltaGRxn[irxn] * RT); + printf(" % -12.4e", (deltaGRxn[irxn] + feFull - feSpecies[kspec]) * RT); + + if (deltaGRxn[irxn] < 0.0) { + if ( molNumSpecies[kspec] > 0.0) { + printf(" growing"); + } else { + printf(" stable"); + } + } else if (deltaGRxn[irxn] > 0.0) { + if ( molNumSpecies[kspec] > 0.0) { + printf(" shrinking"); + } else { + printf(" unstable"); + } + } else { + printf(" balanced"); + } + } + + printf(" \n"); + } + + printf(" "); + vcs_print_line("-", 132); + + } + //==================================================================================================================== // Calculate deltag of formation for all species in a single phase. /* * Calculate deltag of formation for all species in a single diff --git a/Cantera/src/equil/vcs_solve_phaseStability.cpp b/Cantera/src/equil/vcs_solve_phaseStability.cpp index f3e6b35de..73e851416 100644 --- a/Cantera/src/equil/vcs_solve_phaseStability.cpp +++ b/Cantera/src/equil/vcs_solve_phaseStability.cpp @@ -203,7 +203,7 @@ namespace VCSnonideal { */ int VCS_SOLVE::vcs_solve_phaseStability(const int iph, const int ifunc, double &funcVal, - int printLv) { + int printLvl) { int retn = 0; double test = -1.0E-10; int usedZeroedSpecies; @@ -222,11 +222,18 @@ 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_evaluate_speciesType(); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc); + if (printLvl > 3) { + vcs_printSpeciesChemPot(VCS_STATECALC_OLD); + } vcs_deltag(0, true, VCS_STATECALC_OLD); + if (printLvl > 3) { + vcs_printDeltaG(VCS_STATECALC_OLD); + } + vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_Deficient), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc); phasePopPhaseIDs.clear(); iphasePop = vcs_popPhaseID(phasePopPhaseIDs); funcVal = vcs_phaseStabilityTest(iph);