Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
vcs_rank.cpp
Go to the documentation of this file.
1 /*!
2  * @file vcs_rank.cpp
3  * Header file for the internal class that holds the problem.
4  */
5 
6 /*
7  * Copyright (2005) Sandia Corporation. Under the terms of
8  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
9  * U.S. Government retains certain rights in this software.
10  */
11 
14 
15 #include <cstdio>
16 using namespace std;
17 
18 namespace Cantera {
19  static int basisOptMax1(const double * const molNum,
20  const int n) {
21  // int largest = 0;
22 
23  for (int i = 0; i < n; ++i) {
24 
25  if (molNum[i] > -1.0E200 && fabs(molNum[i]) > 1.0E-13) {
26  return i;
27  }
28  }
29  for (int i = 0; i < n; ++i) {
30 
31  if (molNum[i] > -1.0E200) {
32  return i;
33  }
34  }
35  return n-1;
36  }
37 
38  int VCS_SOLVE::vcs_rank(const double * awtmp, size_t numSpecies, const double matrix[], size_t numElemConstraints,
39  std::vector<size_t> &compRes, std::vector<size_t>& elemComp, int * const usedZeroedSpecies) const
40  {
41  warn_deprecated("VCS_SOLVE::vcs_rank", "To be removed after Cantera 2.2");
42  int lindep;
43  size_t j, k, jl, i, l, ml;
44  int numComponents = 0;
45 
46  compRes.clear();
47  elemComp.clear();
48  vector<double> sm(numElemConstraints*numSpecies);
49  vector<double> sa(numSpecies);
50  vector<double> ss(numSpecies);
51 
52  double test = -0.2512345E298;
53  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
54  plogf(" "); for(i=0; i<77; i++) plogf("-"); plogf("\n");
55  plogf(" --- Subroutine vcs_rank called to ");
56  plogf("calculate the rank and independent rows /colums of the following matrix\n");
57  if (m_debug_print_lvl >= 5) {
58  plogf(" --- Species | ");
59  for (j = 0; j < numElemConstraints; j++) {
60  plogf(" ");
61  plogf(" %3d ", j);
62  }
63  plogf("\n");
64  plogf(" --- -----------");
65  for (j = 0; j < numElemConstraints; j++) {
66  plogf("---------");
67  }
68  plogf("\n");
69  for (k = 0; k < numSpecies; k++) {
70  plogf(" --- ");
71  plogf(" %3d ", k);
72  plogf(" |");
73  for (j = 0; j < numElemConstraints; j++) {
74  plogf(" %8.2g", matrix[j*numSpecies + k]);
75  }
76  plogf("\n");
77  }
78  plogf(" ---");
79  plogendl();
80  }
81  }
82  /*
83  * Calculate the maximum value of the number of components possible
84  * It's equal to the minimum of the number of elements and the
85  * number of total species.
86  */
87  int ncTrial = static_cast<int>(std::min(numElemConstraints, numSpecies));
88  numComponents = ncTrial;
89  *usedZeroedSpecies = false;
90 
91  /*
92  * Use a temporary work array for the mole numbers, aw[]
93  */
94  std::vector<double> aw(numSpecies);
95  for (j = 0; j < numSpecies; j++) {
96  aw[j] = awtmp[j];
97  }
98 
99  int jr = -1;
100  /*
101  * Top of a loop of some sort based on the index JR. JR is the
102  * current number of component species found.
103  */
104  do {
105  ++jr;
106  /* - Top of another loop point based on finding a linearly */
107  /* - independent species */
108  do {
109  /*
110  * Search the remaining part of the mole number vector, AW,
111  * for the largest remaining species. Return its identity in K.
112  * The first search criteria is always the largest positive
113  * magnitude of the mole number.
114  */
115  k = basisOptMax1(VCS_DATA_PTR(aw), static_cast<int>(numSpecies));
116 
117  if ((aw[k] != test) && fabs(aw[k]) == 0.0) {
118  *usedZeroedSpecies = true;
119  }
120 
121 
122  if (aw[k] == test) {
123  numComponents = jr;
124 
125  goto L_CLEANUP;
126  }
127  /*
128  * Assign a small negative number to the component that we have
129  * just found, in order to take it out of further consideration.
130  */
131  aw[k] = test;
132  /* *********************************************************** */
133  /* **** CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES ****** */
134  /* *********************************************************** */
135  /*
136  * Modified Gram-Schmidt Method, p. 202 Dalquist
137  * QR factorization of a matrix without row pivoting.
138  */
139  jl = jr;
140  for (j = 0; j < numElemConstraints; ++j) {
141  sm[j + jr*numElemConstraints] = matrix[j*numSpecies + k];
142  }
143  if (jl > 0) {
144  /*
145  * Compute the coefficients of JA column of the
146  * the upper triangular R matrix, SS(J) = R_J_JR
147  * (this is slightly different than Dalquist)
148  * R_JA_JA = 1
149  */
150  for (j = 0; j < jl; ++j) {
151  ss[j] = 0.0;
152  for (i = 0; i < numElemConstraints; ++i) {
153  ss[j] += sm[i + jr* numElemConstraints] * sm[i + j* numElemConstraints];
154  }
155  ss[j] /= sa[j];
156  }
157  /*
158  * Now make the new column, (*,JR), orthogonal to the
159  * previous columns
160  */
161  for (j = 0; j < jl; ++j) {
162  for (l = 0; l < numElemConstraints; ++l) {
163  sm[l + jr*numElemConstraints] -= ss[j] * sm[l + j*numElemConstraints];
164  }
165  }
166  }
167  /*
168  * Find the new length of the new column in Q.
169  * It will be used in the denominator in future row calcs.
170  */
171  sa[jr] = 0.0;
172  for (ml = 0; ml < numElemConstraints; ++ml) {
173  sa[jr] += pow(sm[ml + jr * numElemConstraints], 2);
174  }
175  /* **************************************************** */
176  /* **** IF NORM OF NEW ROW .LT. 1E-3 REJECT ********** */
177  /* **************************************************** */
178  if (sa[jr] < 1.0e-6) lindep = true;
179  else lindep = false;
180  } while(lindep);
181  /* ****************************************** */
182  /* **** REARRANGE THE DATA ****************** */
183  /* ****************************************** */
184  compRes.push_back(k);
185  elemComp.push_back(jr);
186 
187  } while (jr < (ncTrial-1));
188 
189  L_CLEANUP: ;
190 
191  if (numComponents == ncTrial && numElemConstraints == numSpecies) {
192  return numComponents;
193  }
194 
195 
196  int numComponentsR = numComponents;
197 
198  ss.resize(numElemConstraints);
199  sa.resize(numElemConstraints);
200 
201 
202  elemComp.clear();
203 
204  aw.resize(numElemConstraints);
205  for (j = 0; j < numSpecies; j++) {
206  aw[j] = 1.0;
207  }
208 
209  jr = -1;
210 
211  do {
212  ++jr;
213 
214  do {
215 
216  k = basisOptMax1(VCS_DATA_PTR(aw), static_cast<int>(numElemConstraints));
217 
218  if (aw[k] == test) {
219  numComponents = jr;
220  goto LE_CLEANUP;
221  }
222  aw[k] = test;
223 
224 
225  jl = jr;
226  for (j = 0; j < numSpecies; ++j) {
227  sm[j + jr*numSpecies] = matrix[k*numSpecies + j];
228  }
229  if (jl > 0) {
230 
231  for (j = 0; j < jl; ++j) {
232  ss[j] = 0.0;
233  for (i = 0; i < numSpecies; ++i) {
234  ss[j] += sm[i + jr* numSpecies] * sm[i + j* numSpecies];
235  }
236  ss[j] /= sa[j];
237  }
238 
239  for (j = 0; j < jl; ++j) {
240  for (l = 0; l < numSpecies; ++l) {
241  sm[l + jr*numSpecies] -= ss[j] * sm[l + j*numSpecies];
242  }
243  }
244  }
245 
246  sa[jr] = 0.0;
247  for (ml = 0; ml < numSpecies; ++ml) {
248  sa[jr] += pow(sm[ml + jr * numSpecies], 2);
249  }
250 
251  if (sa[jr] < 1.0e-6) lindep = true;
252  else lindep = false;
253  } while(lindep);
254 
255  elemComp.push_back(k);
256 
257  } while (jr < (ncTrial-1));
258  numComponents = jr;
259  LE_CLEANUP: ;
260 
261  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
262  plogf(" --- vcs_rank found rank %d\n", numComponents);
263  if (m_debug_print_lvl >= 5) {
264  if (compRes.size() == elemComp.size()) {
265  printf(" --- compRes elemComp\n");
266  for (int i = 0; i < (int) compRes.size(); i++) {
267  printf(" --- %d %d \n", (int) compRes[i], (int) elemComp[i]);
268  }
269  } else {
270  for (int i = 0; i < (int) compRes.size(); i++) {
271  printf(" --- compRes[%d] = %d \n", (int) i, (int) compRes[i]);
272  }
273  for (int i = 0; i < (int) elemComp.size(); i++) {
274  printf(" --- elemComp[%d] = %d \n", (int) i, (int) elemComp[i]);
275  }
276  }
277  }
278  }
279 
280  if (numComponentsR != numComponents) {
281  printf("vcs_rank ERROR: number of components are different: %d %d\n", numComponentsR, numComponents);
282  throw CanteraError("vcs_rank ERROR:",
283  " logical inconsistency");
284  exit(-1);
285  }
286  return numComponents;
287  }
288 
289 }
#define VCS_DATA_PTR(vvv)
Points to the data in a std::vector<> object.
Definition: vcs_internal.h:18
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:78
Header file for the internal object that holds the vcs equilibrium problem (see Class VCS_SOLVE and E...
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
#define plogendl()
define this Cantera function to replace cout << endl;
Definition: vcs_internal.h:31
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:24
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...