From daa479d4f7ec2f3f8c660a6e2ee720533f862871 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Mon, 14 Apr 2008 17:22:02 +0000 Subject: [PATCH] Changed the algorithm for picking the component basis. It used to be that the components were picked based on the magnitude of their mole numbers. I changed this to the magnitude of the mole numbers multiplied by the # elements in the species. On one problem, a heptane combustion mechanism, this reduced the number of iterations in the VCS algorithm from 1080 to 40. Obviously, the convergence rate within the VCS method depends on the ID of the components. This change should be considered as a trial balloon. It may fall apart on some other problem. I may be forestalling the obvious here. The more complete solution is to formualte the hessian of the rxns involving the major species, and then invert it, using cg probably, to get a better iteration direction. --- Cantera/src/equil/vcs_internal.h | 4 ++- Cantera/src/equil/vcs_prep.cpp | 18 ++++++++++ Cantera/src/equil/vcs_report.cpp | 2 +- Cantera/src/equil/vcs_solve.cpp | 1 + Cantera/src/equil/vcs_solve.h | 7 ++++ Cantera/src/equil/vcs_solve_TP.cpp | 22 ++++++------ Cantera/src/equil/vcs_util.cpp | 58 ++++++++++++++++-------------- 7 files changed, 72 insertions(+), 40 deletions(-) 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) /**************************************************************************