diff --git a/Cantera/src/equil/vcs_defs.h b/Cantera/src/equil/vcs_defs.h index 2f466b692..7627fb398 100644 --- a/Cantera/src/equil/vcs_defs.h +++ b/Cantera/src/equil/vcs_defs.h @@ -151,7 +151,7 @@ namespace VCSnonideal { */ #define VCS_SPECIES_SMALLMS -1 - //! Species lies in a multicomponent phase, with concentration zero + //! Species lies in a multicomponent phase with concentration zero /*! * The species lies in a multicomponent phase which currently doesn't exist. * It concentration is currently zero. @@ -191,10 +191,10 @@ namespace VCSnonideal { */ #define VCS_SPECIES_ZEROEDPHASE -6 - //! Species lies in a multicomponent phase that is active, but its concentration is zero + //! Species lies in a multicomponent phase that is active, but species concentration is zero /*! * The species lies in a multicomponent phase which currently does exist. - * It concentration is currently zero, though the phase exists. Note, this + * It concentration is currently identically zero, though the phase exists. Note, this * is a temporary condition that exists at the start of an equilibrium problem. * The species is soon "birthed" or "deleted". */ diff --git a/Cantera/src/equil/vcs_elem.cpp b/Cantera/src/equil/vcs_elem.cpp index 5e2118f8a..19bef9c9e 100644 --- a/Cantera/src/equil/vcs_elem.cpp +++ b/Cantera/src/equil/vcs_elem.cpp @@ -307,7 +307,7 @@ namespace VCSnonideal { int elType = m_elType[i]; if (elType == VCS_ELEM_TYPE_ABSPOS) { for (kspec = 0; kspec < m_numSpeciesTot; kspec++) { - int irxn = kspec - m_numComponents; + // int irxn = kspec - m_numComponents; if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { double atomComp = m_formulaMatrix[i][kspec]; if (atomComp > 0.0) { @@ -327,15 +327,18 @@ namespace VCSnonideal { if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF) { m_molNumSpecies_old[kspec] = 0.0; if (m_SSPhase[kspec]) { - m_rxnStatus[kspec] = VCS_SPECIES_ZEROEDSS; + // m_rxnStatus[kspec] = VCS_SPECIES_ZEROEDSS; + m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS; } else { - m_rxnStatus[kspec] = VCS_SPECIES_ZEROEDMS; + // m_rxnStatus[kspec] = VCS_SPECIES_ZEROEDMS; + m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDMS; } #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { plogf(" --- vcs_elcorr: Zeroed species %s and changed " "status to %d due to max bounds constraint\n", - m_speciesName[kspec].c_str(), m_rxnStatus[irxn]); + // m_speciesName[kspec].c_str(), m_rxnStatus[irxn]); + m_speciesName[kspec].c_str(), m_speciesStatus[kspec]); } #endif } diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index 7772d82dd..05fbd1ec3 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -75,7 +75,7 @@ namespace VCSnonideal { kspec = m_indexRxnToSpecies[irxn]; - if (m_rxnStatus[irxn] == VCS_SPECIES_ZEROEDPHASE) { + if (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) { m_deltaMolNumSpecies[kspec] = 0.0; #ifdef DEBUG_MODE sprintf(ANOTE, "ZeroedPhase: Phase is artificially zeroed"); @@ -152,10 +152,10 @@ namespace VCSnonideal { * Don't calculate for minor or nonexistent species if * their values are to be decreasing anyway. */ - if ((m_rxnStatus[irxn] != VCS_SPECIES_MAJOR) && (m_deltaGRxn_new[irxn] >= 0.0)) { + if ((m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) && (m_deltaGRxn_new[irxn] >= 0.0)) { #ifdef DEBUG_MODE sprintf(ANOTE,"Skipped: IC = %3d and DG >0: %11.3E", - m_rxnStatus[irxn], m_deltaGRxn_new[irxn]); + m_speciesStatus[kspec], m_deltaGRxn_new[irxn]); if (m_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); plogf(" %12.4E %12.4E %12.4E | %s\n", @@ -301,7 +301,7 @@ namespace VCSnonideal { Vphase = m_VolPhaseList[iph]; Vphase->setTotalMoles(0.0); if (k == kspec) { - m_rxnStatus[irxn] = VCS_SPECIES_ZEROEDSS; + m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS; if (m_SSPhase[kspec] != 1) { printf("vcs_RxnStepSizes:: we shouldn't be here!\n"); std::exit(-1); @@ -418,7 +418,7 @@ namespace VCSnonideal { (void) sprintf(ANOTE, "MultSpec: come alive DG = %11.3E", m_deltaGRxn_new[irxn]); #endif m_deltaMolNumSpecies[kspec] = 1.0e-10; - m_rxnStatus[irxn] = VCS_SPECIES_MAJOR; + m_speciesStatus[kspec] = VCS_SPECIES_MAJOR; --(m_numRxnMinorZeroed); } else { #ifdef DEBUG_MODE @@ -450,10 +450,10 @@ namespace VCSnonideal { * Don't calculate for minor or nonexistent species if * their values are to be decreasing anyway. */ - if (m_rxnStatus[irxn] <= VCS_SPECIES_MINOR && m_deltaGRxn_new[irxn] >= 0.0) { + if (m_speciesStatus[kspec] <= VCS_SPECIES_MINOR && m_deltaGRxn_new[irxn] >= 0.0) { #ifdef DEBUG_MODE sprintf(ANOTE,"Skipped: IC = %3d and DG >0: %11.3E\n", - m_rxnStatus[irxn], m_deltaGRxn_new[irxn]); + m_speciesStatus[kspec], m_deltaGRxn_new[irxn]); plogf(" --- "); plogf("%-12.12s", m_speciesName[kspec].c_str()); plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], ANOTE); diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index d9950ee48..f8983be00 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -179,7 +179,8 @@ namespace VCSnonideal { m_indexRxnToSpecies.resize(nspecies0, 0); /* Initialize all species to be major species */ - m_rxnStatus.resize(nspecies0, 1); + //m_rxnStatus.resize(nspecies0, 1); + m_speciesStatus.resize(nspecies0, 1); m_SSPhase.resize(2*nspecies0, 0); m_phaseID.resize(nspecies0, 0); @@ -644,8 +645,10 @@ namespace VCSnonideal { /* * Define all species to be major species, initially. */ - for (i = 0; i < nspecies; i++) m_rxnStatus[i] = VCS_SPECIES_MAJOR; - + for (i = 0; i < nspecies; i++) { + // m_rxnStatus[i] = VCS_SPECIES_MAJOR; + m_speciesStatus[i] = VCS_SPECIES_MAJOR; + } /* * PhaseID: Fill in the species to phase mapping * -> Check for bad values at the same time. diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index c20e842c0..c31ea4d0f 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -424,6 +424,28 @@ public: * * m_feSpecies(I)(I) = m_SSfeSpecies(I) + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF)) * + * Species in the following categories are treated as "small species" + * + * - VCS_SPECIES_DELETED + * - VCS_SPECIES_ACTIVEBUTZERO + * . + * + * Handling of Small Species: + * ------------------------------ + * For species in multispecies phases which are currently not active, the + * treatment is different. These species are in the following species categories: + * + * - VCS_SPECIES_ZEROEDMS + * - VCS_SPECIES_ZEROEDPHASE + * . + * + * For these species, the ln( ActCoeff[I] X[I]) term is + * dropped altogether. The following equation is used. + * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + Charge[I] * Faraday_dim * phasePhi[iphase]; + * + * * Handling of "Species" Representing Interfacial Voltages * --------------------------------------------------------- * @@ -645,8 +667,14 @@ public: int vcs_TP(int ipr, int ip1, int maxit, double T, double pres); int vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres); - void vcs_fePrep_TP(void); + //! Initialize the chemical potential of single species phases + /*! + * For single species phases, initialize the chemical + * potential with the value of the standard state chemical + * potential. This value doesn't change during the calculation + */ + void vcs_fePrep_TP(); //! Calculation of the total volume and the partial molar volumes /*! @@ -1695,14 +1723,14 @@ public: */ std::vector m_indexRxnToSpecies; - //! Major -Minor status vector for the formation reaction + //! Major -Minor status vector for the species in the problem /*! - * The index for this is rxn. The species that this is refereing + * The index for this is species. The reaction that this is referring * to is * kspec = irxn + m_numComponents * - * formation rxn irxn : 1 -> Major player VCS_SPECIES_MAJOR - * 0 -> Minor player VCS_SPECIES_MINOR + * kspec : 1 -> Major species VCS_SPECIES_MAJOR + * 0 -> Minor species VCS_SPECIES_MINOR * -1 -> Mole number is zero * in inactive phase VCS_SPECIES_ZEROEDPHASE * -2 -> Deleted species in an @@ -1711,9 +1739,9 @@ public: * in a stoich phase - VCS_SPECIES_ZEREODSS * -4 -> Species is deleted * - VCS_SPECIES_DELETED - * -> Length equal to number of non-components + * -> Length equal to number of species */ - std::vector m_rxnStatus; + std::vector m_speciesStatus; //! Mapping from the species number to the phase number std::vector m_phaseID; diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 8d90f6389..d1d75a7e7 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -303,11 +303,11 @@ namespace VCSnonideal { #endif for (irxn = 0; irxn < m_numRxnRdc; ++irxn) { kspec = m_indexRxnToSpecies[irxn]; - m_rxnStatus[irxn] = vcs_species_type(kspec); + m_speciesStatus[kspec] = vcs_species_type(kspec); #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { - if (m_rxnStatus[irxn] != VCS_SPECIES_MINOR) { - switch (m_rxnStatus[irxn]) { + if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { + switch (m_speciesStatus[kspec]) { case VCS_SPECIES_MAJOR: plogf(" --- Major Species : %-s\n", m_speciesName[kspec].c_str()); break; @@ -331,14 +331,14 @@ namespace VCSnonideal { plogf(" --- InterfaceVoltage Species: %-s\n", m_speciesName[kspec].c_str()); break; default: - plogf(" --- Unknown type - ERROR %d\n", m_rxnStatus[irxn]); + plogf(" --- Unknown type - ERROR %d\n", m_speciesStatus[kspec]); plogendl(); std::exit(-1); } } } #endif - if (m_rxnStatus[irxn] != VCS_SPECIES_MAJOR) { + if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) { ++m_numRxnMinorZeroed; } } @@ -544,7 +544,7 @@ namespace VCSnonideal { ANOTE[0] = '\0'; #endif - if (m_rxnStatus[irxn] == VCS_SPECIES_INTERFACIALVOLTAGE) { + if (m_speciesStatus[kspec] == VCS_SPECIES_INTERFACIALVOLTAGE) { /********************************************************************/ /************************ VOLTAGE SPECIES ***************************/ /********************************************************************/ @@ -555,7 +555,7 @@ namespace VCSnonideal { #endif m_deltaMolNumSpecies[kspec] = dx; } - else if (m_rxnStatus[irxn] < VCS_SPECIES_MINOR) { + else if (m_speciesStatus[kspec] < VCS_SPECIES_MINOR) { /********************************************************************/ /********************** ZEROED OUT SPECIES **************************/ /********************************************************************/ @@ -563,7 +563,7 @@ namespace VCSnonideal { #ifdef DEBUG_MODE if (m_debug_print_lvl >= 3) { plogf(" --- %s currently zeroed (SpStatus=%-2d):", - m_speciesName[kspec].c_str(), m_rxnStatus[irxn]); + m_speciesName[kspec].c_str(), m_speciesStatus[kspec]); plogf("%3d DG = %11.4E WT = %11.4E W = %11.4E DS = %11.4E\n", irxn, m_deltaGRxn_new[irxn], m_molNumSpecies_new[kspec], m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec]); @@ -626,7 +626,7 @@ namespace VCSnonideal { plogf("%-12s\n", m_speciesName[kspec].c_str()); } #endif - m_rxnStatus[irxn] = VCS_SPECIES_MAJOR; + m_speciesStatus[kspec] = VCS_SPECIES_MAJOR; MajorSpeciesHaveConverged = false; allMinorZeroedSpecies = false; } else { @@ -636,7 +636,7 @@ namespace VCSnonideal { plogf("%-12s\n", m_speciesName[kspec].c_str()); } #endif - m_rxnStatus[irxn] = VCS_SPECIES_MINOR; + m_speciesStatus[kspec] = VCS_SPECIES_MINOR; } if (m_deltaMolNumSpecies[kspec] > 0.0) { dx = m_deltaMolNumSpecies[kspec] * 0.01; @@ -654,7 +654,7 @@ namespace VCSnonideal { m_deltaMolNumSpecies[kspec] = 0.0; dx = 0.0; } - } else if (m_rxnStatus[irxn] == VCS_SPECIES_MINOR) { + } else if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) { /********************************************************************/ /***************************** MINOR SPECIES ************************/ /********************************************************************/ @@ -671,7 +671,7 @@ namespace VCSnonideal { if (m_debug_print_lvl >= 2) { plogf(" --- "); plogf("%-12s", m_speciesName[kspec].c_str()); plogf("%3d%11.4E%11.4E%11.4E | %s", - m_rxnStatus[irxn], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], + m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], m_deltaMolNumSpecies[kspec], ANOTE); plogendl(); } @@ -757,7 +757,7 @@ namespace VCSnonideal { if (m_debug_print_lvl >= 2) { plogf(" --- "); plogf("%-12s", m_speciesName[kspec].c_str()); plogf("%3d%11.4E%11.4E%11.4E | %s", - m_rxnStatus[irxn], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], + m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], m_deltaMolNumSpecies[kspec], ANOTE); plogendl(); } @@ -873,14 +873,14 @@ namespace VCSnonideal { #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { - if (m_rxnStatus[irxn] >= 0) { + if (m_speciesStatus[kspec] >= 0) { plogf(" --- SS species changed to zeroedss: "); plogf("%-12s", m_speciesName[kspec].c_str()); plogendl(); } } #endif - m_rxnStatus[irxn] = VCS_SPECIES_ZEROEDSS; + m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS; ++m_numRxnMinorZeroed; allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc); @@ -974,7 +974,7 @@ namespace VCSnonideal { m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec]; plogf(" --- "); plogf("%-12.12s", m_speciesName[kspec].c_str()); plogf("%3d%11.4E%11.4E%11.4E | %s", - m_rxnStatus[irxn], m_molNumSpecies_old[kspec], + m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], m_deltaMolNumSpecies[kspec], ANOTE); plogendl(); @@ -1140,7 +1140,7 @@ namespace VCSnonideal { for (kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) { irxn = kspec - m_numComponents; plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); - plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", m_rxnStatus[irxn], + plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_old[kspec]+m_deltaMolNumSpecies[kspec], m_molNumSpecies_new[kspec], m_deltaGRxn_old[irxn], @@ -1190,7 +1190,7 @@ namespace VCSnonideal { l1 = i - m_numComponents; plogf(" --- %-12.12s", m_speciesName[i].c_str()); plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", - m_rxnStatus[l1], m_molNumSpecies_old[i], + m_speciesStatus[i], m_molNumSpecies_old[i], m_molNumSpecies_new[i], m_feSpecies_old[i], m_feSpecies_new[i], m_deltaGRxn_old[l1], m_deltaGRxn_new[l1]); } @@ -1198,7 +1198,7 @@ namespace VCSnonideal { l1 = kspec - m_numComponents; plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", - m_rxnStatus[l1], m_molNumSpecies_old[kspec], + m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], m_feSpecies_old[kspec], m_feSpecies_new[kspec], m_deltaGRxn_old[l1], m_deltaGRxn_new[l1]); } @@ -1444,7 +1444,7 @@ namespace VCSnonideal { break; } #ifdef DEBUG_NOT - if (m_rxnStatus[i] == VCS_SPECIES_ZEROEDMS) { + if (m_speciesStatus[l] == VCS_SPECIES_ZEROEDMS) { if (m_molNumSpecies_old[j] == 0.0) { if (m_stoichCoeffRxnMatrix[i][j] != 0.0) { if (dg[i] < 0.0) { @@ -1509,7 +1509,7 @@ namespace VCSnonideal { } } #ifdef DEBUG_NOT - if (m_rxnStatus[i] == VCS_SPECIES_ZEROEDMS) { + if (m_speciesStatus[l] == VCS_SPECIES_ZEROEDMS) { if (m_molNumSpecies_old[j] == 0.0) { if (m_stoichCoeffRxnMatrix[i][j] != 0.0) { if (dg[i] < 0.0) { @@ -1560,7 +1560,7 @@ namespace VCSnonideal { if (speciesType < VCS_SPECIES_MINOR) { #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { - if (m_rxnStatus[irxn] >= VCS_SPECIES_MINOR) { + if (m_speciesStatus[kspec] >= VCS_SPECIES_MINOR) { plogf(" --- major/minor species is now zeroed out: %s\n", m_speciesName[kspec].c_str()); } @@ -1570,10 +1570,10 @@ namespace VCSnonideal { } else if (speciesType == VCS_SPECIES_MINOR) { #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { - if (m_rxnStatus[irxn] != VCS_SPECIES_MINOR) { - if (m_rxnStatus[irxn] == VCS_SPECIES_MAJOR) { + if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { + if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR) { plogf(" --- Noncomponent turned from major to minor: "); - } else if (m_rxnStatus[irxn] == VCS_SPECIES_COMPONENT) { + } else if (m_speciesStatus[kspec] == VCS_SPECIES_COMPONENT) { plogf(" --- Component turned into a minor species: "); } else { plogf(" --- Zeroed Species turned into a " @@ -1585,12 +1585,12 @@ namespace VCSnonideal { #endif ++m_numRxnMinorZeroed; } else if (speciesType == VCS_SPECIES_MAJOR) { - if (m_rxnStatus[irxn] != VCS_SPECIES_MAJOR) { + if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) { #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { - if (m_rxnStatus[irxn] == VCS_SPECIES_MINOR) { + if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) { plogf(" --- Noncomponent turned from minor to major: "); - } else if (m_rxnStatus[irxn] == VCS_SPECIES_COMPONENT) { + } else if (m_speciesStatus[kspec] == VCS_SPECIES_COMPONENT) { plogf(" --- Component turned into a major: "); } else { plogf(" --- Noncomponent turned from zeroed to major: "); @@ -1598,7 +1598,7 @@ namespace VCSnonideal { plogf("%s\n", m_speciesName[kspec].c_str()); } #endif - m_rxnStatus[irxn] = VCS_SPECIES_MAJOR; + m_speciesStatus[kspec] = VCS_SPECIES_MAJOR; /* * For this special case, we must reevaluate thermo functions */ @@ -1609,7 +1609,7 @@ namespace VCSnonideal { } } } - m_rxnStatus[irxn] = speciesType; + m_speciesStatus[kspec] = speciesType; } /* * This logical variable indicates whether all current @@ -1628,7 +1628,8 @@ namespace VCSnonideal { } #endif for (irxn = 0; irxn < m_numRxnRdc; ++irxn) { - if (m_rxnStatus[irxn] == VCS_SPECIES_MAJOR && (fabs(m_deltaGRxn_new[irxn]) > m_tolmaj)) { + kspec = irxn + m_numComponents; + if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR && (fabs(m_deltaGRxn_new[irxn]) > m_tolmaj)) { if (m_VCount->Its >= maxit) { solveFail = -1; /* @@ -1696,7 +1697,8 @@ namespace VCSnonideal { } #endif for (irxn = 0; irxn < m_numRxnRdc; ++irxn) { - if (m_rxnStatus[irxn] == VCS_SPECIES_MINOR && (fabs(m_deltaGRxn_new[irxn]) > m_tolmin)) { + kspec = irxn + m_numComponents; + if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR && (fabs(m_deltaGRxn_new[irxn]) > m_tolmin)) { if (m_VCount->Its >= maxit) { solveFail = -1; /* @@ -2282,8 +2284,8 @@ namespace VCSnonideal { * Decrement the minor species counter if the current species is * a minor species */ - if (m_rxnStatus[irxn] != VCS_SPECIES_MAJOR) --(m_numRxnMinorZeroed); - m_rxnStatus[irxn] = VCS_SPECIES_DELETED; + if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) --(m_numRxnMinorZeroed); + m_speciesStatus[kspec] = VCS_SPECIES_DELETED; m_deltaGRxn_new[irxn] = 0.0; m_deltaGRxn_old[irxn] = 0.0; m_feSpecies_new[kspec] = 0.0; @@ -2357,7 +2359,8 @@ namespace VCSnonideal { * This routine is responsible for the global data manipulation only. */ void VCS_SOLVE::vcs_reinsert_deleted(int kspec) { - int i, k, irxn = kspec - m_numComponents; + int i, k; + // int irxn = kspec - m_numComponents; int iph = m_phaseID[kspec]; double dx; #ifdef DEBUG_MODE @@ -2372,10 +2375,10 @@ namespace VCSnonideal { */ dx = m_tPhaseMoles_old[iph] * VCS_RELDELETE_SPECIES_CUTOFF * 10.; delta_species(kspec, &dx); - m_rxnStatus[irxn] = VCS_SPECIES_MINOR; + m_speciesStatus[kspec] = VCS_SPECIES_MINOR; if (m_SSPhase[kspec]) { - m_rxnStatus[irxn] = VCS_SPECIES_MAJOR; + m_speciesStatus[kspec] = VCS_SPECIES_MAJOR; --(m_numRxnMinorZeroed); } @@ -2387,8 +2390,8 @@ namespace VCSnonideal { * We may have popped a multispecies phase back * into existence. If we did, we have to check * the other species in that phase. - * Take care of the m_rxnStatus[] flag. - * The value of m_rxnStatus[] must change from + * Take care of the m_speciesStatus[] flag. + * The value of m_speciesStatus[] must change from * VCS_SPECIES_ZEROEDPHASE to VCS_SPECIES_ZEROEDMS * for those other species. */ @@ -2398,8 +2401,8 @@ namespace VCSnonideal { for (k = 0; k < m_numSpeciesTot; k++) { if (m_phaseID[k] == iph) { i = k - m_numComponents; - if (m_rxnStatus[i] != VCS_SPECIES_DELETED) { - m_rxnStatus[i] = VCS_SPECIES_MINOR; + if (m_speciesStatus[k] != VCS_SPECIES_DELETED) { + m_speciesStatus[k] = VCS_SPECIES_MINOR; } } } @@ -2494,7 +2497,7 @@ namespace VCSnonideal { * Change the status flag of the species to that of an * zeroed phase */ - m_rxnStatus[irxn] = VCS_SPECIES_ZEROEDMS; + m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDMS; } } } @@ -2570,7 +2573,7 @@ namespace VCSnonideal { m_molNumSpecies_old[kspec] = 0.0; m_molNumSpecies_new[kspec] = 0.0; m_deltaMolNumSpecies[kspec] = 0.0; - m_rxnStatus[irxn] = VCS_SPECIES_ZEROEDMS; + m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDMS; ++(m_numRxnRdc); ++(m_numSpeciesRdc); @@ -3957,9 +3960,8 @@ namespace VCSnonideal { for (k = 0; k < nkk; k++) { kspec = Vphase->spGlobalIndexVCS(k); if (kspec >= m_numComponents) { - int irxn = kspec - m_numComponents; if (!do_deleted && - (m_rxnStatus[irxn] == VCS_SPECIES_DELETED)) { + (m_speciesStatus[kspec] == VCS_SPECIES_DELETED)) { continue; } } @@ -4283,8 +4285,8 @@ namespace VCSnonideal { /* ************************************************ */ if (ll < 0) { for (irxn = 0; irxn < m_numRxnRdc; ++irxn) { - if (m_rxnStatus[irxn] != VCS_SPECIES_MINOR) { - kspec = m_indexRxnToSpecies[irxn]; + kspec = m_indexRxnToSpecies[irxn]; + if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { iphase = m_phaseID[kspec]; if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { #ifdef DEBUG_MODE @@ -4330,8 +4332,8 @@ namespace VCSnonideal { /* ************************************************ */ } else if (ll > 0) { for (irxn = 0; irxn < m_numRxnRdc; ++irxn) { - if (m_rxnStatus[irxn] == VCS_SPECIES_MINOR) { - kspec = m_indexRxnToSpecies[irxn]; + kspec = m_indexRxnToSpecies[irxn]; + if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) { iphase = m_phaseID[kspec]; if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { #ifdef DEBUG_MODE @@ -4425,9 +4427,10 @@ namespace VCSnonideal { int irxn; if (m_numRxnRdc <= 0) return 0.0; for (irxn = 0, tmp = 0.0; irxn < m_numRxnRdc; ++irxn) { - if (m_rxnStatus[irxn] == VCS_SPECIES_MAJOR || m_rxnStatus[irxn] == VCS_SPECIES_MINOR || + int kspec = irxn + m_numComponents; + if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR || m_speciesStatus[kspec] == VCS_SPECIES_MINOR || dgLocal[irxn] < 0.0) { - if (m_rxnStatus[irxn] != VCS_SPECIES_ZEROEDMS) { + if (m_speciesStatus[kspec] != VCS_SPECIES_ZEROEDMS) { tmp += dgLocal[irxn] * dgLocal[irxn]; } } @@ -4637,7 +4640,8 @@ namespace VCSnonideal { /* ************************************************* */ if (l < 0) { for (irxn = 0; irxn < m_numRxnRdc; ++irxn) { - if (m_rxnStatus[irxn] != VCS_SPECIES_MINOR) { + kspec = irxn + m_numComponents; + if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) { icase = 0; deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]]; dtmp_ptr = m_stoichCoeffRxnMatrix[irxn]; @@ -4676,7 +4680,8 @@ namespace VCSnonideal { /* **** MINORS AND ZEROED SPECIES ****************** */ /* ************************************************* */ for (irxn = 0; irxn < m_numRxnRdc; ++irxn) { - if (m_rxnStatus[irxn] <= VCS_SPECIES_MINOR) { + kspec = irxn + m_numComponents; + if (m_speciesStatus[kspec] <= VCS_SPECIES_MINOR) { icase = 0; deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]]; dtmp_ptr = m_stoichCoeffRxnMatrix[irxn]; @@ -5033,7 +5038,7 @@ namespace VCSnonideal { if (m_useActCoeffJac) { vcs_switch2D(m_dLnActCoeffdMolNum.baseDataAddr(), k1, k2); } - + SWAP(m_speciesStatus[k1], m_speciesStatus[k2], j); /* * Handle the index pointer in the phase structures */ @@ -5064,7 +5069,6 @@ namespace VCSnonideal { SWAP(m_deltaGRxn_new[i1], m_deltaGRxn_new[i2], t1); SWAP(m_deltaGRxn_old[i1], m_deltaGRxn_old[i2], t1); SWAP(m_deltaGRxn_tmp[i1], m_deltaGRxn_tmp[i2], t1); - SWAP(m_rxnStatus[i1], m_rxnStatus[i2], j); /* * We don't want to swap ir[], because the values of ir should