16static const double USEDBEFORE = -1;
19 vector<size_t>& orderVectorSpecies,
20 vector<size_t>& orderVectorElements,
21 vector<double>& formRxnMatrix)
27 size_t nspecies = mphase->
nSpecies();
30 if (orderVectorElements.size() < ne) {
31 orderVectorElements.resize(ne);
32 iota(orderVectorElements.begin(), orderVectorElements.end(), 0);
36 if (orderVectorSpecies.size() != nspecies) {
37 orderVectorSpecies.resize(nspecies);
38 iota(orderVectorSpecies.begin(), orderVectorSpecies.end(), 0);
44 writelog(
" --- Subroutine BASOPT called to ");
45 writelog(
"calculate the number of components and ");
46 writelog(
"evaluate the formation matrix\n");
50 writelog(
" --- Formula Matrix used in BASOPT calculation\n");
52 for (
size_t j = 0; j < ne; j++) {
53 size_t jj = orderVectorElements[j];
57 for (
size_t k = 0; k < nspecies; k++) {
58 size_t kk = orderVectorSpecies[k];
59 writelog(
" --- {:>11.11s} | {:4d} |",
61 for (
size_t j = 0; j < ne; j++) {
62 size_t jj = orderVectorElements[j];
63 double num = mphase->
nAtoms(kk,jj);
75 size_t nComponents = std::min(ne, nspecies);
76 size_t nNonComponents = nspecies - nComponents;
79 *usedZeroedSpecies =
false;
82 vector<double> molNum(nspecies,0.0);
87 vector<double> ss(ne, 0.0);
88 vector<double> sa(ne, 0.0);
89 if (formRxnMatrix.size() < nspecies*ne) {
90 formRxnMatrix.resize(nspecies*ne, 0.0);
94 vector<double> molNumBase = molNum;
100 while (jr < nComponents) {
108 size_t kk = max_element(molNum.begin(), molNum.end()) - molNum.begin();
110 for (j = 0; j < nspecies; j++) {
111 if (orderVectorSpecies[j] == kk) {
117 throw CanteraError(
"BasisOptimize",
"orderVectorSpecies contains an error");
120 if (molNum[kk] == 0.0) {
121 *usedZeroedSpecies =
true;
124 if (molNum[kk] == USEDBEFORE) {
126 nNonComponents = nspecies - nComponents;
132 molSave = molNum[kk];
133 molNum[kk] = USEDBEFORE;
140 for (j = 0; j < ne; ++j) {
141 size_t jj = orderVectorElements[j];
142 sm(j, jr) = mphase->
nAtoms(kk,jj);
148 for (j = 0; j < jl; ++j) {
150 for (
size_t i = 0; i < ne; ++i) {
151 ss[j] += sm(i, jr) * sm(i, j);
158 for (j = 0; j < jl; ++j) {
159 for (
size_t i = 0; i < ne; ++i) {
160 sm(i, jr) -= ss[j] * sm(i, j);
168 for (
size_t ml = 0; ml < ne; ++ml) {
169 sa[jr] += pow(sm(ml, jr), 2);
173 if (sa[jr] > 1.0e-6) {
181 size_t kk = orderVectorSpecies[k];
183 size_t jj = orderVectorSpecies[jr];
186 writelogf(
"(%9.2g) as component %3d\n", molNum[jj], jr);
188 std::swap(orderVectorSpecies[jr], orderVectorSpecies[k]);
230 sm.
resize(nComponents, nComponents);
231 for (
size_t k = 0; k < nComponents; ++k) {
232 size_t kk = orderVectorSpecies[k];
233 for (
size_t j = 0; j < nComponents; ++j) {
234 size_t jj = orderVectorElements[j];
235 sm(j, k) = mphase->
nAtoms(kk, jj);
239 for (
size_t i = 0; i < nNonComponents; ++i) {
240 size_t k = nComponents + i;
241 size_t kk = orderVectorSpecies[k];
242 for (
size_t j = 0; j < nComponents; ++j) {
243 size_t jj = orderVectorElements[j];
244 formRxnMatrix[j + i * ne] = - mphase->
nAtoms(kk, jj);
249 solve(sm, formRxnMatrix.
data(), nNonComponents, ne);
253 writelogf(
" --- Number of Components = %d\n", nComponents);
256 for (
size_t k = 0; k < nComponents; k++) {
257 size_t kk = orderVectorSpecies[k];
260 writelog(
"\n --- Components Moles: ");
261 for (
size_t k = 0; k < nComponents; k++) {
262 size_t kk = orderVectorSpecies[k];
265 writelog(
"\n --- NonComponent | Moles | ");
266 for (
size_t i = 0; i < nComponents; i++) {
267 size_t kk = orderVectorSpecies[i];
272 for (
size_t i = 0; i < nNonComponents; i++) {
273 size_t k = i + nComponents;
274 size_t kk = orderVectorSpecies[k];
279 for (
size_t j = 0; j < nComponents; j++) {
280 writelogf(
" %6.2f", - formRxnMatrix[j + i * ne]);
292void ElemRearrange(
size_t nComponents,
const vector<double>& elementAbundances,
294 vector<size_t>& orderVectorSpecies,
295 vector<size_t>& orderVectorElements)
299 size_t nspecies = mphase->
nSpecies();
304 writelog(
" --- Subroutine ElemRearrange() called to ");
305 writelog(
"check stoich. coefficient matrix\n");
306 writelog(
" --- and to rearrange the element ordering once\n");
310 if (orderVectorElements.size() < nelements) {
311 orderVectorElements.resize(nelements);
312 for (
size_t j = 0; j < nelements; j++) {
313 orderVectorElements[j] = j;
319 if (orderVectorSpecies.size() != nspecies) {
320 orderVectorSpecies.resize(nspecies);
321 for (
size_t k = 0; k < nspecies; k++) {
322 orderVectorSpecies[k] = k;
329 vector<double> eAbund(nelements,0.0);
330 if (elementAbundances.size() != nelements) {
331 for (
size_t j = 0; j < nelements; j++) {
333 for (
size_t k = 0; k < nspecies; k++) {
334 eAbund[j] += fabs(mphase->
nAtoms(k, j));
338 copy(elementAbundances.begin(), elementAbundances.end(),
342 vector<double> sa(nelements,0.0);
343 vector<double> ss(nelements,0.0);
344 vector<double> sm(nelements*nelements,0.0);
349 while (jr < nComponents) {
352 size_t k = nelements;
359 for (
size_t ielem = jr; ielem < nelements; ielem++) {
360 kk = orderVectorElements[ielem];
361 if (eAbund[kk] != USEDBEFORE && eAbund[kk] > 0.0) {
366 for (
size_t ielem = jr; ielem < nelements; ielem++) {
367 kk = orderVectorElements[ielem];
368 if (eAbund[kk] != USEDBEFORE) {
374 if (k == nelements) {
378 writelogf(
"Error exit: returning with nComponents = %d\n", jr);
380 throw CanteraError(
"ElemRearrange",
"Required number of elements not found.");
385 eAbund[kk] = USEDBEFORE;
398 for (
size_t j = 0; j < nComponents; ++j) {
399 size_t jj = orderVectorSpecies[j];
400 kk = orderVectorElements[k];
401 sm[j + jr*nComponents] = mphase->
nAtoms(jj,kk);
407 for (
size_t j = 0; j < jl; ++j) {
409 for (
size_t i = 0; i < nComponents; ++i) {
410 ss[j] += sm[i + jr*nComponents] * sm[i + j*nComponents];
417 for (
size_t j = 0; j < jl; ++j) {
418 for (
size_t i = 0; i < nComponents; ++i) {
419 sm[i + jr*nComponents] -= ss[j] * sm[i + j*nComponents];
427 for (
size_t ml = 0; ml < nComponents; ++ml) {
428 double tmp = sm[ml + jr*nComponents];
432 if (sa[jr] > 1.0e-6) {
439 size_t kk = orderVectorElements[k];
443 kk = orderVectorElements[jr];
447 std::swap(orderVectorElements[jr], orderVectorElements[k]);
Headers for the MultiPhase object that is used to set up multiphase equilibrium problems (see Chemica...
vector< double > & data()
Return a reference to the data vector.
Base class for exceptions thrown by Cantera classes.
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
A class for multiphase mixtures.
double nAtoms(const size_t kGlob, const size_t mGlob) const
Returns the Number of atoms of global element mGlob in global species kGlob.
size_t nSpecies() const
Number of species, summed over all phases.
void getMoles(double *molNum) const
Get the mole numbers of all species in the multiphase object.
string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
size_t nElements() const
Number of elements.
string elementName(size_t m) const
Returns the name of the global element m.
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
void ElemRearrange(size_t nComponents, const vector< double > &elementAbundances, MultiPhase *mphase, vector< size_t > &orderVectorSpecies, vector< size_t > &orderVectorElements)
Handles the potential rearrangement of the constraint equations represented by the Formula Matrix.
size_t BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, MultiPhase *mphase, vector< size_t > &orderVectorSpecies, vector< size_t > &orderVectorElements, vector< double > &formRxnMatrix)
Choose the optimum basis of species for the equilibrium calculations.
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
int BasisOptimize_print_lvl
External int that is used to turn on debug printing for the BasisOptimize program.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...