Cantera  3.0.0
Loading...
Searching...
No Matches
BasisOptimize.cpp
Go to the documentation of this file.
1/**
2 * @file BasisOptimize.cpp Functions which calculation optimized basis of the
3 * stoichiometric coefficient matrix (see /ref equil functions)
4 */
5
6// This file is part of Cantera. See License.txt in the top-level directory or
7// at https://cantera.org/license.txt for license and copyright information.
8
11#include "cantera/base/global.h"
12
13namespace Cantera
14{
16static const double USEDBEFORE = -1;
17
18size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
19 vector<size_t>& orderVectorSpecies,
20 vector<size_t>& orderVectorElements,
21 vector<double>& formRxnMatrix)
22{
23 // Get the total number of elements defined in the multiphase object
24 size_t ne = mphase->nElements();
25
26 // Get the total number of species in the multiphase object
27 size_t nspecies = mphase->nSpecies();
28
29 // Perhaps, initialize the element ordering
30 if (orderVectorElements.size() < ne) {
31 orderVectorElements.resize(ne);
32 iota(orderVectorElements.begin(), orderVectorElements.end(), 0);
33 }
34
35 // Perhaps, initialize the species ordering
36 if (orderVectorSpecies.size() != nspecies) {
37 orderVectorSpecies.resize(nspecies);
38 iota(orderVectorSpecies.begin(), orderVectorSpecies.end(), 0);
39 }
40
41 if (BasisOptimize_print_lvl >= 1) {
42 writelog(" ");
43 writeline('-', 77);
44 writelog(" --- Subroutine BASOPT called to ");
45 writelog("calculate the number of components and ");
46 writelog("evaluate the formation matrix\n");
48 writelog(" ---\n");
49
50 writelog(" --- Formula Matrix used in BASOPT calculation\n");
51 writelog(" --- Species | Order | ");
52 for (size_t j = 0; j < ne; j++) {
53 size_t jj = orderVectorElements[j];
54 writelog(" {:>4.4s}({:1d})", mphase->elementName(jj), j);
55 }
56 writelog("\n");
57 for (size_t k = 0; k < nspecies; k++) {
58 size_t kk = orderVectorSpecies[k];
59 writelog(" --- {:>11.11s} | {:4d} |",
60 mphase->speciesName(kk), k);
61 for (size_t j = 0; j < ne; j++) {
62 size_t jj = orderVectorElements[j];
63 double num = mphase->nAtoms(kk,jj);
64 writelogf("%6.1g ", num);
65 }
66 writelog("\n");
67 }
68 writelog(" --- \n");
69 }
70 }
71
72 // Calculate the maximum value of the number of components possible. It's
73 // equal to the minimum of the number of elements and the number of total
74 // species.
75 size_t nComponents = std::min(ne, nspecies);
76 size_t nNonComponents = nspecies - nComponents;
77
78 // Set this return variable to false
79 *usedZeroedSpecies = false;
80
81 // Create an array of mole numbers
82 vector<double> molNum(nspecies,0.0);
83 mphase->getMoles(molNum.data());
84
85 // Other workspace
86 DenseMatrix sm(ne, ne);
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);
91 }
92
93 // For debugging purposes keep an unmodified copy of the array.
94 vector<double> molNumBase = molNum;
95 double molSave = 0.0;
96 size_t jr = 0;
97
98 // Top of a loop of some sort based on the index JR. JR is the current
99 // number of component species found.
100 while (jr < nComponents) {
101 // Top of another loop point based on finding a linearly independent
102 // species
103 size_t k = npos;
104 while (true) {
105 // Search the remaining part of the mole number vector, molNum for
106 // the largest remaining species. Return its identity. kk is the raw
107 // number. k is the orderVectorSpecies index.
108 size_t kk = max_element(molNum.begin(), molNum.end()) - molNum.begin();
109 size_t j;
110 for (j = 0; j < nspecies; j++) {
111 if (orderVectorSpecies[j] == kk) {
112 k = j;
113 break;
114 }
115 }
116 if (j == nspecies) {
117 throw CanteraError("BasisOptimize", "orderVectorSpecies contains an error");
118 }
119
120 if (molNum[kk] == 0.0) {
121 *usedZeroedSpecies = true;
122 }
123 // If the largest molNum is negative, then we are done.
124 if (molNum[kk] == USEDBEFORE) {
125 nComponents = jr;
126 nNonComponents = nspecies - nComponents;
127 break;
128 }
129
130 // Assign a small negative number to the component that we have
131 // just found, in order to take it out of further consideration.
132 molSave = molNum[kk];
133 molNum[kk] = USEDBEFORE;
134
135 // CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES
136
137 // Modified Gram-Schmidt Method, p. 202 Dalquist
138 // QR factorization of a matrix without row pivoting.
139 size_t jl = jr;
140 for (j = 0; j < ne; ++j) {
141 size_t jj = orderVectorElements[j];
142 sm(j, jr) = mphase->nAtoms(kk,jj);
143 }
144 if (jl > 0) {
145 // Compute the coefficients of JA column of the the upper
146 // triangular R matrix, SS(J) = R_J_JR (this is slightly
147 // different than Dalquist) R_JA_JA = 1
148 for (j = 0; j < jl; ++j) {
149 ss[j] = 0.0;
150 for (size_t i = 0; i < ne; ++i) {
151 ss[j] += sm(i, jr) * sm(i, j);
152 }
153 ss[j] /= sa[j];
154 }
155
156 // Now make the new column, (*,JR), orthogonal to the previous
157 // columns
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);
161 }
162 }
163 }
164
165 // Find the new length of the new column in Q.
166 // It will be used in the denominator in future row calcs.
167 sa[jr] = 0.0;
168 for (size_t ml = 0; ml < ne; ++ml) {
169 sa[jr] += pow(sm(ml, jr), 2);
170 }
171
172 // IF NORM OF NEW ROW .LT. 1E-3 REJECT
173 if (sa[jr] > 1.0e-6) {
174 break;
175 }
176 }
177
178 // REARRANGE THE DATA
179 if (jr != k) {
180 if (BasisOptimize_print_lvl >= 1) {
181 size_t kk = orderVectorSpecies[k];
182 writelogf(" --- %-12.12s", mphase->speciesName(kk));
183 size_t jj = orderVectorSpecies[jr];
184 writelogf("(%9.2g) replaces %-12.12s",
185 molSave, mphase->speciesName(jj));
186 writelogf("(%9.2g) as component %3d\n", molNum[jj], jr);
187 }
188 std::swap(orderVectorSpecies[jr], orderVectorSpecies[k]);
189 }
190
191 // If we haven't found enough components, go back and find some more
192 jr++;
193 }
194
195 if (! doFormRxn) {
196 return nComponents;
197 }
198
199 // EVALUATE THE STOICHIOMETRY
200 //
201 // Formulate the matrix problem for the stoichiometric
202 // coefficients. CX + B = 0
203 //
204 // C will be an nc x nc matrix made up of the formula vectors for the
205 // components. Each component's formula vector is a column. The rows are the
206 // elements.
207 //
208 // n RHS's will be solved for. Thus, B is an nc x n matrix.
209 //
210 // BIG PROBLEM 1/21/99:
211 //
212 // This algorithm makes the assumption that the first nc rows of the formula
213 // matrix aren't rank deficient. However, this might not be the case. For
214 // example, assume that the first element in FormulaMatrix[] is argon.
215 // Assume that no species in the matrix problem actually includes argon.
216 // Then, the first row in sm[], below will be identically zero. bleh.
217 //
218 // What needs to be done is to perform a rearrangement of the ELEMENTS ->
219 // that is, rearrange, FormulaMatrix, sp, and gai, such that the first nc
220 // elements form in combination with the nc components create an invertible
221 // sm[]. not a small project, but very doable.
222 //
223 // An alternative would be to turn the matrix problem below into an ne x nc
224 // problem, and do QR elimination instead of Gauss-Jordan elimination.
225 //
226 // Note the rearrangement of elements need only be done once in the problem.
227 // It's actually very similar to the top of this program with ne being the
228 // species and nc being the elements!!
229
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);
236 }
237 }
238
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);
245 }
246 }
247
248 // // Use LU factorization to calculate the reaction matrix
249 solve(sm, formRxnMatrix.data(), nNonComponents, ne);
250
251 if (BasisOptimize_print_lvl >= 1) {
252 writelog(" ---\n");
253 writelogf(" --- Number of Components = %d\n", nComponents);
254 writelog(" --- Formula Matrix:\n");
255 writelog(" --- Components: ");
256 for (size_t k = 0; k < nComponents; k++) {
257 size_t kk = orderVectorSpecies[k];
258 writelogf(" %3d (%3d) ", k, kk);
259 }
260 writelog("\n --- Components Moles: ");
261 for (size_t k = 0; k < nComponents; k++) {
262 size_t kk = orderVectorSpecies[k];
263 writelogf("%-11.3g", molNumBase[kk]);
264 }
265 writelog("\n --- NonComponent | Moles | ");
266 for (size_t i = 0; i < nComponents; i++) {
267 size_t kk = orderVectorSpecies[i];
268 writelogf("%-11.10s", mphase->speciesName(kk));
269 }
270 writelog("\n");
271
272 for (size_t i = 0; i < nNonComponents; i++) {
273 size_t k = i + nComponents;
274 size_t kk = orderVectorSpecies[k];
275 writelogf(" --- %3d (%3d) ", k, kk);
276 writelogf("%-10.10s", mphase->speciesName(kk));
277 writelogf("|%10.3g|", molNumBase[kk]);
278 // Print the negative of formRxnMatrix[]; it's easier to interpret.
279 for (size_t j = 0; j < nComponents; j++) {
280 writelogf(" %6.2f", - formRxnMatrix[j + i * ne]);
281 }
282 writelog("\n");
283 }
284 writelog(" ");
285 writeline('-', 77);
286 }
287
288 return nComponents;
289} // basopt()
290
291
292void ElemRearrange(size_t nComponents, const vector<double>& elementAbundances,
293 MultiPhase* mphase,
294 vector<size_t>& orderVectorSpecies,
295 vector<size_t>& orderVectorElements)
296{
297 size_t nelements = mphase->nElements();
298 // Get the total number of species in the multiphase object
299 size_t nspecies = mphase->nSpecies();
300
301 if (BasisOptimize_print_lvl > 0) {
302 writelog(" ");
303 writeline('-', 77);
304 writelog(" --- Subroutine ElemRearrange() called to ");
305 writelog("check stoich. coefficient matrix\n");
306 writelog(" --- and to rearrange the element ordering once\n");
307 }
308
309 // Perhaps, initialize the element ordering
310 if (orderVectorElements.size() < nelements) {
311 orderVectorElements.resize(nelements);
312 for (size_t j = 0; j < nelements; j++) {
313 orderVectorElements[j] = j;
314 }
315 }
316
317 // Perhaps, initialize the species ordering. However, this is dangerous, as
318 // this ordering is assumed to yield the component species for the problem
319 if (orderVectorSpecies.size() != nspecies) {
320 orderVectorSpecies.resize(nspecies);
321 for (size_t k = 0; k < nspecies; k++) {
322 orderVectorSpecies[k] = k;
323 }
324 }
325
326 // If the elementAbundances aren't input, just create a fake one based on
327 // summing the column of the stoich matrix. This will force elements with
328 // zero species to the end of the element ordering.
329 vector<double> eAbund(nelements,0.0);
330 if (elementAbundances.size() != nelements) {
331 for (size_t j = 0; j < nelements; j++) {
332 eAbund[j] = 0.0;
333 for (size_t k = 0; k < nspecies; k++) {
334 eAbund[j] += fabs(mphase->nAtoms(k, j));
335 }
336 }
337 } else {
338 copy(elementAbundances.begin(), elementAbundances.end(),
339 eAbund.begin());
340 }
341
342 vector<double> sa(nelements,0.0);
343 vector<double> ss(nelements,0.0);
344 vector<double> sm(nelements*nelements,0.0);
345
346 // Top of a loop of some sort based on the index JR. JR is the current
347 // number independent elements found.
348 size_t jr = 0;
349 while (jr < nComponents) {
350 // Top of another loop point based on finding a linearly independent
351 // element
352 size_t k = nelements;
353 while (true) {
354 // Search the element vector. We first locate elements that are
355 // present in any amount. Then, we locate elements that are not
356 // present in any amount. Return its identity in K.
357 k = nelements;
358 size_t kk;
359 for (size_t ielem = jr; ielem < nelements; ielem++) {
360 kk = orderVectorElements[ielem];
361 if (eAbund[kk] != USEDBEFORE && eAbund[kk] > 0.0) {
362 k = ielem;
363 break;
364 }
365 }
366 for (size_t ielem = jr; ielem < nelements; ielem++) {
367 kk = orderVectorElements[ielem];
368 if (eAbund[kk] != USEDBEFORE) {
369 k = ielem;
370 break;
371 }
372 }
373
374 if (k == nelements) {
375 // When we are here, there is an error usually.
376 // We haven't found the number of elements necessary.
377 if (BasisOptimize_print_lvl > 0) {
378 writelogf("Error exit: returning with nComponents = %d\n", jr);
379 }
380 throw CanteraError("ElemRearrange", "Required number of elements not found.");
381 }
382
383 // Assign a large negative number to the element that we have
384 // just found, in order to take it out of further consideration.
385 eAbund[kk] = USEDBEFORE;
386
387 // CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX
388 // LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX
389
390 // Modified Gram-Schmidt Method, p. 202 Dalquist
391 // QR factorization of a matrix without row pivoting.
392 size_t jl = jr;
393
394 // Fill in the row for the current element, k, under consideration
395 // The row will contain the Formula matrix value for that element
396 // with respect to the vector of component species. (note j and k
397 // indices are flipped compared to the previous routine)
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);
402 }
403 if (jl > 0) {
404 // Compute the coefficients of JA column of the the upper
405 // triangular R matrix, SS(J) = R_J_JR (this is slightly
406 // different than Dalquist) R_JA_JA = 1
407 for (size_t j = 0; j < jl; ++j) {
408 ss[j] = 0.0;
409 for (size_t i = 0; i < nComponents; ++i) {
410 ss[j] += sm[i + jr*nComponents] * sm[i + j*nComponents];
411 }
412 ss[j] /= sa[j];
413 }
414
415 // Now make the new column, (*,JR), orthogonal to the
416 // previous columns
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];
420 }
421 }
422 }
423
424 // Find the new length of the new column in Q.
425 // It will be used in the denominator in future row calcs.
426 sa[jr] = 0.0;
427 for (size_t ml = 0; ml < nComponents; ++ml) {
428 double tmp = sm[ml + jr*nComponents];
429 sa[jr] += tmp * tmp;
430 }
431 // IF NORM OF NEW ROW .LT. 1E-6 REJECT
432 if (sa[jr] > 1.0e-6) {
433 break;
434 }
435 }
436 // REARRANGE THE DATA
437 if (jr != k) {
438 if (BasisOptimize_print_lvl > 0) {
439 size_t kk = orderVectorElements[k];
440 writelog(" --- ");
441 writelogf("%-2.2s", mphase->elementName(kk));
442 writelog("replaces ");
443 kk = orderVectorElements[jr];
444 writelogf("%-2.2s", mphase->elementName(kk));
445 writelogf(" as element %3d\n", jr);
446 }
447 std::swap(orderVectorElements[jr], orderVectorElements[k]);
448 }
449
450 // If we haven't found enough components, go back and find some more
451 jr++;
452 };
453}
454
455}
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.
Definition Array.h:218
Base class for exceptions thrown by Cantera classes.
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition DenseMatrix.h:55
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
A class for multiphase mixtures.
Definition MultiPhase.h:57
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.
Definition MultiPhase.h:123
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.
Definition MultiPhase.h:97
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.
Definition global.h:195
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:175
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition ct_defs.h:195
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...