Cantera  2.0
vcs_rank.cpp
Go to the documentation of this file.
1
2 /*!
3  * @file vcs_rank.cpp
4  * Header file for the internal class that holds the problem.
5  */
6 /*
7  * $Id: vcs_solve.cpp 735 2011-07-25 14:44:41Z hkmoffa$
8  */
9 /*
10  * Copywrite (2005) Sandia Corporation. Under the terms of
11  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
12  * U.S. Government retains certain rights in this software.
13  */
14
15
18 #include "cantera/equil/vcs_prob.h"
19
21 #include "cantera/equil/vcs_SpeciesProperties.h"
22 #include "cantera/equil/vcs_species_thermo.h"
23
24 #include "cantera/base/clockWC.h"
26
27 #include <string>
28 #include <cstdio>
29 #include "math.h"
30 using namespace std;
31
32 namespace VCSnonideal {
33  //====================================================================================================================
34  static int basisOptMax1(const double * const molNum,
35  const int n) {
36  // int largest = 0;
37
38  for (int i = 0; i < n; ++i) {
39
40  if (molNum[i] > -1.0E200 && fabs(molNum[i]) > 1.0E-13) {
41  return i;
42  }
43  }
44  for (int i = 0; i < n; ++i) {
45
46  if (molNum[i] > -1.0E200) {
47  return i;
48  }
49  }
50  return n-1;
51  }
52  //====================================================================================================================
53  // Calculate the rank of a matrix and return the rows and columns that will generate an independent basis
54  // for that rank
55  /*
56  * Choose the optimum component species basis for the calculations, finding the rank and
57  * set of linearly independent rows for that calculation.
58  * Then find the set of linearly indepedent element columns that can support that rank.
59  * This is done by taking the transpose of the matrix and redoing the same calculation.
60  * (there may be a better way to do this. I don't know.)
61  *
62  *
63  * Input
64  * ---------
65  *
66  * @param awtmp Vector of mole numbers which will be used to construct a
67  * ranking for how to pick the basis species. This is largely ignored
68  * here.
69  *
70  * @param numSpecies Number of species. This is the number of rows in the matrix.
71  *
72  * @param matrix Matrix. This is the formula matrix. Nominally, the rows are species, while
73  * the columns are element compositions. However, this routine
74  * is totally general, so that the rows and columns can be anything.
75  *
76  * @param numElemConstraints Number of element constraints
77  *
78  * Output
79  * ---------
80  * @param usedZeroedSpecies = If true, then a species with a zero concentration
81  * was used as a component.
82  *
83  *
84  * @param compRes Vector of rows which are linearly independent. (these are the components)
85  *
86  * @param elemComp Vector of columns which are linearly independent (These are the actionable element
87  * constraints).
88  *
89  * @return Returns number of components. This is the rank of the matrix
90  */
91  int VCS_SOLVE::vcs_rank(const double * awtmp, size_t numSpecies, const double matrix[], size_t numElemConstraints,
92  std::vector<size_t> &compRes, std::vector<size_t>& elemComp, int * const usedZeroedSpecies) const
93  {
94
95  int lindep;
96  size_t j, k, jl, i, l, ml;
97  int numComponents = 0;
98
99  compRes.clear();
100  elemComp.clear();
101  vector<double> sm(numElemConstraints*numSpecies);
102  vector<double> sa(numSpecies);
103  vector<double> ss(numSpecies);
104
105  double test = -0.2512345E298;
106 #ifdef DEBUG_MODE
107  if (m_debug_print_lvl >= 2) {
108  plogf(" "); for(i=0; i<77; i++) plogf("-"); plogf("\n");
109  plogf(" --- Subroutine vcs_rank called to ");
110  plogf("calculate the rank and independent rows /colums of the following matrix\n");
111  if (m_debug_print_lvl >= 5) {
112  plogf(" --- Species | ");
113  for (j = 0; j < numElemConstraints; j++) {
114  plogf(" ");
115  plogf(" %3d ", j);
116  }
117  plogf("\n");
118  plogf(" --- -----------");
119  for (j = 0; j < numElemConstraints; j++) {
120  plogf("---------");
121  }
122  plogf("\n");
123  for (k = 0; k < numSpecies; k++) {
124  plogf(" --- ");
125  plogf(" %3d ", k);
126  plogf(" |");
127  for (j = 0; j < numElemConstraints; j++) {
128  plogf(" %8.2g", matrix[j*numSpecies + k]);
129  }
130  plogf("\n");
131  }
132  plogf(" ---");
133  plogendl();
134  }
135  }
136 #endif
137
138  /*
139  * Calculate the maximum value of the number of components possible
140  * It's equal to the minimum of the number of elements and the
141  * number of total species.
142  */
143  int ncTrial = std::min(numElemConstraints, numSpecies);
144  numComponents = ncTrial;
145  *usedZeroedSpecies = false;
146
147  /*
148  * Use a temporary work array for the mole numbers, aw[]
149  */
150  std::vector<double> aw(numSpecies);
151  for (j = 0; j < numSpecies; j++) {
152  aw[j] = awtmp[j];
153  }
154
155  int jr = -1;
156  /*
157  * Top of a loop of some sort based on the index JR. JR is the
158  * current number of component species found.
159  */
160  do {
161  ++jr;
162  /* - Top of another loop point based on finding a linearly */
163  /* - independent species */
164  do {
165  /*
166  * Search the remaining part of the mole number vector, AW,
167  * for the largest remaining species. Return its identity in K.
168  * The first search criteria is always the largest positive
169  * magnitude of the mole number.
170  */
171  k = basisOptMax1(VCS_DATA_PTR(aw), numSpecies);
172
173  if ((aw[k] != test) && fabs(aw[k]) == 0.0) {
174  *usedZeroedSpecies = true;
175  }
176
177
178  if (aw[k] == test) {
179  numComponents = jr;
180
181  goto L_CLEANUP;
182  }
183  /*
184  * Assign a small negative number to the component that we have
185  * just found, in order to take it out of further consideration.
186  */
187  aw[k] = test;
188  /* *********************************************************** */
189  /* **** CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES ****** */
190  /* *********************************************************** */
191  /*
192  * Modified Gram-Schmidt Method, p. 202 Dalquist
193  * QR factorization of a matrix without row pivoting.
194  */
195  jl = jr;
196  for (j = 0; j < numElemConstraints; ++j) {
197  sm[j + jr*numElemConstraints] = matrix[j*numSpecies + k];
198  }
199  if (jl > 0) {
200  /*
201  * Compute the coefficients of JA column of the
202  * the upper triangular R matrix, SS(J) = R_J_JR
203  * (this is slightly different than Dalquist)
204  * R_JA_JA = 1
205  */
206  for (j = 0; j < jl; ++j) {
207  ss[j] = 0.0;
208  for (i = 0; i < numElemConstraints; ++i) {
209  ss[j] += sm[i + jr* numElemConstraints] * sm[i + j* numElemConstraints];
210  }
211  ss[j] /= sa[j];
212  }
213  /*
214  * Now make the new column, (*,JR), orthogonal to the
215  * previous columns
216  */
217  for (j = 0; j < jl; ++j) {
218  for (l = 0; l < numElemConstraints; ++l) {
219  sm[l + jr*numElemConstraints] -= ss[j] * sm[l + j*numElemConstraints];
220  }
221  }
222  }
223  /*
224  * Find the new length of the new column in Q.
225  * It will be used in the denominator in future row calcs.
226  */
227  sa[jr] = 0.0;
228  for (ml = 0; ml < numElemConstraints; ++ml) {
229  sa[jr] += SQUARE(sm[ml + jr * numElemConstraints]);
230  }
231  /* **************************************************** */
232  /* **** IF NORM OF NEW ROW .LT. 1E-3 REJECT ********** */
233  /* **************************************************** */
234  if (sa[jr] < 1.0e-6) lindep = true;
235  else lindep = false;
236  } while(lindep);
237  /* ****************************************** */
238  /* **** REARRANGE THE DATA ****************** */
239  /* ****************************************** */
240  compRes.push_back(k);
241  elemComp.push_back(jr);
242
243  } while (jr < (ncTrial-1));
244
245  L_CLEANUP: ;
246
247  if (numComponents == ncTrial && numElemConstraints == numSpecies) {
248  return numComponents;
249  }
250
251
252  int numComponentsR = numComponents;
253
254  ss.resize(numElemConstraints);
255  sa.resize(numElemConstraints);
256
257
258  elemComp.clear();
259
260  aw.resize(numElemConstraints);
261  for (j = 0; j < numSpecies; j++) {
262  aw[j] = 1.0;
263  }
264
265  jr = -1;
266
267  do {
268  ++jr;
269
270  do {
271
272  k = basisOptMax1(VCS_DATA_PTR(aw), numElemConstraints);
273
274  if (aw[k] == test) {
275  numComponents = jr;
276  goto LE_CLEANUP;
277  }
278  aw[k] = test;
279
280
281  jl = jr;
282  for (j = 0; j < numSpecies; ++j) {
283  sm[j + jr*numSpecies] = matrix[k*numSpecies + j];
284  }
285  if (jl > 0) {
286
287  for (j = 0; j < jl; ++j) {
288  ss[j] = 0.0;
289  for (i = 0; i < numSpecies; ++i) {
290  ss[j] += sm[i + jr* numSpecies] * sm[i + j* numSpecies];
291  }
292  ss[j] /= sa[j];
293  }
294
295  for (j = 0; j < jl; ++j) {
296  for (l = 0; l < numSpecies; ++l) {
297  sm[l + jr*numSpecies] -= ss[j] * sm[l + j*numSpecies];
298  }
299  }
300  }
301
302  sa[jr] = 0.0;
303  for (ml = 0; ml < numSpecies; ++ml) {
304  sa[jr] += SQUARE(sm[ml + jr * numSpecies]);
305  }
306
307  if (sa[jr] < 1.0e-6) lindep = true;
308  else lindep = false;
309  } while(lindep);
310
311  elemComp.push_back(k);
312
313  } while (jr < (ncTrial-1));
314  numComponents = jr;
315  LE_CLEANUP: ;
316
317 #ifdef DEBUG_MODE
318  if (m_debug_print_lvl >= 2) {
319  plogf(" --- vcs_rank found rank %d\n", numComponents);
320  if (m_debug_print_lvl >= 5) {
321  if (compRes.size() == elemComp.size()) {
322  printf(" --- compRes elemComp\n");
323  for (int i = 0; i < (int) compRes.size(); i++) {
324  printf(" --- %d %d \n", (int) compRes[i], (int) elemComp[i]);
325  }
326  } else {
327  for (int i = 0; i < (int) compRes.size(); i++) {
328  printf(" --- compRes[%d] = %d \n", (int) i, (int) compRes[i]);
329  }
330  for (int i = 0; i < (int) elemComp.size(); i++) {
331  printf(" --- elemComp[%d] = %d \n", (int) i, (int) elemComp[i]);
332  }
333  }
334  }
335  }
336 #endif
337
338  if (numComponentsR != numComponents) {
339  printf("vcs_rank ERROR: number of components are different: %d %d\n", numComponentsR, numComponents);
340  throw Cantera::CanteraError("vcs_rank ERROR:",
341  " logical inconsistency");
342  exit(-1);
343  }
344  return numComponents;
345  }
346
347
348 }