diff --git a/Cantera/src/equil/vcs_internal.h b/Cantera/src/equil/vcs_internal.h index 07d8e7fe2..bb7743873 100644 --- a/Cantera/src/equil/vcs_internal.h +++ b/Cantera/src/equil/vcs_internal.h @@ -429,12 +429,14 @@ namespace VCSnonideal { //! Finds the location of the maximum component in a double vector /*! * @param x pointer to a vector of doubles + * @param xSize pointer to a vector of doubles used as a multiplier + * to x[] * @param j lowest index to search from * @param n highest index to search from * @return Return index of the greatest value on X(i) searched * j <= i < n */ - int vcs_amax(const double *x, int j, int n); + int vcs_optMax(const double *x, const double *xSize, int j, int n); //! Returns the maximum integer in a list /*! diff --git a/Cantera/src/equil/vcs_prep.cpp b/Cantera/src/equil/vcs_prep.cpp index 3031a46ea..efe5b3c88 100644 --- a/Cantera/src/equil/vcs_prep.cpp +++ b/Cantera/src/equil/vcs_prep.cpp @@ -20,6 +20,7 @@ #include "vcs_internal.h" #include "vcs_prob.h" #include "vcs_VolPhase.h" +#include "vcs_SpeciesProperties.h" namespace VCSnonideal { @@ -144,6 +145,23 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) ir[i] = m_numElemConstraints + i; } + for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) { + int pID = PhaseID[kspec]; + int spPhIndex = indPhSp[kspec]; + vcs_VolPhase *vPhase = VPhaseList[pID]; + vcs_SpeciesProperties *spProp = vPhase->ListSpeciesPtr[spPhIndex]; + double sz = 0.0; + int eSize = spProp->FormulaMatrixCol.size(); + for (int e = 0; e < eSize; e++) { + sz += fabs(spProp->FormulaMatrixCol[e]); + } + if (sz > 0.0) { + m_spSize[kspec] = sz; + } else { + m_spSize[kspec] = 1.0; + } + } + /* ***************************************************** */ /* **** DETERMINE THE NUMBER OF COMPONENTS ************* */ /* ***************************************************** */ diff --git a/Cantera/src/equil/vcs_report.cpp b/Cantera/src/equil/vcs_report.cpp index 341e6b14e..82c878152 100644 --- a/Cantera/src/equil/vcs_report.cpp +++ b/Cantera/src/equil/vcs_report.cpp @@ -70,7 +70,7 @@ int VCS_SOLVE::vcs_report(int iconv) * the magnitude of the mole fraction vector. */ for (l = m_numComponents; l < m_numSpeciesRdc; ++l) { - k = vcs_amax(VCS_DATA_PTR(xy), l, m_numSpeciesRdc); + k = vcs_optMax(VCS_DATA_PTR(xy), 0, l, m_numSpeciesRdc); if (k != l) { vcsUtil_dsw(VCS_DATA_PTR(xy), k, l); vcsUtil_isw(VCS_DATA_PTR(sortindex), k, l); diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index c9748b6da..05ff5c747 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -107,6 +107,7 @@ namespace VCSnonideal { sc.resize(nspecies0, nelements, 0.0); scSize.resize(nspecies0, 0.0); + m_spSize.resize(nspecies0, 1.0); m_gibbsSpecies.resize(nspecies0, 0.0); ff.resize(nspecies0, 0.0); diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index de5d887e7..f48d9f4d5 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -401,6 +401,13 @@ public: */ std::vector scSize; + //! total size of the species + /*! + * This is used as a multiplier to the mole number in figuring out which + * species should be components. + */ + std::vector m_spSize; + //! Standard state chemical potentials for species K at the current //! temperature and pressure. /*! diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index fc57e78a2..88718bd77 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -1384,7 +1384,7 @@ namespace VCSnonideal { */ dofast = (m_numComponents != 1); for (i = 1; i < m_numComponents; ++i) { - if (soln[i - 1] < soln[i]) { + if ((soln[i - 1] * m_spSize[i-1]) < (soln[i] * m_spSize[i])) { dofast = FALSE; break; } @@ -1394,7 +1394,7 @@ namespace VCSnonideal { for (i = 0; i < m_numRxnRdc; ++i) { l = ir[i]; for (j = m_numComponents - 1; j >= 0; j--) { - if (soln[l] > soln[j]) { + if ((soln[l] * m_spSize[l]) > (soln[j]* m_spSize[j] * 1.01)) { if (sc[i][j] != 0.0) { #ifdef DEBUG_MODE if (vcs_debug_print_lvl >= 2) { @@ -1436,7 +1436,7 @@ namespace VCSnonideal { for (i = 0; i < m_numRxnRdc; ++i) { l = ir[i]; for (j = 0; j < m_numComponents; ++j) { - if (soln[l] > soln[j]) { + if ((soln[l] * m_spSize[l]) > (soln[j] * m_spSize[j] * 1.01)) { if (sc[i][j] != 0.0) { #ifdef DEBUG_MODE if (vcs_debug_print_lvl >= 2) { @@ -3333,7 +3333,7 @@ namespace VCSnonideal { plogf(" --- Species | "); for (j = 0; j < m_numElemConstraints; j++) { plogf(" "); - vcs_print_stringTrunc(ElName[j].c_str(), 4, 1); + vcs_print_stringTrunc(ElName[j].c_str(), 8, 1); } plogf("\n"); for (k = 0; k < m_numSpeciesTot; k++) { @@ -3341,11 +3341,11 @@ namespace VCSnonideal { vcs_print_stringTrunc(SpName[k].c_str(), 11, 1); plogf(" | "); for (j = 0; j < m_numElemConstraints; j++) { - plogf("%5.1g", FormulaMatrix[j][k]); + plogf(" %8.2g", FormulaMatrix[j][k]); } plogf("\n"); } - plogf("\n"); + plogendl(); } } #endif @@ -3388,7 +3388,7 @@ namespace VCSnonideal { * The first search criteria is always the largest positive * magnitude of the mole number. */ - k = vcs_amax(aw, jr, m_numSpeciesTot); + k = vcs_optMax(aw, VCS_DATA_PTR(m_spSize), jr, m_numSpeciesTot); /* * The fun really starts when you have run out of species that have a significant * concentration. It becomes extremely important to make a good choice of which @@ -4593,6 +4593,7 @@ namespace VCSnonideal { SWAP(SpeciesUnknownType[k1], SpeciesUnknownType[k2], j); SWAP(wt[k1], wt[k2], t1); SWAP(ff[k1], ff[k2], t1); + SWAP(m_spSize[k1], m_spSize[k2], t1); SWAP(m_gibbsSpecies[k1], m_gibbsSpecies[k2], t1); SWAP(ds[k1], ds[k2], t1); SWAP(fel[k1], fel[k2], t1); @@ -4624,7 +4625,7 @@ namespace VCSnonideal { if (ifunc) { /* - * Find the noncomponent indecises for the two species + * Find the Rxn indecises corresponding to the two species */ i1 = k1 - m_numComponents; i2 = k2 - m_numComponents; @@ -4646,6 +4647,7 @@ namespace VCSnonideal { } SWAP(dg[i1], dg[i2], t1); SWAP(dgl[i1], dgl[i2], t1); + SWAP(m_deltaGRxn_tmp[i1], m_deltaGRxn_tmp[i2], t1); SWAP(spStatus[i1], spStatus[i2], j); /* @@ -4656,9 +4658,7 @@ namespace VCSnonideal { */ } } /* vcs_switch_pos() ********************************************************/ - /*****************************************************************************/ - /*****************************************************************************/ - /*****************************************************************************/ + static void print_space(int num) { int j; diff --git a/Cantera/src/equil/vcs_util.cpp b/Cantera/src/equil/vcs_util.cpp index 517ce1487..ffbb3c342 100644 --- a/Cantera/src/equil/vcs_util.cpp +++ b/Cantera/src/equil/vcs_util.cpp @@ -162,40 +162,44 @@ namespace VCSnonideal { std::copy(vec_from.begin(), vec_from.begin() + length, vec_to.begin()); } #endif - /*****************************************************************************/ - /*****************************************************************************/ - /*****************************************************************************/ - int vcs_amax(const double *x, int j, int n) - - /************************************************************************** - * - * vcs_amax: - * - * Finds the location of the maximum component in a double vector - * INPUT - * x(*) - Vector to search - * j <= i < n : i is the range of indecises to search in X(*) - * - * RETURN - * return index of the greatest value on X(*) searched - ***************************************************************************/ - { + /* + * + * Finds the location of the maximum component in a double vector + * INPUT + * x(*) - Vector to search + * xSize(*) if nonnull, this is the multiplier vector to be + * multiplied into x(*) before making the decision. + * j <= i < n : i is the range of indecises to search in X(*) + * + * RETURN + * return index of the greatest value on X(*) searched + */ + int vcs_optMax(const double *x, const double * xSize, int j, int n) { int i; int largest = j; double big = x[j]; - for (i = j + 1; i < n; ++i) { - if (x[i] > big) { - largest = i; - big = x[i]; + if (xSize) { + assert(xSize[j] > 0.0); + big *= xSize[j]; + for (i = j + 1; i < n; ++i) { + assert(xSize[i] > 0.0); + if ((x[i]*xSize[i]) > big) { + largest = i; + big = x[i]*xSize[i]; + } + } + } else { + for (i = j + 1; i < n; ++i) { + if (x[i] > big) { + largest = i; + big = x[i]; + } } } return largest; - } /* vcs_amax() **************************************************************/ - /*****************************************************************************/ - /*****************************************************************************/ - /*****************************************************************************/ - + } + int vcs_max_int(const int *vector, int length) /**************************************************************************