diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp index 1b9b474e3..849628646 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp @@ -1371,7 +1371,7 @@ namespace VCSnonideal { */ if (tMoles > 0.0) { for (k = 0; k < nSpPhase; k++) { - int kTa = VolPhase->IndSpecies[k]; + int kTa = VolPhase->spGlobalIndexVCS(k); vprob->mf[kTa] = vprob->w[kTa] / tMoles; } } else { @@ -1380,7 +1380,7 @@ namespace VCSnonideal { * But, will start with this. */ for (k = 0; k < nSpPhase; k++) { - int kTa = VolPhase->IndSpecies[k]; + int kTa = VolPhase->spGlobalIndexVCS(k); vprob->mf[kTa]= 1.0 / (double) nSpPhase; } } @@ -1507,7 +1507,8 @@ namespace VCSnonideal { kT++; } if (volPhase->phiVarIndex() >= 0) { - int kglob = volPhase->IndSpecies[volPhase->phiVarIndex()]; + int kphi = volPhase->phiVarIndex(); + int kglob = volPhase->spGlobalIndexVCS(kphi); vprob->w[kglob] = tPhase->electricPotential(); } volPhase->setMolesFromVCS(VCS_STATECALC_OLD, VCS_DATA_PTR(vprob->w)); diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index 8202bd214..6849fae63 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -1239,6 +1239,33 @@ namespace VCSnonideal { #endif m_existence = existence; } - /**********************************************************************/ + /**********************************************************************/ + + // Return the Global VCS index of the kth species in the phase + /* + * @param spIndex local species index (0 to the number of species + * in the phase) + * + * @return Returns the VCS_SOLVE species index of the that species + * This changes as rearrangements are carried out. + */ + int vcs_VolPhase::spGlobalIndexVCS(const int spIndex) const { + return IndSpecies[spIndex]; + } + /**********************************************************************/ + + //! set the Global VCS index of the kth species in the phase + /*! + * @param spIndex local species index (0 to the number of species + * in the phase) + * + * @return Returns the VCS_SOLVE species index of the that species + * This changes as rearrangements are carried out. + */ + void vcs_VolPhase::setSpGlobalIndexVCS(const int spIndex, + const int spGlobalIndex) { + IndSpecies[spIndex] = spGlobalIndex; + } + /**********************************************************************/ } diff --git a/Cantera/src/equil/vcs_VolPhase.h b/Cantera/src/equil/vcs_VolPhase.h index 65495d16f..8171fd489 100644 --- a/Cantera/src/equil/vcs_VolPhase.h +++ b/Cantera/src/equil/vcs_VolPhase.h @@ -431,6 +431,27 @@ namespace VCSnonideal { */ void setExistence(const int existence); + //! Return the Global VCS index of the kth species in the phase + /*! + * @param spIndex local species index (0 to the number of species + * in the phase) + * + * @return Returns the VCS_SOLVE species index of the that species + * This changes as rearrangements are carried out. + */ + int spGlobalIndexVCS(const int spIndex) const; + + + //! set the Global VCS index of the kth species in the phase + /*! + * @param spIndex local species index (0 to the number of species + * in the phase) + * + * @return Returns the VCS_SOLVE species index of the that species + * This changes as rearrangements are carried out. + */ + void setSpGlobalIndexVCS(const int spIndex, const int spGlobalIndex); + private: //! Evaluate the activity coefficients at the current conditions @@ -663,6 +684,7 @@ namespace VCSnonideal { */ int Activity_Coeff_Model; + private: //! Index into the species vectors /*! * Maps the phase species number into the global species number. @@ -671,7 +693,6 @@ namespace VCSnonideal { */ std::vector IndSpecies; - private: //! Vector of Species structures for the species belonging to this phase /*! * The index into this vector is the species index within the phase. diff --git a/Cantera/src/equil/vcs_prob.cpp b/Cantera/src/equil/vcs_prob.cpp index 611ddb261..bb794a061 100644 --- a/Cantera/src/equil/vcs_prob.cpp +++ b/Cantera/src/equil/vcs_prob.cpp @@ -322,7 +322,7 @@ namespace VCSnonideal { Vphase = VPhaseList[iphase]; Vphase->setState_TP(T, PresPA); for (int kindex = 0; kindex < Vphase->NVolSpecies; kindex++) { - int kglob = Vphase->IndSpecies[kindex]; + int kglob = Vphase->spGlobalIndexVCS(kindex); plogf("%16s ", SpName[kglob].c_str()); if (kindex == 0) { plogf("%16s", Vphase->PhaseName.c_str()); @@ -454,7 +454,7 @@ namespace VCSnonideal { * Tell the phase object about the current position of the * species within the global species vector */ - volPhase->IndSpecies[k] = kT; + volPhase->setSpGlobalIndexVCS(k, kT); return kT; } diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index a96a688cc..df1518362 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -725,7 +725,7 @@ namespace VCSnonideal { Vphase = m_VolPhaseList[iph]; for (int k = 0; k < Vphase->NVolSpecies; k++) { vcs_SpeciesProperties *sProp = Vphase->speciesProperty(k); - int kT = Vphase->IndSpecies[k]; + int kT = Vphase->spGlobalIndexVCS(k); sProp->SpeciesThermo = m_speciesThermoList[kT]; } } @@ -746,10 +746,10 @@ namespace VCSnonideal { * So SpecLnMnaught[iSolvent] = 0.0, and the * loop below starts at 1, not 0. */ - int iSolvent = Vphase->IndSpecies[0]; + int iSolvent = Vphase->spGlobalIndexVCS(0); double mnaught = m_wtSpecies[iSolvent] / 1000.; for (int k = 1; k < Vphase->NVolSpecies; k++) { - int kspec = Vphase->IndSpecies[k]; + int kspec = Vphase->spGlobalIndexVCS(k); m_actConventionSpecies[kspec] = Vphase->m_activityConvention; m_lnMnaughtSpecies[kspec] = log(mnaught); } @@ -959,24 +959,16 @@ namespace VCSnonideal { for (int iph = 0; iph < pub->NPhase; iph++) { vcs_VolPhase *pubPhase = pub->VPhaseList[iph]; vcs_VolPhase *vPhase = m_VolPhaseList[iph]; - //pubPhase->setExistence(vPhase->exists()); - // Note pubPhase is not the same as vPhase, since they contain - // different indexing into the solution vector. - // pubPhase->TMoles = vPhase->TMoles; pubPhase->TMolesInert = vPhase->TMolesInert; pubPhase->setTotalMoles(vPhase->TotalMoles()); pubPhase->setElectricPotential(vPhase->electricPotential()); double sumMoles = pubPhase->TMolesInert; pubPhase->setMoleFractions(VCS_DATA_PTR(vPhase->moleFractions())); for (int k = 0; k < pubPhase->NVolSpecies; k++) { - kT = pubPhase->IndSpecies[k]; - //pubPhase->SS0ChemicalPotential[k] = vPhase->SS0ChemicalPotential[k]; - //pubPhase->StarMolarVol[k] = vPhase->StarMolarVol[k]; - //pubPhase->PartialMolarVol[k] = vPhase->PartialMolarVol[k]; - //pubPhase->ActCoeff[k] = vPhase->ActCoeff[k]; + kT = pubPhase->spGlobalIndexVCS(k); if (pubPhase->phiVarIndex() == k) { - k1 = vPhase->IndSpecies[k]; + k1 = vPhase->spGlobalIndexVCS(k); double tmp = m_molNumSpecies_old[k1]; if (! vcs_doubleEqual( pubPhase->electricPotential() , tmp)) { plogf("We have an inconsistency in voltage, %g, %g\n", diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index a72e5eeb3..3e619f18e 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -2658,7 +2658,7 @@ namespace VCSnonideal { } int irxn, kspec; if (Vphase->SingleSpecies) { - kspec = Vphase->IndSpecies[0]; + kspec = Vphase->spGlobalIndexVCS(0); irxn = kspec + m_numComponents; if(m_deltaGRxn_old[irxn] < 0.0) { return true; @@ -2668,7 +2668,7 @@ namespace VCSnonideal { double phaseDG = 1.0; for (int kk = 0; kk < Vphase->NVolSpecies; kk++) { - kspec = Vphase->IndSpecies[kk]; + kspec = Vphase->spGlobalIndexVCS(kk); irxn = kspec + m_numComponents; if (m_deltaGRxn_old[irxn] > 50.0) m_deltaGRxn_old[irxn] = 50.0; if (m_deltaGRxn_old[irxn] < -50.0) m_deltaGRxn_old[irxn] = -50.0; @@ -3457,7 +3457,7 @@ namespace VCSnonideal { if (! Vphase->SingleSpecies) { double sum = 0.0; for (k = 0; k < Vphase->NVolSpecies; k++) { - kspec = Vphase->IndSpecies[k]; + kspec = Vphase->spGlobalIndexVCS(k); if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { sum += molNumSpecies[kspec]; } @@ -3471,7 +3471,7 @@ namespace VCSnonideal { if (lneed) { double poly = 0.0; for (k = 0; k < Vphase->NVolSpecies; k++) { - kspec = Vphase->IndSpecies[k]; + kspec = Vphase->spGlobalIndexVCS(k); irxn = kspec - m_numComponents; // We may need to look at deltaGRxn for components! if (irxn >= 0) { @@ -3486,7 +3486,7 @@ namespace VCSnonideal { * the phase will come back into existence. */ for (k = 0; k < Vphase->NVolSpecies; k++) { - kspec = Vphase->IndSpecies[k]; + kspec = Vphase->spGlobalIndexVCS(k); irxn = kspec - m_numComponents; if (irxn >= 0) { deltaGRxn[irxn] = 1.0 - poly; @@ -4390,7 +4390,7 @@ namespace VCSnonideal { #endif double tMoles = TPhInertMoles[iph]; for (k = 0; k < nkk; k++) { - kspec = Vphase->IndSpecies[k]; + kspec = Vphase->spGlobalIndexVCS(k); tMoles += molNum[kspec]; } double tlogMoles = 0.0; @@ -4405,7 +4405,7 @@ namespace VCSnonideal { double Faraday_phi = m_Faraday_dim * phi; for (k = 0; k < nkk; k++) { - kspec = Vphase->IndSpecies[k]; + kspec = Vphase->spGlobalIndexVCS(k); if (kspec >= m_numComponents) { int irxn = kspec - m_numComponents; if (!do_deleted && @@ -5009,17 +5009,19 @@ namespace VCSnonideal { kp1 = m_speciesLocalPhaseIndex[k1]; kp2 = m_speciesLocalPhaseIndex[k2]; #ifdef DEBUG_MODE - if (pv1->IndSpecies[kp1] != k1) { + if (pv1->spGlobalIndexVCS(kp1) != k1) { plogf("Indexing error in program\n"); exit(-1); } - if (pv2->IndSpecies[kp2] != k2) { + if (pv2->spGlobalIndexVCS(kp2) != k2) { plogf("Indexing error in program\n"); exit(-1); } #endif - pv1->IndSpecies[kp1] = k2; - pv2->IndSpecies[kp2] = k1; + pv1->setSpGlobalIndexVCS(kp1, k2); + pv2->setSpGlobalIndexVCS(kp2, k1); + //pv1->IndSpecies[kp1] = k2; + //pv2->IndSpecies[kp2] = k1; vcsUtil_stsw(m_speciesName, k1, k2); SWAP(m_molNumSpecies_old[k1], m_molNumSpecies_old[k2], t1); @@ -5156,7 +5158,7 @@ namespace VCSnonideal { * Single species Phase */ if (vPhase->SingleSpecies) { - kspec = vPhase->IndSpecies[0]; + kspec = vPhase->spGlobalIndexVCS(0); #ifdef DEBUG_MODE if (iphase != m_phaseID[kspec]) { plogf("vcs_deltag_Phase index error\n"); diff --git a/Cantera/src/equil/vcs_species_thermo.cpp b/Cantera/src/equil/vcs_species_thermo.cpp index d84fdd979..6b8436821 100644 --- a/Cantera/src/equil/vcs_species_thermo.cpp +++ b/Cantera/src/equil/vcs_species_thermo.cpp @@ -391,7 +391,7 @@ double VCS_SOLVE::vcs_Gxs_phase_calc(vcs_VolPhase *Vphase, double *mf_PO) VCS_SPECIES_THERMO *ts_ptr; if (Vphase->Activity_Coeff_Model != VCS_AC_CONSTANT) { for (kspec = 0; kspec < Vphase->NVolSpecies; kspec++) { - kglob = Vphase->IndSpecies[kspec]; + kglob = Vphase->spGlobalIndexVCS(kspec); ts_ptr = m_speciesThermoList[kglob]; ac = ts_ptr->eval_ac(kspec); Gxs += mf_PO[kspec] * log(ac);