diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index 30eee5b42..9e14c546f 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -146,7 +146,11 @@ namespace VCSnonideal { if (&b != this) { int old_num = NVolSpecies; - m_owningSolverObject = b.m_owningSolverObject; + // Note: we comment this out for the assignment operator + // specifically, because it isn't true for the assignment + // operator but is true for a copy constructor + // m_owningSolverObject = b.m_owningSolverObject; + VP_ID = b.VP_ID; Domain_ID = b.Domain_ID; SingleSpecies = b.SingleSpecies; @@ -688,7 +692,9 @@ namespace VCSnonideal { * in all of the phases in a VCS problem. Only the * entries for the current phase are filled in. */ - void vcs_VolPhase::sendToVCS_ActCoeff(double * const AC) const { + void vcs_VolPhase::sendToVCS_ActCoeff(const int stateCalc, + double * const AC) { + updateFromVCS_MoleNumbers(stateCalc); if (!m_UpToDate_AC) { evaluateActCoeff(); } diff --git a/Cantera/src/equil/vcs_VolPhase.h b/Cantera/src/equil/vcs_VolPhase.h index 29d6a8fb0..7371e3135 100644 --- a/Cantera/src/equil/vcs_VolPhase.h +++ b/Cantera/src/equil/vcs_VolPhase.h @@ -226,7 +226,7 @@ namespace VCSnonideal { * in all of the phases in a VCS problem. Only the * entries for the current phase are filled in. */ - void sendToVCS_ActCoeff(double * const AC) const; + void sendToVCS_ActCoeff(const int stateCalc, double * const AC); //! set the electric potential of the phase /*! diff --git a/Cantera/src/equil/vcs_report.cpp b/Cantera/src/equil/vcs_report.cpp index 96edd0725..6c8ff892f 100644 --- a/Cantera/src/equil/vcs_report.cpp +++ b/Cantera/src/equil/vcs_report.cpp @@ -203,22 +203,21 @@ namespace VCSnonideal { } plogf(" | DG/RT Rxn |\n"); print_line("-", m_numComponents*10 + 45); - for (i = 0; i < m_numRxnTot; i++) { - int kspec = m_indexRxnToSpecies[i]; + for (int irxn = 0; irxn < m_numRxnTot; irxn++) { + int kspec = m_indexRxnToSpecies[irxn]; plogf(" %3d ", kspec); plogf("%-10.10s", m_speciesName[kspec].c_str()); plogf("|%10.3g |", m_molNumSpecies_old[kspec]); for (j = 0; j < m_numComponents; j++) { - plogf(" %6.2f", m_stoichCoeffRxnMatrix[i][j]); + plogf(" %6.2f", m_stoichCoeffRxnMatrix[irxn][j]); } - // Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve - plogf(" |%10.3g |", m_deltaGRxn_new[kspec]); + plogf(" |%10.3g |", m_deltaGRxn_new[irxn]); plogf("\n"); } print_line("-", m_numComponents*10 + 45); plogf("\n"); - /* + /* * ------------------ TABLE OF PHASE INFORMATION --------------------- */ std::vector gaPhase(m_numElemConstraints, 0.0); @@ -262,7 +261,7 @@ namespace VCSnonideal { plogf(" %10.3g", gaPhase[j]); gaTPhase[j] += gaPhase[j]; } - gibbsPhase = vcs_GibbsPhase(iphase, VCS_DATA_PTR(m_molNumSpecies_old), + gibbsPhase = vcs_GibbsPhase(iphase, VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old)); gibbsTotal += gibbsPhase; plogf(" | %18.11E |\n", gibbsPhase); @@ -313,8 +312,12 @@ namespace VCSnonideal { vcs_printChemPotUnits(m_VCS_UnitsFormat); plogf(")\n"); plogf(" Name TKMoles StandStateChemPot " - " ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)\n"); - print_line("-", 115); + " ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)"); +#ifdef DEBUG_MODE + plogf("| (MolNum ChemPot)|"); +#endif + plogf("\n"); + print_line("-", 147); for (i = 0; i < nspecies; ++i) { l = sortindex[i]; int pid = m_phaseID[l]; @@ -346,12 +349,22 @@ namespace VCSnonideal { exit(-1); } plogf(" %12.4E |", m_feSpecies_old[l]); - if (m_lnMnaughtSpecies[l] != 0.0) { - plogf(" (%14.7E)", - m_lnMnaughtSpecies[l]); + if( m_lnMnaughtSpecies[l] != 0.0) { + plogf("(%11.5E)", - m_lnMnaughtSpecies[l]); + } else { + plogf(" "); } + +#ifdef DEBUG_MODE + plogf("| %20.13E |", m_feSpecies_old[l] * m_molNumSpecies_old[l]); +#endif plogf("\n"); } - print_line("-", 115); +#ifdef DEBUG_MODE + for (i = 0; i < 125; i++) plogf(" "); + plogf("%20.13E\n", g); +#endif + print_line("-", 147); /* * ------------- TABLE OF SOLUTION COUNTERS -------------------------- diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index 8e33bc259..317df3602 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -372,13 +372,14 @@ namespace VCSnonideal { * Note, this is a dangerous routine that leaves the underlying objects in * an unknown state. */ - double VCS_SOLVE::deltaG_Recalc_Rxn(const int irxn, const double *const molNum, + double VCS_SOLVE::deltaG_Recalc_Rxn(const int stateCalc, + const int irxn, const double *const molNum, double * const ac, double * const mu_i) { int kspec = irxn + m_numComponents; int *pp_ptr = m_phaseParticipation[irxn]; for (int iphase = 0; iphase < m_numPhases; iphase++) { if (pp_ptr[iphase]) { - vcs_chemPotPhase(iphase, molNum, ac, mu_i); + vcs_chemPotPhase(stateCalc, iphase, molNum, ac, mu_i); } } double deltaG = mu_i[kspec]; @@ -424,7 +425,9 @@ namespace VCSnonideal { /* * Calculate the deltaG value at the dx = 0.0 point */ - double deltaGOrig = deltaG_Recalc_Rxn(irxn, molNumBase, acBase, + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); + double deltaGOrig = deltaG_Recalc_Rxn(VCS_STATECALC_OLD, + irxn, molNumBase, acBase, VCS_DATA_PTR(m_feSpecies_old)); double forig = fabs(deltaGOrig) + 1.0E-15; if (deltaGOrig > 0.0) { @@ -461,8 +464,10 @@ namespace VCSnonideal { m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx_orig; molSum += molNumBase[k]; } + vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); - double deltaG1 = deltaG_Recalc_Rxn(irxn, VCS_DATA_PTR(m_molNumSpecies_new), + double deltaG1 = deltaG_Recalc_Rxn(VCS_STATECALC_NEW, + irxn, VCS_DATA_PTR(m_molNumSpecies_new), ac, VCS_DATA_PTR(m_feSpecies_new)); /* @@ -500,7 +505,9 @@ namespace VCSnonideal { for (k = 0; k < m_numComponents; k++) { m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx; } - double deltaG = deltaG_Recalc_Rxn(irxn, VCS_DATA_PTR(m_molNumSpecies_new), + vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); + double deltaG = deltaG_Recalc_Rxn(VCS_STATECALC_NEW, + irxn, VCS_DATA_PTR(m_molNumSpecies_new), ac, VCS_DATA_PTR(m_feSpecies_new)); /* * If deltaG hasn't switched signs when going the full distance @@ -524,6 +531,7 @@ namespace VCSnonideal { } finalize: + vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); if (its >= MAXITS) { #ifdef DEBUG_MODE sprintf(ANOTE,"Rxn reduced to zero step size from %g to %g (MAXITS)", diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 0ddf84cb4..84da9fa94 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -341,7 +341,7 @@ public: * (VCS species order) * */ - void vcs_chemPotPhase(const int iph, const double *const molNum, + void vcs_chemPotPhase(const int stateCalc, const int iph, const double *const molNum, double * const ac, double * const mu_i, const bool do_deleted = false); @@ -1301,7 +1301,8 @@ private: * * @return Returns the dimensionless deltaG of the reaction */ - double deltaG_Recalc_Rxn(const int irxn, const double *const molNum, + double deltaG_Recalc_Rxn(const int stateCalc, + const int irxn, const double *const molNum, double * const ac, double * const mu_i); //! Delete memory that isn't just resizeable STL containers @@ -1494,16 +1495,13 @@ public: //std::vector wt; std::vector m_molNumSpecies_new; - //! Delta G(I) for the noncomponent species in the mechanism. + //! Delta G(irxn) for the noncomponent species in the mechanism. /*! - * Computed by the subroutine DELTAG. DG is the free - * energy change for the reaction which - * forms species K from the - * component species. This vector has length - * equal to the number of noncomponent - * species in the mechanism. It starts with - * the first current noncomponent species - * in the mechanism. + * Computed by the subroutine deltaG. m_deltaGRxn is the free + * energy change for the reaction which forms species K from the + * component species. This vector has length equal to the number + * of noncomponent species in the mechanism. It starts with + * the first current noncomponent species in the mechanism. */ std::vector m_deltaGRxn_new; diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index a3ae00397..3f545f748 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -926,11 +926,28 @@ namespace VCSnonideal { (doPhaseDeleteIph == -1) && (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) { double dx_old = dx; +#ifdef DEBUG_HKM_NOT + std::vector feSpecies_tmp(m_numSpeciesTot); + feSpecies_tmp = m_feSpecies_old; + std::vector molNumSpecies_tmp(m_numSpeciesTot); + molNumSpecies_tmp = m_molNumSpecies_old; +#endif #ifdef DEBUG_MODE dx = vcs_line_search(irxn, dx_old, ANOTE); #else dx = vcs_line_search(irxn, dx_old); #endif +#ifdef DEBUG_HKM_NOT + for (int kk = 0; kk < m_numSpeciesTot; kk++) { + if (feSpecies_tmp[kk] != m_feSpecies_old[kk]) { + printf("we are here\n"); + } + if (molNumSpecies_tmp[kk] != m_molNumSpecies_old[kk]) { + printf("we are here\n"); + } + } +#endif + vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); } m_deltaMolNumSpecies[kspec] = dx; @@ -1942,14 +1959,14 @@ namespace VCSnonideal { * Store the final Delta G values for each non-component species * in the species slot rather than the reaction slot */ - kspec = m_numSpeciesTot; - i = m_numRxnTot; - for (irxn = 0; irxn < m_numRxnTot; ++irxn) { - --kspec; - --i; - m_deltaGRxn_new[kspec] = m_deltaGRxn_new[i]; - } - vcs_dzero(VCS_DATA_PTR(m_deltaGRxn_new), m_numComponents); + // kspec = m_numSpeciesTot; + // i = m_numRxnTot; + //for (irxn = 0; irxn < m_numRxnTot; ++irxn) { + // --kspec; + // --i; + // m_deltaGRxn_new[kspec] = m_deltaGRxn_new[i]; + //} + // vcs_dzero(VCS_DATA_PTR(m_deltaGRxn_new), m_numComponents); /* * Evaluate the final mole fractions * storring them in wt[] @@ -2182,13 +2199,18 @@ namespace VCSnonideal { m_molNumSpecies_old[kspec] += dx; int iph = m_phaseID[kspec]; m_tPhaseMoles_old[iph] += dx; + vcs_setFlagsVolPhase(iph, false, VCS_STATECALC_OLD); + for (j = 0; j < m_numComponents; ++j) { - iph = m_phaseID[j]; tmp = sc_irxn[j] * dx; - m_molNumSpecies_old[j] += tmp; - m_tPhaseMoles_old[iph] += tmp; - if (m_molNumSpecies_old[j] < 0.0) { - m_molNumSpecies_old[j] = 0.0; + if (tmp != 0.0) { + iph = m_phaseID[j]; + m_molNumSpecies_old[j] += tmp; + m_tPhaseMoles_old[iph] += tmp; + vcs_setFlagsVolPhase(iph, false, VCS_STATECALC_OLD); + if (m_molNumSpecies_old[j] < 0.0) { + m_molNumSpecies_old[j] = 0.0; + } } } } @@ -2936,7 +2958,7 @@ namespace VCSnonideal { * only step is being carried out, then we don't need to * update the minor noncomponents. */ - vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); + vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesRdc); /* @@ -4388,14 +4410,14 @@ namespace VCSnonideal { * (VCS species order) * */ - void VCS_SOLVE::vcs_chemPotPhase(const int iph, const double *const molNum, + void VCS_SOLVE::vcs_chemPotPhase(const int stateCalc, + const int iph, const double *const molNum, double * const ac, double * const mu_i, const bool do_deleted) { vcs_VolPhase *Vphase = m_VolPhaseList[iph]; int nkk = Vphase->NVolSpecies; int k, kspec; - int stateCalc = VCS_STATECALC_OLD; #ifdef DEBUG_MODE //if (m_debug_print_lvl >= 2) { @@ -4414,7 +4436,7 @@ namespace VCSnonideal { } Vphase->setMolesFromVCS(stateCalc, molNum); - Vphase->sendToVCS_ActCoeff(ac); + Vphase->sendToVCS_ActCoeff(stateCalc, ac); double phi = Vphase->electricPotential(); double Faraday_phi = m_Faraday_dim * phi; @@ -4695,7 +4717,7 @@ namespace VCSnonideal { Vphase = m_VolPhaseList[iphase]; if (!Vphase->SingleSpecies) { Vphase->setMolesFromVCS(stateCalc, molNum); - Vphase->sendToVCS_ActCoeff(VCS_DATA_PTR(actCoeff_ptr)); + Vphase->sendToVCS_ActCoeff(stateCalc, VCS_DATA_PTR(actCoeff_ptr)); } m_phasePhi[iphase] = Vphase->electricPotential(); m_phaseACAreCurrent[iphase] = 1;