Cantera  2.2.1
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  */
7
8 using namespace Cantera;
9 using namespace std;
10
11 namespace Cantera
12 {
14 }
15
16 //! Print a string within a given space limit.
17 /*!
18  * This routine limits the amount of the string that will be printed to a
19  * maximum of "space" characters.
20  * @param str String -> must be null terminated.
21  * @param space space limit for the printing.
22  * @param alignment 0 centered
23  * 1 right aligned
24  * 2 left aligned
25  */
26 static void print_stringTrunc(const char* str, int space, int alignment);
27
28 size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn,
29  MultiPhase* mphase, std::vector<size_t>& orderVectorSpecies,
30  std::vector<size_t>& orderVectorElements,
31  vector_fp& formRxnMatrix)
32 {
33  size_t j, jj, k=0, kk, l, i, jl, ml;
34  bool lindep;
35  std::string ename;
36  std::string sname;
37  /*
38  * Get the total number of elements defined in the multiphase object
39  */
40  size_t ne = mphase->nElements();
41  /*
42  * Get the total number of species in the multiphase object
43  */
44  size_t nspecies = mphase->nSpecies();
45  doublereal tmp;
46  doublereal const USEDBEFORE = -1;
47
48  /*
49  * Perhaps, initialize the element ordering
50  */
51  if (orderVectorElements.size() < ne) {
52  orderVectorElements.resize(ne);
53  for (j = 0; j < ne; j++) {
54  orderVectorElements[j] = j;
55  }
56  }
57
58  /*
59  * Perhaps, initialize the species ordering
60  */
61  if (orderVectorSpecies.size() != nspecies) {
62  orderVectorSpecies.resize(nspecies);
63  for (k = 0; k < nspecies; k++) {
64  orderVectorSpecies[k] = k;
65  }
66  }
67
68  if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl >= 1) {
69  writelog(" ");
70  for (i=0; i<77; i++) {
71  writelog("-");
72  }
73  writelog("\n");
74  writelog(" --- Subroutine BASOPT called to ");
75  writelog("calculate the number of components and ");
76  writelog("evaluate the formation matrix\n");
77  if (BasisOptimize_print_lvl > 0) {
78  writelog(" ---\n");
79
80  writelog(" --- Formula Matrix used in BASOPT calculation\n");
81  writelog(" --- Species | Order | ");
82  for (j = 0; j < ne; j++) {
83  jj = orderVectorElements[j];
84  writelog(" ");
85  ename = mphase->elementName(jj);
86  print_stringTrunc(ename.c_str(), 4, 1);
87  writelogf("(%1d)", j);
88  }
89  writelog("\n");
90  for (k = 0; k < nspecies; k++) {
91  kk = orderVectorSpecies[k];
92  writelog(" --- ");
93  sname = mphase->speciesName(kk);
94  print_stringTrunc(sname.c_str(), 11, 1);
95  writelogf(" | %4d |", k);
96  for (j = 0; j < ne; j++) {
97  jj = orderVectorElements[j];
98  double num = mphase->nAtoms(kk,jj);
99  writelogf("%6.1g ", num);
100  }
101  writelog("\n");
102  }
103  writelog(" --- \n");
104  }
105  }
106
107  /*
108  * Calculate the maximum value of the number of components possible
109  * It's equal to the minimum of the number of elements and the
110  * number of total species.
111  */
112  size_t nComponents = std::min(ne, nspecies);
113  size_t nNonComponents = nspecies - nComponents;
114  /*
115  * Set this return variable to false
116  */
117  *usedZeroedSpecies = false;
118
119  /*
120  * Create an array of mole numbers
121  */
122  vector_fp molNum(nspecies,0.0);
123  mphase->getMoles(DATA_PTR(molNum));
124
125  /*
126  * Other workspace
127  */
128  vector_fp sm(ne*ne, 0.0);
129  vector_fp ss(ne, 0.0);
130  vector_fp sa(ne, 0.0);
131  if (formRxnMatrix.size() < nspecies*ne) {
132  formRxnMatrix.resize(nspecies*ne, 0.0);
133  }
134
135  /*
136  * For debugging purposes keep an unmodified copy of the array.
137  */
138  vector_fp molNumBase;
139  if (DEBUG_MODE_ENABLED) {
140  molNumBase = molNum;
141  }
142  double molSave = 0.0;
143
144  size_t jr = npos;
145  /*
146  * Top of a loop of some sort based on the index JR. JR is the
147  * current number of component species found.
148  */
149  do {
150  ++jr;
151  /* - Top of another loop point based on finding a linearly */
152  /* - independent species */
153  do {
154  /*
155  * Search the remaining part of the mole number vector, molNum
156  * for the largest remaining species. Return its identity.
157  * kk is the raw number. k is the orderVectorSpecies index.
158  */
159  kk = max_element(molNum.begin(), molNum.end()) - molNum.begin();
160  for (j = 0; j < nspecies; j++) {
161  if (orderVectorSpecies[j] == kk) {
162  k = j;
163  break;
164  }
165  }
166  if (j == nspecies) {
167  throw CanteraError("BasisOptimize", "orderVectorSpecies contains an error");
168  }
169
170  if (molNum[kk] == 0.0) {
171  *usedZeroedSpecies = true;
172  }
173  /*
174  * If the largest molNum is negative, then we are done.
175  */
176  if (molNum[kk] == USEDBEFORE) {
177  nComponents = jr;
178  nNonComponents = nspecies - nComponents;
179  break;
180  }
181  /*
182  * Assign a small negative number to the component that we have
183  * just found, in order to take it out of further consideration.
184  */
185 #ifdef DEBUG_MODE
186  molSave = molNum[kk];
187 #endif
188  molNum[kk] = USEDBEFORE;
189
190  /* *********************************************************** */
191  /* **** CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES ****** */
192  /* *********************************************************** */
193  /*
194  * Modified Gram-Schmidt Method, p. 202 Dalquist
195  * QR factorization of a matrix without row pivoting.
196  */
197  jl = jr;
198  for (j = 0; j < ne; ++j) {
199  jj = orderVectorElements[j];
200  sm[j + jr*ne] = mphase->nAtoms(kk,jj);
201  }
202  if (jl > 0) {
203  /*
204  * Compute the coefficients of JA column of the
205  * the upper triangular R matrix, SS(J) = R_J_JR
206  * (this is slightly different than Dalquist)
207  * R_JA_JA = 1
208  */
209  for (j = 0; j < jl; ++j) {
210  ss[j] = 0.0;
211  for (i = 0; i < ne; ++i) {
212  ss[j] += sm[i + jr*ne] * sm[i + j*ne];
213  }
214  ss[j] /= sa[j];
215  }
216  /*
217  * Now make the new column, (*,JR), orthogonal to the
218  * previous columns
219  */
220  for (j = 0; j < jl; ++j) {
221  for (l = 0; l < ne; ++l) {
222  sm[l + jr*ne] -= ss[j] * sm[l + j*ne];
223  }
224  }
225  }
226  /*
227  * Find the new length of the new column in Q.
228  * It will be used in the denominator in future row calcs.
229  */
230  sa[jr] = 0.0;
231  for (ml = 0; ml < ne; ++ml) {
232  tmp = sm[ml + jr*ne];
233  sa[jr] += tmp * tmp;
234  }
235  /* **************************************************** */
236  /* **** IF NORM OF NEW ROW .LT. 1E-3 REJECT ********** */
237  /* **************************************************** */
238  if (sa[jr] < 1.0e-6) {
239  lindep = true;
240  } else {
241  lindep = false;
242  }
243  } while (lindep);
244  /* ****************************************** */
245  /* **** REARRANGE THE DATA ****************** */
246  /* ****************************************** */
247  if (jr != k) {
248  if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl >= 1) {
249  kk = orderVectorSpecies[k];
250  sname = mphase->speciesName(kk);
251  writelogf(" --- %-12.12s", sname.c_str());
252  jj = orderVectorSpecies[jr];
253  ename = mphase->speciesName(jj);
254  writelogf("(%9.2g) replaces %-12.12s", molSave, ename.c_str());
255  writelogf("(%9.2g) as component %3d\n", molNum[jj], jr);
256  }
257  std::swap(orderVectorSpecies[jr], orderVectorSpecies[k]);
258  }
259
260  /*
261  * If we haven't found enough components, go back
262  * and find some more. (nc -1 is used below, because
263  * jr is counted from 0, via the C convention.
264  */
265  } while (jr < (nComponents-1));
266
267
268  if (! doFormRxn) {
269  return nComponents;
270  }
271
272  /* ****************************************************** */
273  /* **** EVALUATE THE STOICHIOMETRY ********************** */
274  /* ****************************************************** */
275  /*
276  * Formulate the matrix problem for the stoichiometric
277  * coefficients. CX + B = 0
278  * C will be an nc x nc matrix made up of the formula
279  * vectors for the components. Each component's formula
280  * vector is a column. The rows are the elements.
281  * n RHS's will be solved for. Thus, B is an nc x n
282  * matrix.
283  *
284  * BIG PROBLEM 1/21/99:
285  *
286  * This algorithm makes the assumption that the
287  * first nc rows of the formula matrix aren't rank deficient.
288  * However, this might not be the case. For example, assume
289  * that the first element in FormulaMatrix[] is argon. Assume that
290  * no species in the matrix problem actually includes argon.
291  * Then, the first row in sm[], below will be identically
292  * zero. bleh.
293  * What needs to be done is to perform a rearrangement
294  * of the ELEMENTS -> i.e. rearrange, FormulaMatrix, sp, and gai, such
295  * that the first nc elements form in combination with the
296  * nc components create an invertible sm[]. not a small
297  * project, but very doable.
298  * An alternative would be to turn the matrix problem
299  * below into an ne x nc problem, and do QR elimination instead
300  * of Gauss-Jordan elimination.
301  * Note the rearrangement of elements need only be done once
302  * in the problem. It's actually very similar to the top of
303  * this program with ne being the species and nc being the
304  * elements!!
305  */
306  for (k = 0; k < nComponents; ++k) {
307  kk = orderVectorSpecies[k];
308  for (j = 0; j < nComponents; ++j) {
309  jj = orderVectorElements[j];
310  sm[j + k*ne] = mphase->nAtoms(kk, jj);
311  }
312  }
313
314  for (i = 0; i < nNonComponents; ++i) {
315  k = nComponents + i;
316  kk = orderVectorSpecies[k];
317  for (j = 0; j < nComponents; ++j) {
318  jj = orderVectorElements[j];
319  formRxnMatrix[j + i * ne] = - mphase->nAtoms(kk, jj);
320  }
321  }
322  // Use LU factorization to calculate the reaction matrix
323  int info;
324  vector_int ipiv(nComponents);
325  ct_dgetrf(nComponents, nComponents, &sm[0], ne, &ipiv[0], info);
326  if (info) {
327  throw CanteraError("BasisOptimize", "factorization returned an error condition");
328  }
329  ct_dgetrs(ctlapack::NoTranspose, nComponents, nNonComponents, &sm[0], ne,
330  &ipiv[0], &formRxnMatrix[0], ne, info);
331
332  if (DEBUG_MODE_ENABLED && Cantera::BasisOptimize_print_lvl >= 1) {
333  writelog(" ---\n");
334  writelogf(" --- Number of Components = %d\n", nComponents);
335  writelog(" --- Formula Matrix:\n");
336  writelog(" --- Components: ");
337  for (k = 0; k < nComponents; k++) {
338  kk = orderVectorSpecies[k];
339  writelogf(" %3d (%3d) ", k, kk);
340  }
341  writelog("\n --- Components Moles: ");
342  for (k = 0; k < nComponents; k++) {
343  kk = orderVectorSpecies[k];
344  writelogf("%-11.3g", molNumBase[kk]);
345  }
346  writelog("\n --- NonComponent | Moles | ");
347  for (i = 0; i < nComponents; i++) {
348  kk = orderVectorSpecies[i];
349  sname = mphase->speciesName(kk);
350  writelogf("%-11.10s", sname.c_str());
351  }
352  writelog("\n");
353
354  for (i = 0; i < nNonComponents; i++) {
355  k = i + nComponents;
356  kk = orderVectorSpecies[k];
357  writelogf(" --- %3d (%3d) ", k, kk);
358  sname = mphase->speciesName(kk);
359  writelogf("%-10.10s", sname.c_str());
360  writelogf("|%10.3g|", molNumBase[kk]);
361  /*
362  * Print the negative of formRxnMatrix[]; it's easier to interpret.
363  */
364  for (j = 0; j < nComponents; j++) {
365  writelogf(" %6.2f", - formRxnMatrix[j + i * ne]);
366  }
367  writelog("\n");
368  }
369  writelog(" ");
370  for (i=0; i<77; i++) {
371  writelog("-");
372  }
373  writelog("\n");
374  }
375
376  return nComponents;
377 } /* basopt() ************************************************************/
378
379 static void print_stringTrunc(const char* str, int space, int alignment)
380
381 /***********************************************************************
382  * vcs_print_stringTrunc():
383  *
384  * Print a string within a given space limit. This routine
385  * limits the amount of the string that will be printed to a
386  * maximum of "space" characters.
387  *
388  * str = String -> must be null terminated.
389  * space = space limit for the printing.
390  * alignment = 0 centered
391  * 1 right aligned
392  * 2 left aligned
393  ***********************************************************************/
394 {
395  int i, ls=0, rs=0;
396  int len = static_cast<int>(strlen(str));
397  if ((len) >= space) {
398  for (i = 0; i < space; i++) {
399  writelogf("%c", str[i]);
400  }
401  } else {
402  if (alignment == 1) {
403  ls = space - len;
404  } else if (alignment == 2) {
405  rs = space - len;
406  } else {
407  ls = (space - len) / 2;
408  rs = space - len - ls;
409  }
410  if (ls != 0) {
411  for (i = 0; i < ls; i++) {
412  writelog(" ");
413  }
414  }
415  writelogf("%s", str);
416  if (rs != 0) {
417  for (i = 0; i < rs; i++) {
418  writelog(" ");
419  }
420  }
421  }
422 }
423
424 size_t Cantera::ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
425  MultiPhase* mphase,
426  std::vector<size_t>& orderVectorSpecies,
427  std::vector<size_t>& orderVectorElements)
428 {
429  size_t j, k, l, i, jl, ml, jr, ielem, jj, kk=0;
430
431  bool lindep = false;
432  size_t nelements = mphase->nElements();
433  std::string ename;
434  /*
435  * Get the total number of species in the multiphase object
436  */
437  size_t nspecies = mphase->nSpecies();
438
439  double test = -1.0E10;
440  if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) {
441  writelog(" ");
442  for (i=0; i<77; i++) {
443  writelog("-");
444  }
445  writelog("\n");
446  writelog(" --- Subroutine ElemRearrange() called to ");
447  writelog("check stoich. coefficient matrix\n");
448  writelog(" --- and to rearrange the element ordering once\n");
449  }
450
451  /*
452  * Perhaps, initialize the element ordering
453  */
454  if (orderVectorElements.size() < nelements) {
455  orderVectorElements.resize(nelements);
456  for (j = 0; j < nelements; j++) {
457  orderVectorElements[j] = j;
458  }
459  }
460
461  /*
462  * Perhaps, initialize the species ordering. However, this is
463  * dangerous, as this ordering is assumed to yield the
464  * component species for the problem
465  */
466  if (orderVectorSpecies.size() != nspecies) {
467  orderVectorSpecies.resize(nspecies);
468  for (k = 0; k < nspecies; k++) {
469  orderVectorSpecies[k] = k;
470  }
471  }
472
473  /*
474  * If the elementAbundances aren't input, just create a fake one
475  * based on summing the column of the stoich matrix.
476  * This will force elements with zero species to the
477  * end of the element ordering.
478  */
479  vector_fp eAbund(nelements,0.0);
480  if (elementAbundances.size() != nelements) {
481  for (j = 0; j < nelements; j++) {
482  eAbund[j] = 0.0;
483  for (k = 0; k < nspecies; k++) {
484  eAbund[j] += fabs(mphase->nAtoms(k, j));
485  }
486  }
487  } else {
488  copy(elementAbundances.begin(), elementAbundances.end(),
489  eAbund.begin());
490  }
491
492  vector_fp sa(nelements,0.0);
493  vector_fp ss(nelements,0.0);
494  vector_fp sm(nelements*nelements,0.0);
495
496  /*
497  * Top of a loop of some sort based on the index JR. JR is the
498  * current number independent elements found.
499  */
500  jr = npos;
501  do {
502  ++jr;
503  /*
504  * Top of another loop point based on finding a linearly
505  * independent element
506  */
507  do {
508  /*
509  * Search the element vector. We first locate elements that
510  * are present in any amount. Then, we locate elements that
511  * are not present in any amount.
512  * Return its identity in K.
513  */
514  k = nelements;
515  for (ielem = jr; ielem < nelements; ielem++) {
516  kk = orderVectorElements[ielem];
517  if (eAbund[kk] != test && eAbund[kk] > 0.0) {
518  k = ielem;
519  break;
520  }
521  }
522  for (ielem = jr; ielem < nelements; ielem++) {
523  kk = orderVectorElements[ielem];
524  if (eAbund[kk] != test) {
525  k = ielem;
526  break;
527  }
528  }
529
530  if (k == nelements) {
531  // When we are here, there is an error usually.
532  // We haven't found the number of elements necessary.
533  if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) {
534  writelogf("Error exit: returning with nComponents = %d\n", jr);
535  }
536  throw CanteraError("ElemRearrange", "Required number of elements not found.");
537  }
538
539  /*
540  * Assign a large negative number to the element that we have
541  * just found, in order to take it out of further consideration.
542  */
543  eAbund[kk] = test;
544
545  /* *********************************************************** */
546  /* **** CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX */
547  /* **** LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX ****** */
548  /* *********************************************************** */
549  /*
550  * Modified Gram-Schmidt Method, p. 202 Dalquist
551  * QR factorization of a matrix without row pivoting.
552  */
553  jl = jr;
554  /*
555  * Fill in the row for the current element, k, under consideration
556  * The row will contain the Formula matrix value for that element
557  * with respect to the vector of component species.
558  * (note j and k indices are flipped compared to the previous routine)
559  */
560  for (j = 0; j < nComponents; ++j) {
561  jj = orderVectorSpecies[j];
562  kk = orderVectorElements[k];
563  sm[j + jr*nComponents] = mphase->nAtoms(jj,kk);
564  }
565  if (jl > 0) {
566  /*
567  * Compute the coefficients of JA column of the
568  * the upper triangular R matrix, SS(J) = R_J_JR
569  * (this is slightly different than Dalquist)
570  * R_JA_JA = 1
571  */
572  for (j = 0; j < jl; ++j) {
573  ss[j] = 0.0;
574  for (i = 0; i < nComponents; ++i) {
575  ss[j] += sm[i + jr*nComponents] * sm[i + j*nComponents];
576  }
577  ss[j] /= sa[j];
578  }
579  /*
580  * Now make the new column, (*,JR), orthogonal to the
581  * previous columns
582  */
583  for (j = 0; j < jl; ++j) {
584  for (l = 0; l < nComponents; ++l) {
585  sm[l + jr*nComponents] -= ss[j] * sm[l + j*nComponents];
586  }
587  }
588  }
589
590  /*
591  * Find the new length of the new column in Q.
592  * It will be used in the denominator in future row calcs.
593  */
594  sa[jr] = 0.0;
595  for (ml = 0; ml < nComponents; ++ml) {
596  double tmp = sm[ml + jr*nComponents];
597  sa[jr] += tmp * tmp;
598  }
599  /* **************************************************** */
600  /* **** IF NORM OF NEW ROW .LT. 1E-6 REJECT ********** */
601  /* **************************************************** */
602  if (sa[jr] < 1.0e-6) {
603  lindep = true;
604  } else {
605  lindep = false;
606  }
607  } while (lindep);
608  /* ****************************************** */
609  /* **** REARRANGE THE DATA ****************** */
610  /* ****************************************** */
611  if (jr != k) {
612  if (DEBUG_MODE_ENABLED && BasisOptimize_print_lvl > 0) {
613  kk = orderVectorElements[k];
614  ename = mphase->elementName(kk);
615  writelog(" --- ");
616  writelogf("%-2.2s", ename.c_str());
617  writelog("replaces ");
618  kk = orderVectorElements[jr];
619  ename = mphase->elementName(kk);
620  writelogf("%-2.2s", ename.c_str());
621  writelogf(" as element %3d\n", jr);
622  }
623  std::swap(orderVectorElements[jr], orderVectorElements[k]);
624  }
625
626  /*
627  * If we haven't found enough components, go back
628  * and find some more. (nc -1 is used below, because
629  * jr is counted from 0, via the C convention.
630  */
631  } while (jr < (nComponents-1));
632  return nComponents;
633 }
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
std::string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
Definition: MultiPhase.cpp:904
Headers for the MultiPhase object that is used to set up multiphase equilibrium problems (see Equilib...
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:159
A class for multiphase mixtures.
Definition: MultiPhase.h:53
void writelogf(const char *fmt,...)
Write a formatted message to the screen.
Definition: global.cpp:38
std::string elementName(size_t m) const
Returns the name of the global element m.
Definition: MultiPhase.cpp:875
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
void getMoles(doublereal *molNum) const
Get the mole numbers of all species in the multiphase object.
Definition: MultiPhase.cpp:441
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:909
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
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
static void print_stringTrunc(const char *str, int space, int alignment)
Print a string within a given space limit.
void writelog(const std::string &msg)
Write a message to the screen.
Definition: global.cpp:33
size_t nSpecies() const
Number of species, summed over all phases.
Definition: MultiPhase.h:124
size_t 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...
size_t nElements() const
Number of elements.
Definition: MultiPhase.h:98