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

Fixed MultiPhaseEquil::getComponents when there are more elements than species

Browse source 
When the atomic composition matrix is rank deficient, rows of zeros
will be formed as the elimination proceeds. These rows are now moved
to the bottom of the matrix as they are encountered, so the final
matrix is upper triangular.
This commit is contained in:
Ray Speth 2012-04-12 21:41:31 +00:00
parent 22806f2896
commit 3331c7f3ad
1 changed files with 38 additions and 14 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

52
src/equil/MultiPhaseEquil.cpp
View file

@ -413,8 +413,6 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
m_order[k] = order[k];
}
}
doublereal tmp;
index_t itmp;
index_t nRows = m_nel;
index_t nColumns = m_nsp;
@ -427,15 +425,47 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
}
}
// Do Gauss elimination
// Do Gaussian elimination
for (m = 0; m < nRows; m++) {
// Check for rows that are zero
bool isZeroRow = true;
for (k = m; k < nColumns; k++) {
if (fabs(m_A(m,k)) > sqrt(TINY)) {
isZeroRow = false;
break;
}
}
if (isZeroRow) {
// Find the last non-zero row
index_t n = nRows - 1;
bool foundSwapCandidate = false;
for (; n > m; n--) {
for (k = m; k < nColumns; k++) {
if (fabs(m_A(n,k)) > sqrt(TINY)) {
foundSwapCandidate = true;
break;
}
}
if (foundSwapCandidate) {
break;
}
}
if (m != n) {
// Swap this row with the last non-zero row
for (k = 0; k < nColumns; k++) {
std::swap(m_A(n,k), m_A(m,k));
}
} else {
// All remaining rows are zero. Elimination is complete.
break;
}
}
// If a pivot is zero, exchange columns. This occurs when
// a species has an elemental composition that is not
// linearly independent of the component species that have
// already been assigned
if (m_A(m,m) == 0.0) {
if (m < nColumns && m_A(m,m) == 0.0) {
// First, we need to find a good candidate for a
// component species to swap in for the one that has
// zero pivot. It must contain element m, be linearly
@ -457,16 +487,11 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
// Now exchange the column with zero pivot with the
// column for this major species
for (size_t n = 0; n < nRows; n++) {
tmp = m_A(n,m);
m_A(n, m) = m_A(n, kmax);
m_A(n, kmax) = tmp;
std::swap(m_A(n, m), m_A(n, kmax));
}
// exchange the species labels on the columns
itmp = m_order[m];
m_order[m] = m_order[kmax];
m_order[kmax] = itmp;
std::swap(m_order[m], m_order[kmax]);
}
// scale row m so that the diagonal element is unity
@ -485,7 +510,6 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
}
}
// The left m_nel columns of A are now upper-diagonal. Now
// reduce the m_nel columns to diagonal form by back-solving
for (m = nRows-1; m > 0; m--) {
@ -499,7 +523,7 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
}
}
// create stoichometric coefficient matrix.
// create stoichiometric coefficient matrix.
for (size_t n = 0; n < m_nsp; n++) {
if (n < m_nel)
for (k = 0; k < nFree(); k++) {