Cantera  2.3.0
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 http://www.cantera.org/license.txt for license and copyright information.
8 
10 
11 using namespace std;
12 
13 namespace Cantera
14 {
16 static const double USEDBEFORE = -1;
17 
18 size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
19  std::vector<size_t>& orderVectorSpecies,
20  std::vector<size_t>& orderVectorElements,
21  vector_fp& 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");
47  if (BasisOptimize_print_lvl > 0) {
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_fp molNum(nspecies,0.0);
83  mphase->getMoles(molNum.data());
84 
85  // Other workspace
86  DenseMatrix sm(ne, ne);
87  vector_fp ss(ne, 0.0);
88  vector_fp 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_fp 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  // i.e. 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 
292 void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
293  MultiPhase* mphase,
294  std::vector<size_t>& orderVectorSpecies,
295  std::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_fp 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_fp sa(nelements,0.0);
343  vector_fp ss(nelements,0.0);
344  vector_fp 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 }
size_t BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, MultiPhase *mphase, std::vector< size_t > &orderVectorSpecies, std::vector< size_t > &orderVectorElements, vector_fp &formRxnMatrix)
Choose the optimum basis of species for the equilibrium calculations.
int BasisOptimize_print_lvl
External int that is used to turn on debug printing for the BasisOptimze program. ...
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:179
size_t nSpecies() const
Number of species, summed over all phases.
Definition: MultiPhase.h:124
STL namespace.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
void ElemRearrange(size_t nComponents, const vector_fp &elementAbundances, MultiPhase *mphase, std::vector< size_t > &orderVectorSpecies, std::vector< size_t > &orderVectorElements)
Handles the potential rearrangement of the constraint equations represented by the Formula Matrix...
Headers for the MultiPhase object that is used to set up multiphase equilibrium problems (see Equilfu...
A class for multiphase mixtures.
Definition: MultiPhase.h:57
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
Definition: DenseMatrix.cpp:69
doublereal nAtoms(const size_t kGlob, const size_t mGlob) const
Returns the Number of atoms of global element mGlob in global species kGlob.
Definition: MultiPhase.cpp:764
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
vector_fp & data()
Return a reference to the data vector.
Definition: Array.h:299
void getMoles(doublereal *molNum) const
Get the mole numbers of all species in the multiphase object.
Definition: MultiPhase.cpp:370
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:200
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
std::string elementName(size_t m) const
Returns the name of the global element m.
Definition: MultiPhase.cpp:730
Namespace for the Cantera kernel.
Definition: application.cpp:29
size_t nElements() const
Number of elements.
Definition: MultiPhase.h:98
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition: DenseMatrix.h:72
std::string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
Definition: MultiPhase.cpp:759