research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Changed the algorithm for picking the component basis.

Browse source 
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.
This commit is contained in:
Harry Moffat 2008-04-14 17:22:02 +00:00
parent 6c7e375020
commit daa479d4f7
7 changed files with 72 additions and 40 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
Cantera/src/equil/vcs_internal.h
View file

@ -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
/*!

18
Cantera/src/equil/vcs_prep.cpp
View file

@ -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 ************* */
/* ***************************************************** */

2
Cantera/src/equil/vcs_report.cpp
View file

@ -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);

1
Cantera/src/equil/vcs_solve.cpp
View file

@ -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);

7
Cantera/src/equil/vcs_solve.h
View file

@ -401,6 +401,13 @@ public:
*/
std::vector<double> 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<double> m_spSize;
//! Standard state chemical potentials for species K at the current
//! temperature and pressure.
/*!

22
Cantera/src/equil/vcs_solve_TP.cpp
View file

@ -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;

58
Cantera/src/equil/vcs_util.cpp
View file

@ -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)
/**************************************************************************