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