25ba149aab
Moved the external libraries to separate library files so that libcantera.a just contains its own namespace externals. Fixed several errors in the equilibrium program that occurred during the port. (int to size_t issues). Moved some equilibrium program headers to the include file system, so that it can link with equilibrium program. Worked on Cantera.mak. Needs more work. Fixed an issue with the Residual virtual base classes within numerics. They didn't inherit due to int to size_t migration. This caused numerous test problems to fail (issue with backwards compatibility - do we want it and how much do we want?). Added csvdiff back so that it's available for shell environment runtests.
894 lines
30 KiB
C++
894 lines
30 KiB
C++
/**
|
|
* @file BasisOptimize.cpp
|
|
* Functions which calculation optimized basis of the
|
|
* stoichiometric coefficient matrix (see /ref equil functions)
|
|
*/
|
|
#include "cantera/base/ct_defs.h"
|
|
#include "cantera/base/global.h"
|
|
#include "cantera/thermo/ThermoPhase.h"
|
|
#include "cantera/equil/MultiPhase.h"
|
|
|
|
#include <cstring>
|
|
|
|
using namespace Cantera;
|
|
using namespace std;
|
|
|
|
#ifdef DEBUG_MODE
|
|
namespace Cantera
|
|
{
|
|
int BasisOptimize_print_lvl = 0;
|
|
}
|
|
|
|
//! Print a string within a given space limit. This routine limits the amount of the string that will be printed to a
|
|
//! maximum of "space" characters.
|
|
/*!
|
|
*
|
|
* @param str String -> must be null terminated.
|
|
* @param space space limit for the printing.
|
|
* @param alignment 0 centered
|
|
* 1 right aligned
|
|
* 2 left aligned
|
|
*/
|
|
static void print_stringTrunc(const char* str, int space, int alignment);
|
|
#endif
|
|
|
|
|
|
//! Finds the location of the maximum component in a double vector INPUT
|
|
/*!
|
|
* @param x Vector to search
|
|
* @param j j <= i < n : i is the range of indices to search in X(*)
|
|
* @param n Length of the vector
|
|
*
|
|
* @return index of the greatest value on X(*) searched
|
|
*/
|
|
static size_t amax(double* x, size_t j, size_t n);
|
|
|
|
//! Invert an nxn matrix and solve m rhs's
|
|
/*!
|
|
*
|
|
* Solve C X + B = 0;
|
|
*
|
|
* This routine uses Gauss elimination and is optimized for the solution
|
|
* of lots of rhs's.
|
|
* A crude form of row pivoting is used here.
|
|
*
|
|
* @param c C is the matrix to be inverted
|
|
* @param idem first dimension in the calling routine
|
|
* idem >= n must be true
|
|
* @param n number of rows and columns in the matrix
|
|
* @param b rhs of the matrix problem
|
|
* @param m number of rhs to be solved for
|
|
*
|
|
* c[i+j*idem] = c_i_j = Matrix to be inverted: i = row number
|
|
* j = column number
|
|
* b[i+j*idem] = b_i_j = vectors of rhs's: i = row number
|
|
* j = column number
|
|
* (each column is a new rhs)
|
|
*
|
|
* @return Retuns the value
|
|
* 1 : Matrix is singular
|
|
* 0 : solution is OK
|
|
*
|
|
* The solution is returned in the matrix b.
|
|
*/
|
|
static int mlequ(double* c, size_t idem, size_t n, double* b, size_t m);
|
|
|
|
/*
|
|
* Choose the optimum basis for the calculations. This is done by
|
|
* choosing the species with the largest mole fraction
|
|
* not currently a linear combination of the previous components.
|
|
* Then, calculate the stoichiometric coefficient matrix for that
|
|
* basis.
|
|
*
|
|
* Calculates the identity of the component species in the mechanism.
|
|
* Rearranges the solution data to put the component data at the
|
|
* front of the species list.
|
|
*
|
|
* Then, calculates SC(J,I) the formation reactions for all noncomponent
|
|
* species in the mechanism.
|
|
*
|
|
* Input
|
|
* ---------
|
|
* mphase Pointer to the multiphase object. Contains the
|
|
* species mole fractions, which are used to pick the
|
|
* current optimal species component basis.
|
|
* orderVectorElement
|
|
* Order vector for the elements. The element rows
|
|
* in the formula matrix are
|
|
* rearranged according to this vector.
|
|
* orderVectorSpecies
|
|
* Order vector for the species. The species are
|
|
* rearranged according to this formula. The first
|
|
* nCompoments of this vector contain the calculated
|
|
* species components on exit.
|
|
* doFormRxn If true, the routine calculates the formation
|
|
* reaction matrix based on the calculated
|
|
* component species. If false, this step is skipped.
|
|
*
|
|
* Output
|
|
* ---------
|
|
* usedZeroedSpecies = If true, then a species with a zero concentration
|
|
* was used as a component. The problem may be
|
|
* converged.
|
|
* formRxnMatrix
|
|
*
|
|
* Return
|
|
* --------------
|
|
* returns the number of components.
|
|
*
|
|
*
|
|
*/
|
|
size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn,
|
|
MultiPhase* mphase, std::vector<size_t>& orderVectorSpecies,
|
|
std::vector<size_t>& orderVectorElements,
|
|
vector_fp& formRxnMatrix)
|
|
{
|
|
|
|
size_t j, jj, k=0, kk, l, i, jl, ml;
|
|
bool lindep;
|
|
std::string ename;
|
|
std::string sname;
|
|
/*
|
|
* Get the total number of elements defined in the multiphase object
|
|
*/
|
|
size_t ne = mphase->nElements();
|
|
/*
|
|
* Get the total number of species in the multiphase object
|
|
*/
|
|
size_t nspecies = mphase->nSpecies();
|
|
doublereal tmp;
|
|
doublereal const USEDBEFORE = -1;
|
|
|
|
/*
|
|
* Perhaps, initialize the element ordering
|
|
*/
|
|
if (orderVectorElements.size() < ne) {
|
|
orderVectorElements.resize(ne);
|
|
for (j = 0; j < ne; j++) {
|
|
orderVectorElements[j] = j;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Perhaps, initialize the species ordering
|
|
*/
|
|
if (orderVectorSpecies.size() != nspecies) {
|
|
orderVectorSpecies.resize(nspecies);
|
|
for (k = 0; k < nspecies; k++) {
|
|
orderVectorSpecies[k] = k;
|
|
}
|
|
}
|
|
|
|
#ifdef DEBUG_MODE
|
|
double molSave = 0.0;
|
|
if (BasisOptimize_print_lvl >= 1) {
|
|
writelog(" ");
|
|
for (i=0; i<77; i++) {
|
|
writelog("-");
|
|
}
|
|
writelog("\n");
|
|
writelog(" --- Subroutine BASOPT called to ");
|
|
writelog("calculate the number of components and ");
|
|
writelog("evaluate the formation matrix\n");
|
|
if (BasisOptimize_print_lvl > 0) {
|
|
writelog(" ---\n");
|
|
|
|
writelog(" --- Formula Matrix used in BASOPT calculation\n");
|
|
writelog(" --- Species | Order | ");
|
|
for (j = 0; j < ne; j++) {
|
|
jj = orderVectorElements[j];
|
|
writelog(" ");
|
|
ename = mphase->elementName(jj);
|
|
print_stringTrunc(ename.c_str(), 4, 1);
|
|
writelogf("(%1d)", j);
|
|
}
|
|
writelog("\n");
|
|
for (k = 0; k < nspecies; k++) {
|
|
kk = orderVectorSpecies[k];
|
|
writelog(" --- ");
|
|
sname = mphase->speciesName(kk);
|
|
print_stringTrunc(sname.c_str(), 11, 1);
|
|
writelogf(" | %4d |", k);
|
|
for (j = 0; j < ne; j++) {
|
|
jj = orderVectorElements[j];
|
|
double num = mphase->nAtoms(kk,jj);
|
|
writelogf("%6.1g ", num);
|
|
}
|
|
writelog("\n");
|
|
}
|
|
writelog(" --- \n");
|
|
}
|
|
}
|
|
#endif
|
|
|
|
/*
|
|
* Calculate the maximum value of the number of components possible
|
|
* It's equal to the minimum of the number of elements and the
|
|
* number of total species.
|
|
*/
|
|
size_t nComponents = std::min(ne, nspecies);
|
|
size_t nNonComponents = nspecies - nComponents;
|
|
/*
|
|
* Set this return variable to false
|
|
*/
|
|
*usedZeroedSpecies = false;
|
|
|
|
/*
|
|
* Create an array of mole numbers
|
|
*/
|
|
vector_fp molNum(nspecies,0.0);
|
|
mphase->getMoles(DATA_PTR(molNum));
|
|
|
|
/*
|
|
* Other workspace
|
|
*/
|
|
vector_fp sm(ne*ne, 0.0);
|
|
vector_fp ss(ne, 0.0);
|
|
vector_fp sa(ne, 0.0);
|
|
if (formRxnMatrix.size() < nspecies*ne) {
|
|
formRxnMatrix.resize(nspecies*ne, 0.0);
|
|
}
|
|
|
|
#ifdef DEBUG_MODE
|
|
/*
|
|
* For debugging purposes keep an unmodified copy of the array.
|
|
*/
|
|
vector_fp molNumBase(molNum);
|
|
#endif
|
|
|
|
|
|
size_t jr = npos;
|
|
/*
|
|
* Top of a loop of some sort based on the index JR. JR is the
|
|
* current number of component species found.
|
|
*/
|
|
do {
|
|
++jr;
|
|
/* - Top of another loop point based on finding a linearly */
|
|
/* - independent species */
|
|
do {
|
|
/*
|
|
* Search the remaining part of the mole number vector, molNum
|
|
* for the largest remaining species. Return its identity.
|
|
* kk is the raw number. k is the orderVectorSpecies index.
|
|
*/
|
|
kk = amax(DATA_PTR(molNum), 0, nspecies);
|
|
for (j = 0; j < nspecies; j++) {
|
|
if (orderVectorSpecies[j] == kk) {
|
|
k = j;
|
|
break;
|
|
}
|
|
}
|
|
if (j == nspecies) {
|
|
throw CanteraError("BasisOptimize", "orderVectorSpecies contains an error");
|
|
}
|
|
|
|
if (molNum[kk] == 0.0) {
|
|
*usedZeroedSpecies = true;
|
|
}
|
|
/*
|
|
* If the largest molNum is negative, then we are done.
|
|
*/
|
|
if (molNum[kk] == USEDBEFORE) {
|
|
nComponents = jr;
|
|
nNonComponents = nspecies - nComponents;
|
|
goto L_END_LOOP;
|
|
}
|
|
/*
|
|
* Assign a small negative number to the component that we have
|
|
* just found, in order to take it out of further consideration.
|
|
*/
|
|
#ifdef DEBUG_MODE
|
|
molSave = molNum[kk];
|
|
#endif
|
|
molNum[kk] = USEDBEFORE;
|
|
|
|
/* *********************************************************** */
|
|
/* **** CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES ****** */
|
|
/* *********************************************************** */
|
|
/*
|
|
* Modified Gram-Schmidt Method, p. 202 Dalquist
|
|
* QR factorization of a matrix without row pivoting.
|
|
*/
|
|
jl = jr;
|
|
for (j = 0; j < ne; ++j) {
|
|
jj = orderVectorElements[j];
|
|
sm[j + jr*ne] = mphase->nAtoms(kk,jj);
|
|
}
|
|
if (jl > 0) {
|
|
/*
|
|
* Compute the coefficients of JA column of the
|
|
* the upper triangular R matrix, SS(J) = R_J_JR
|
|
* (this is slightly different than Dalquist)
|
|
* R_JA_JA = 1
|
|
*/
|
|
for (j = 0; j < jl; ++j) {
|
|
ss[j] = 0.0;
|
|
for (i = 0; i < ne; ++i) {
|
|
ss[j] += sm[i + jr*ne] * sm[i + j*ne];
|
|
}
|
|
ss[j] /= sa[j];
|
|
}
|
|
/*
|
|
* Now make the new column, (*,JR), orthogonal to the
|
|
* previous columns
|
|
*/
|
|
for (j = 0; j < jl; ++j) {
|
|
for (l = 0; l < ne; ++l) {
|
|
sm[l + jr*ne] -= ss[j] * sm[l + j*ne];
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* Find the new length of the new column in Q.
|
|
* It will be used in the denominator in future row calcs.
|
|
*/
|
|
sa[jr] = 0.0;
|
|
for (ml = 0; ml < ne; ++ml) {
|
|
tmp = sm[ml + jr*ne];
|
|
sa[jr] += tmp * tmp;
|
|
}
|
|
/* **************************************************** */
|
|
/* **** IF NORM OF NEW ROW .LT. 1E-3 REJECT ********** */
|
|
/* **************************************************** */
|
|
if (sa[jr] < 1.0e-6) {
|
|
lindep = true;
|
|
} else {
|
|
lindep = false;
|
|
}
|
|
} while (lindep);
|
|
/* ****************************************** */
|
|
/* **** REARRANGE THE DATA ****************** */
|
|
/* ****************************************** */
|
|
if (jr != k) {
|
|
#ifdef DEBUG_MODE
|
|
if (BasisOptimize_print_lvl >= 1) {
|
|
kk = orderVectorSpecies[k];
|
|
sname = mphase->speciesName(kk);
|
|
writelogf(" --- %-12.12s", sname.c_str());
|
|
jj = orderVectorSpecies[jr];
|
|
ename = mphase->speciesName(jj);
|
|
writelogf("(%9.2g) replaces %-12.12s", molSave, ename.c_str());
|
|
writelogf("(%9.2g) as component %3d\n", molNum[jj], jr);
|
|
}
|
|
#endif
|
|
std::swap(orderVectorSpecies[jr], orderVectorSpecies[k]);
|
|
}
|
|
/* - entry point from up above */
|
|
L_END_LOOP:
|
|
;
|
|
/*
|
|
* If we haven't found enough components, go back
|
|
* and find some more. (nc -1 is used below, because
|
|
* jr is counted from 0, via the C convention.
|
|
*/
|
|
} while (jr < (nComponents-1));
|
|
|
|
|
|
if (! doFormRxn) {
|
|
return nComponents;
|
|
}
|
|
|
|
/* ****************************************************** */
|
|
/* **** EVALUATE THE STOICHIOMETRY ********************** */
|
|
/* ****************************************************** */
|
|
/*
|
|
* Formulate the matrix problem for the stoichiometric
|
|
* coefficients. CX + B = 0
|
|
* C will be an nc x nc matrix made up of the formula
|
|
* vectors for the components. Each component's formula
|
|
* vector is a column. The rows are the elements.
|
|
* n rhs's will be solved for. Thus, B is an nc x n
|
|
* matrix.
|
|
*
|
|
* BIG PROBLEM 1/21/99:
|
|
*
|
|
* This algorithm makes the assumption that the
|
|
* first nc rows of the formula matrix aren't rank deficient.
|
|
* However, this might not be the case. For example, assume
|
|
* that the first element in FormulaMatrix[] is argon. Assume that
|
|
* no species in the matrix problem actually includes argon.
|
|
* Then, the first row in sm[], below will be indentically
|
|
* zero. bleh.
|
|
* What needs to be done is to perform a rearrangement
|
|
* of the ELEMENTS -> i.e. rearrange, FormulaMatrix, sp, and gai, such
|
|
* that the first nc elements form in combination with the
|
|
* nc components create an invertible sm[]. not a small
|
|
* project, but very doable.
|
|
* An alternative would be to turn the matrix problem
|
|
* below into an ne x nc problem, and do QR elimination instead
|
|
* of Gauss-Jordon elimination.
|
|
* Note the rearrangement of elements need only be done once
|
|
* in the problem. It's actually very similar to the top of
|
|
* this program with ne being the species and nc being the
|
|
* elements!!
|
|
*/
|
|
for (k = 0; k < nComponents; ++k) {
|
|
kk = orderVectorSpecies[k];
|
|
for (j = 0; j < nComponents; ++j) {
|
|
jj = orderVectorElements[j];
|
|
sm[j + k*ne] = mphase->nAtoms(kk, jj);
|
|
}
|
|
}
|
|
|
|
for (i = 0; i < nNonComponents; ++i) {
|
|
k = nComponents + i;
|
|
kk = orderVectorSpecies[k];
|
|
for (j = 0; j < nComponents; ++j) {
|
|
jj = orderVectorElements[j];
|
|
formRxnMatrix[j + i * ne] = mphase->nAtoms(kk, jj);
|
|
}
|
|
}
|
|
/*
|
|
* Use Gauss-Jordon block elimination to calculate
|
|
* the reaction matrix
|
|
*/
|
|
int ierr = mlequ(DATA_PTR(sm), ne, nComponents, DATA_PTR(formRxnMatrix), nNonComponents);
|
|
if (ierr == 1) {
|
|
writelog("ERROR: mlequ returned an error condition\n");
|
|
throw CanteraError("basopt", "mlequ returned an error condition");
|
|
}
|
|
|
|
#ifdef DEBUG_MODE
|
|
if (Cantera::BasisOptimize_print_lvl >= 1) {
|
|
writelog(" ---\n");
|
|
writelogf(" --- Number of Components = %d\n", nComponents);
|
|
writelog(" --- Formula Matrix:\n");
|
|
writelog(" --- Components: ");
|
|
for (k = 0; k < nComponents; k++) {
|
|
kk = orderVectorSpecies[k];
|
|
writelogf(" %3d (%3d) ", k, kk);
|
|
}
|
|
writelog("\n --- Components Moles: ");
|
|
for (k = 0; k < nComponents; k++) {
|
|
kk = orderVectorSpecies[k];
|
|
writelogf("%-11.3g", molNumBase[kk]);
|
|
}
|
|
writelog("\n --- NonComponent | Moles | ");
|
|
for (i = 0; i < nComponents; i++) {
|
|
kk = orderVectorSpecies[i];
|
|
sname = mphase->speciesName(kk);
|
|
writelogf("%-11.10s", sname.c_str());
|
|
}
|
|
writelog("\n");
|
|
|
|
for (i = 0; i < nNonComponents; i++) {
|
|
k = i + nComponents;
|
|
kk = orderVectorSpecies[k];
|
|
writelogf(" --- %3d (%3d) ", k, kk);
|
|
sname = mphase->speciesName(kk);
|
|
writelogf("%-10.10s", sname.c_str());
|
|
writelogf("|%10.3g|", molNumBase[kk]);
|
|
/*
|
|
* Print the negative of formRxnMatrix[]; it's easier to interpret.
|
|
*/
|
|
for (j = 0; j < nComponents; j++) {
|
|
writelogf(" %6.2f", - formRxnMatrix[j + i * ne]);
|
|
}
|
|
writelog("\n");
|
|
}
|
|
writelog(" ");
|
|
for (i=0; i<77; i++) {
|
|
writelog("-");
|
|
}
|
|
writelog("\n");
|
|
}
|
|
#endif
|
|
|
|
return nComponents;
|
|
} /* basopt() ************************************************************/
|
|
|
|
|
|
|
|
#ifdef DEBUG_MODE
|
|
static void print_stringTrunc(const char* str, int space, int alignment)
|
|
|
|
/***********************************************************************
|
|
* vcs_print_stringTrunc():
|
|
*
|
|
* Print a string within a given space limit. This routine
|
|
* limits the amount of the string that will be printed to a
|
|
* maximum of "space" characters.
|
|
*
|
|
* str = String -> must be null terminated.
|
|
* space = space limit for the printing.
|
|
* alignment = 0 centered
|
|
* 1 right aligned
|
|
* 2 left aligned
|
|
***********************************************************************/
|
|
{
|
|
int i, ls=0, rs=0;
|
|
int len = strlen(str);
|
|
if ((len) >= space) {
|
|
for (i = 0; i < space; i++) {
|
|
writelogf("%c", str[i]);
|
|
}
|
|
} else {
|
|
if (alignment == 1) {
|
|
ls = space - len;
|
|
} else if (alignment == 2) {
|
|
rs = space - len;
|
|
} else {
|
|
ls = (space - len) / 2;
|
|
rs = space - len - ls;
|
|
}
|
|
if (ls != 0) {
|
|
for (i = 0; i < ls; i++) {
|
|
writelog(" ");
|
|
}
|
|
}
|
|
writelogf("%s", str);
|
|
if (rs != 0) {
|
|
for (i = 0; i < rs; i++) {
|
|
writelog(" ");
|
|
}
|
|
}
|
|
}
|
|
}
|
|
#endif
|
|
|
|
/*
|
|
* Finds the location of the maximum component in a double vector
|
|
* INPUT
|
|
* x(*) - Vector to search
|
|
* j <= i < n : i is the range of indices to search in X(*)
|
|
*
|
|
* RETURN
|
|
* return index of the greatest value on X(*) searched
|
|
*/
|
|
static size_t amax(double* x, size_t j, size_t n)
|
|
{
|
|
size_t largest = j;
|
|
double big = x[j];
|
|
for (size_t i = j + 1; i < n; ++i) {
|
|
if (x[i] > big) {
|
|
largest = i;
|
|
big = x[i];
|
|
}
|
|
}
|
|
return largest;
|
|
}
|
|
|
|
/*
|
|
* vcs_mlequ:
|
|
*
|
|
* Invert an nxn matrix and solve m rhs's
|
|
*
|
|
* Solve C X + B = 0;
|
|
*
|
|
* This routine uses Gauss elimination and is optimized for the solution
|
|
* of lots of rhs's.
|
|
* A crude form of row pivoting is used here.
|
|
*
|
|
*
|
|
* c[i+j*idem] = c_i_j = Matrix to be inverted: i = row number
|
|
* j = column number
|
|
* b[i+j*idem] = b_i_j = vectors of rhs's: i = row number
|
|
* j = column number
|
|
* (each column is a new rhs)
|
|
* n = number of rows and columns in the matrix
|
|
* m = number of rhs to be solved for
|
|
* idem = first dimension in the calling routine
|
|
* idem >= n must be true
|
|
*
|
|
* Return Value
|
|
* 1 : Matrix is singular
|
|
* 0 : solution is OK
|
|
*
|
|
* The solution is returned in the matrix b.
|
|
*/
|
|
static int mlequ(double* c, size_t idem, size_t n, double* b, size_t m)
|
|
{
|
|
size_t i, j, k, l;
|
|
double R;
|
|
|
|
/*
|
|
* Loop over the rows
|
|
* -> At the end of each loop, the only nonzero entry in the column
|
|
* will be on the diagonal. We can therfore just invert the
|
|
* diagonal at the end of the program to solve the equation system.
|
|
*/
|
|
for (i = 0; i < n; ++i) {
|
|
if (c[i + i * idem] == 0.0) {
|
|
/*
|
|
* Do a simple form of row pivoting to find a non-zero pivot
|
|
*/
|
|
for (k = i + 1; k < n; ++k) {
|
|
if (c[k + i * idem] != 0.0) {
|
|
goto FOUND_PIVOT;
|
|
}
|
|
}
|
|
#ifdef DEBUG_MODE
|
|
writelogf("vcs_mlequ ERROR: Encountered a zero column: %d\n", i);
|
|
#endif
|
|
return 1;
|
|
FOUND_PIVOT:
|
|
;
|
|
for (j = 0; j < n; ++j) {
|
|
c[i + j * idem] += c[k + j * idem];
|
|
}
|
|
for (j = 0; j < m; ++j) {
|
|
b[i + j * idem] += b[k + j * idem];
|
|
}
|
|
}
|
|
|
|
for (l = 0; l < n; ++l) {
|
|
if (l != i && c[l + i * idem] != 0.0) {
|
|
R = c[l + i * idem] / c[i + i * idem];
|
|
c[l + i * idem] = 0.0;
|
|
for (j = i+1; j < n; ++j) {
|
|
c[l + j * idem] -= c[i + j * idem] * R;
|
|
}
|
|
for (j = 0; j < m; ++j) {
|
|
b[l + j * idem] -= b[i + j * idem] * R;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* The negative in the last expression is due to the form of B upon
|
|
* input
|
|
*/
|
|
for (i = 0; i < n; ++i) {
|
|
for (j = 0; j < m; ++j) {
|
|
b[i + j * idem] = -b[i + j * idem] / c[i + i*idem];
|
|
}
|
|
}
|
|
return 0;
|
|
} /* mlequ() *************************************************************/
|
|
|
|
|
|
/*
|
|
*
|
|
* ElemRearrange:
|
|
*
|
|
* This subroutine handles the rearrangement of the constraint
|
|
* equations represented by the Formula Matrix. Rearrangement is only
|
|
* necessary when the number of components is less than the number of
|
|
* elements. For this case, some constraints can never be satisfied
|
|
* exactly, because the range space represented by the Formula
|
|
* Matrix of the components can't span the extra space. These
|
|
* constraints, which are out of the range space of the component
|
|
* Formula matrix entries, are migrated to the back of the Formula
|
|
* matrix.
|
|
*
|
|
* A prototypical example is an extra element column in
|
|
* FormulaMatrix[],
|
|
* which is identically zero. For example, let's say that argon is
|
|
* has an element column in FormulaMatrix[], but no species in the
|
|
* mechanism
|
|
* actually contains argon. Then, nc < ne. Unless the entry for
|
|
* desired element abundance vector for Ar is zero, then this
|
|
* element abundance constraint can never be satisfied. The
|
|
* constraint vector is not in the range space of the formula
|
|
* matrix.
|
|
* Also, without perturbation
|
|
* of FormulaMatrix[], BasisOptimize[] would produce a zero pivot
|
|
* because the matrix
|
|
* would be singular (unless the argon element column was already the
|
|
* last column of FormulaMatrix[].
|
|
* This routine borrows heavily from BasisOptimize algorithm. It
|
|
* finds nc constraints which span the range space of the Component
|
|
* Formula matrix, and assigns them as the first nc components in the
|
|
* formula matrix. This guarantees that BasisOptimize has a
|
|
* nonsingular matrix to invert.
|
|
*/
|
|
size_t Cantera::ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
|
|
MultiPhase* mphase,
|
|
std::vector<size_t>& orderVectorSpecies,
|
|
std::vector<size_t>& orderVectorElements)
|
|
{
|
|
|
|
size_t j, k, l, i, jl, ml, jr, ielem, jj, kk=0;
|
|
|
|
bool lindep = false;
|
|
size_t nelements = mphase->nElements();
|
|
std::string ename;
|
|
/*
|
|
* Get the total number of species in the multiphase object
|
|
*/
|
|
size_t nspecies = mphase->nSpecies();
|
|
|
|
double test = -1.0E10;
|
|
#ifdef DEBUG_MODE
|
|
if (BasisOptimize_print_lvl > 0) {
|
|
writelog(" ");
|
|
for (i=0; i<77; i++) {
|
|
writelog("-");
|
|
}
|
|
writelog("\n");
|
|
writelog(" --- Subroutine ElemRearrange() called to ");
|
|
writelog("check stoich. coefficient matrix\n");
|
|
writelog(" --- and to rearrange the element ordering once\n");
|
|
}
|
|
#endif
|
|
|
|
/*
|
|
* Perhaps, initialize the element ordering
|
|
*/
|
|
if (orderVectorElements.size() < nelements) {
|
|
orderVectorElements.resize(nelements);
|
|
for (j = 0; j < nelements; j++) {
|
|
orderVectorElements[j] = j;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Perhaps, initialize the species ordering. However, this is
|
|
* dangerous, as this ordering is assumed to yield the
|
|
* component species for the problem
|
|
*/
|
|
if (orderVectorSpecies.size() != nspecies) {
|
|
orderVectorSpecies.resize(nspecies);
|
|
for (k = 0; k < nspecies; k++) {
|
|
orderVectorSpecies[k] = k;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* If the elementAbundances aren't input, just create a fake one
|
|
* based on summing the column of the stoich matrix.
|
|
* This will force elements with zero species to the
|
|
* end of the element ordering.
|
|
*/
|
|
vector_fp eAbund(nelements,0.0);
|
|
if (elementAbundances.size() != nelements) {
|
|
for (j = 0; j < nelements; j++) {
|
|
eAbund[j] = 0.0;
|
|
for (k = 0; k < nspecies; k++) {
|
|
eAbund[j] += fabs(mphase->nAtoms(k, j));
|
|
}
|
|
}
|
|
} else {
|
|
copy(elementAbundances.begin(), elementAbundances.end(),
|
|
eAbund.begin());
|
|
}
|
|
|
|
vector_fp sa(nelements,0.0);
|
|
vector_fp ss(nelements,0.0);
|
|
vector_fp sm(nelements*nelements,0.0);
|
|
|
|
/*
|
|
* Top of a loop of some sort based on the index JR. JR is the
|
|
* current number independent elements found.
|
|
*/
|
|
jr = npos;
|
|
do {
|
|
++jr;
|
|
/*
|
|
* Top of another loop point based on finding a linearly
|
|
* independent element
|
|
*/
|
|
do {
|
|
/*
|
|
* Search the element vector. We first locate elements that
|
|
* are present in any amount. Then, we locate elements that
|
|
* are not present in any amount.
|
|
* Return its identity in K.
|
|
*/
|
|
k = nelements;
|
|
for (ielem = jr; ielem < nelements; ielem++) {
|
|
kk = orderVectorElements[ielem];
|
|
if (eAbund[kk] != test && eAbund[kk] > 0.0) {
|
|
k = ielem;
|
|
break;
|
|
}
|
|
}
|
|
for (ielem = jr; ielem < nelements; ielem++) {
|
|
kk = orderVectorElements[ielem];
|
|
if (eAbund[kk] != test) {
|
|
k = ielem;
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (k == nelements) {
|
|
// When we are here, there is an error usually.
|
|
// We haven't found the number of elements necessary.
|
|
// This is signalled by returning jr != nComponents.
|
|
#ifdef DEBUG_MODE
|
|
if (BasisOptimize_print_lvl > 0) {
|
|
writelogf("Error exit: returning with nComponents = %d\n", jr);
|
|
}
|
|
#endif
|
|
return jr;
|
|
}
|
|
|
|
/*
|
|
* Assign a large negative number to the element that we have
|
|
* just found, in order to take it out of further consideration.
|
|
*/
|
|
eAbund[kk] = test;
|
|
|
|
/* *********************************************************** */
|
|
/* **** CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX */
|
|
/* **** LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX ****** */
|
|
/* *********************************************************** */
|
|
/*
|
|
* Modified Gram-Schmidt Method, p. 202 Dalquist
|
|
* QR factorization of a matrix without row pivoting.
|
|
*/
|
|
jl = jr;
|
|
/*
|
|
* Fill in the row for the current element, k, under consideration
|
|
* The row will contain the Formula matrix value for that element
|
|
* with respect to the vector of component species.
|
|
* (note j and k indices are flipped compared to the previous routine)
|
|
*/
|
|
for (j = 0; j < nComponents; ++j) {
|
|
jj = orderVectorSpecies[j];
|
|
kk = orderVectorElements[k];
|
|
sm[j + jr*nComponents] = mphase->nAtoms(jj,kk);
|
|
}
|
|
if (jl > 0) {
|
|
/*
|
|
* Compute the coefficients of JA column of the
|
|
* the upper triangular R matrix, SS(J) = R_J_JR
|
|
* (this is slightly different than Dalquist)
|
|
* R_JA_JA = 1
|
|
*/
|
|
for (j = 0; j < jl; ++j) {
|
|
ss[j] = 0.0;
|
|
for (i = 0; i < nComponents; ++i) {
|
|
ss[j] += sm[i + jr*nComponents] * sm[i + j*nComponents];
|
|
}
|
|
ss[j] /= sa[j];
|
|
}
|
|
/*
|
|
* Now make the new column, (*,JR), orthogonal to the
|
|
* previous columns
|
|
*/
|
|
for (j = 0; j < jl; ++j) {
|
|
for (l = 0; l < nComponents; ++l) {
|
|
sm[l + jr*nComponents] -= ss[j] * sm[l + j*nComponents];
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Find the new length of the new column in Q.
|
|
* It will be used in the denominator in future row calcs.
|
|
*/
|
|
sa[jr] = 0.0;
|
|
for (ml = 0; ml < nComponents; ++ml) {
|
|
double tmp = sm[ml + jr*nComponents];
|
|
sa[jr] += tmp * tmp;
|
|
}
|
|
/* **************************************************** */
|
|
/* **** IF NORM OF NEW ROW .LT. 1E-6 REJECT ********** */
|
|
/* **************************************************** */
|
|
if (sa[jr] < 1.0e-6) {
|
|
lindep = true;
|
|
} else {
|
|
lindep = false;
|
|
}
|
|
} while (lindep);
|
|
/* ****************************************** */
|
|
/* **** REARRANGE THE DATA ****************** */
|
|
/* ****************************************** */
|
|
if (jr != k) {
|
|
#ifdef DEBUG_MODE
|
|
if (BasisOptimize_print_lvl > 0) {
|
|
kk = orderVectorElements[k];
|
|
ename = mphase->elementName(kk);
|
|
writelog(" --- ");
|
|
writelogf("%-2.2s", ename.c_str());
|
|
writelog("replaces ");
|
|
kk = orderVectorElements[jr];
|
|
ename = mphase->elementName(kk);
|
|
writelogf("%-2.2s", ename.c_str());
|
|
writelogf(" as element %3d\n", jr);
|
|
}
|
|
#endif
|
|
std::swap(orderVectorElements[jr], orderVectorElements[k]);
|
|
}
|
|
|
|
/*
|
|
* If we haven't found enough components, go back
|
|
* and find some more. (nc -1 is used below, because
|
|
* jr is counted from 0, via the C convention.
|
|
*/
|
|
} while (jr < (nComponents-1));
|
|
return nComponents;
|
|
} /* vcs_elem_rearrange() ****************************************************/
|
|
|