From 1e68f099d551c72df402f43b398e69b47d476aac Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 23 Oct 2008 00:37:00 +0000 Subject: [PATCH] Changed the vector m_rxnStatus[] into the vectur m_speciesStatus[]. This is the same information, but now its over kspec instead of irxn. This clarifies the code, and allows for the specification of the status of component species. For example, component species can now be specified to have VCS_SPECIES_ZEROEDMS. --- Cantera/src/equil/vcs_defs.h | 6 +- Cantera/src/equil/vcs_elem.cpp | 11 ++- Cantera/src/equil/vcs_rxnadj.cpp | 14 ++-- Cantera/src/equil/vcs_solve.cpp | 9 ++- Cantera/src/equil/vcs_solve.h | 42 +++++++++-- Cantera/src/equil/vcs_solve_TP.cpp | 114 +++++++++++++++-------------- 6 files changed, 117 insertions(+), 79 deletions(-) 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