Cantera  2.2.1
vcs_elem_rearrange.cpp
Go to the documentation of this file.
1 /**
2  * @file vcs_elem_rearrange.cpp
3  * Contains implementations for rearranging the element columns, and
4  * it contains the algorithm for choosing the rearrangement.
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
16
17 namespace Cantera
18 {
19
20 int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
21  double* const sm, double* const ss)
22 {
23  size_t ncomponents = m_numComponents;
24  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
25  plogf(" ");
26  for (size_t i=0; i<77; i++) {
27  plogf("-");
28  }
29  plogf("\n");
30  plogf(" --- Subroutine elem_rearrange() called to ");
31  plogf("check stoich. coefficient matrix\n");
32  plogf(" --- and to rearrange the element ordering once");
33  plogendl();
34  }
35
36  /*
37  * Use a temporary work array for the element numbers
38  * Also make sure the value of test is unique.
39  */
40  bool lindep = false;
41  double test = -1.0E10;
42  do {
43  lindep = false;
44  for (size_t i = 0; i < m_numElemConstraints; ++i) {
45  test -= 1.0;
46  aw[i] = m_elemAbundancesGoal[i];
47  if (test == aw[i]) {
48  lindep = true;
49  }
50  }
51  } while (lindep);
52
53  /*
54  * Top of a loop of some sort based on the index JR. JR is the
55  * current number independent elements found.
56  */
57  size_t jr = npos;
58  do {
59  ++jr;
60  size_t k;
61  /*
62  * Top of another loop point based on finding a linearly
63  * independent species
64  */
65  do {
66  /*
67  * Search the remaining part of the mole fraction vector, AW,
68  * for the largest remaining species. Return its identity in K.
69  */
71  for (size_t ielem = jr; ielem < m_numElemConstraints; ielem++) {
72  if (m_elementActive[ielem]) {
73  if (aw[ielem] != test) {
74  k = ielem;
75  break;
76  }
77  }
78  }
79  if (k == m_numElemConstraints) {
80  throw CanteraError("vcs_elem_rearrange",
81  "Shouldn't be here. Algorithm misfired.");
82  }
83
84  /*
85  * Assign a large negative number to the element that we have
86  * just found, in order to take it out of further consideration.
87  */
88  aw[k] = test;
89
90  /* *********************************************************** */
91  /* **** CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX */
92  /* **** LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX ****** */
93  /* *********************************************************** */
94  /*
95  * Modified Gram-Schmidt Method, p. 202 Dalquist
96  * QR factorization of a matrix without row pivoting.
97  */
98  size_t jl = jr;
99  /*
100  * Fill in the row for the current element, k, under consideration
101  * The row will contain the Formula matrix value for that element
102  * from the current component.
103  */
104  for (size_t j = 0; j < ncomponents; ++j) {
105  sm[j + jr*ncomponents] = m_formulaMatrix(j,k);
106  }
107  if (jl > 0) {
108  /*
109  * Compute the coefficients of JA column of the
110  * the upper triangular R matrix, SS(J) = R_J_JR
111  * (this is slightly different than Dalquist)
112  * R_JA_JA = 1
113  */
114  for (size_t j = 0; j < jl; ++j) {
115  ss[j] = 0.0;
116  for (size_t i = 0; i < ncomponents; ++i) {
117  ss[j] += sm[i + jr*ncomponents] * sm[i + j*ncomponents];
118  }
119  ss[j] /= sa[j];
120  }
121  /*
122  * Now make the new column, (*,JR), orthogonal to the
123  * previous columns
124  */
125  for (size_t j = 0; j < jl; ++j) {
126  for (size_t l = 0; l < ncomponents; ++l) {
127  sm[l + jr*ncomponents] -= ss[j] * sm[l + j*ncomponents];
128  }
129  }
130  }
131
132  /*
133  * Find the new length of the new column in Q.
134  * It will be used in the denominator in future row calcs.
135  */
136  sa[jr] = 0.0;
137  for (size_t ml = 0; ml < ncomponents; ++ml) {
138  sa[jr] += pow(sm[ml + jr*ncomponents], 2);
139  }
140  /* **************************************************** */
141  /* **** IF NORM OF NEW ROW .LT. 1E-6 REJECT ********** */
142  /* **************************************************** */
143  if (sa[jr] < 1.0e-6) {
144  lindep = true;
145  } else {
146  lindep = false;
147  }
148  } while (lindep);
149  /* ****************************************** */
150  /* **** REARRANGE THE DATA ****************** */
151  /* ****************************************** */
152  if (jr != k) {
153  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
154  plogf(" --- ");
155  plogf("%-2.2s", (m_elementName[k]).c_str());
156  plogf("(%9.2g) replaces ", m_elemAbundancesGoal[k]);
157  plogf("%-2.2s", (m_elementName[jr]).c_str());
158  plogf("(%9.2g) as element %3d", m_elemAbundancesGoal[jr], jr);
159  plogendl();
160  }
161  vcs_switch_elem_pos(jr, k);
162  std::swap(aw[jr], aw[k]);
163  }
164
165  /*
166  * If we haven't found enough components, go back
167  * and find some more. (nc -1 is used below, because
168  * jr is counted from 0, via the C convention.
169  */
170  } while (jr < (ncomponents-1));
171  return VCS_SUCCESS;
172 }
173
174 void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos)
175 {
176  if (ipos == jpos) {
177  return;
178  }
180  "vcs_switch_elem_pos",
181  "inappropriate args: " + int2str(ipos) + " " + int2str(jpos));
182  /*
183  * Change the element Global Index list in each vcs_VolPhase object
184  * to reflect the switch in the element positions.
185  */
186  for (size_t iph = 0; iph < m_numPhases; iph++) {
187  vcs_VolPhase* volPhase = m_VolPhaseList[iph];
188  for (size_t e = 0; e < volPhase->nElemConstraints(); e++) {
189  if (volPhase->elemGlobalIndex(e) == ipos) {
190  volPhase->setElemGlobalIndex(e, jpos);
191  }
192  if (volPhase->elemGlobalIndex(e) == jpos) {
193  volPhase->setElemGlobalIndex(e, ipos);
194  }
195  }
196  }
197  std::swap(m_elemAbundancesGoal[ipos], m_elemAbundancesGoal[jpos]);
198  std::swap(m_elemAbundances[ipos], m_elemAbundances[jpos]);
199  std::swap(m_elementMapIndex[ipos], m_elementMapIndex[jpos]);
200  std::swap(m_elType[ipos], m_elType[jpos]);
201  std::swap(m_elementActive[ipos], m_elementActive[jpos]);
202  for (size_t j = 0; j < m_numSpeciesTot; ++j) {
203  std::swap(m_formulaMatrix(j,ipos), m_formulaMatrix(j,jpos));
204  }
205  std::swap(m_elementName[ipos], m_elementName[jpos]);
206 }
207 }
size_t nElemConstraints() const
Returns the number of element constraints.
int vcs_elem_rearrange(double *const aw, double *const sa, double *const sm, double *const ss)
Rearrange the constraint equations represented by the Formula Matrix so that the operational ones are...
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:39
std::vector< int > m_elType
Type of the element constraint.
Definition: vcs_solve.h:1860
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
std::vector< vcs_VolPhase * > m_VolPhaseList
Array of Phase Structures. Length = number of phases.
Definition: vcs_solve.h:1869
std::vector< size_t > m_elementMapIndex
Index vector that keeps track of the rearrangement of the elements.
Definition: vcs_solve.h:1808
#define AssertThrowMsg(expr, procedure, message)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:301
std::vector< std::string > m_elementName
Vector of strings containing the element names.
Definition: vcs_solve.h:1845
Header file for the internal object that holds the vcs equilibrium problem (see Class VCS_SOLVE and E...
Header for the object representing each phase within vcs.
#define VCS_SUCCESS
Definition: vcs_defs.h:21
std::vector< double > m_elemAbundances
Element abundances vector.
Definition: vcs_solve.h:1686
size_t m_numElemConstraints
Number of element constraints in the problem.
Definition: vcs_solve.h:1499
size_t m_numSpeciesTot
Total number of species in the problems.
Definition: vcs_solve.h:1493
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
int m_debug_print_lvl
Debug printing lvl.
Definition: vcs_solve.h:1995
size_t m_numComponents
Number of components calculated for the problem.
Definition: vcs_solve.h:1502
Phase information and Phase calculations for vcs.
Definition: vcs_VolPhase.h:86
Array2D m_formulaMatrix
Formula matrix for the problem.
Definition: vcs_solve.h:1534
std::vector< int > m_elementActive
Specifies whether an element constraint is active.
Definition: vcs_solve.h:1866
#define plogendl()
define this Cantera function to replace cout << endl;
Definition: vcs_internal.h:31
Contains declarations for string manipulation functions within Cantera.
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:24
size_t elemGlobalIndex(const size_t e) const
Returns the global index of the local element index for the phase.
size_t m_numPhases
Number of Phases in the problem.
Definition: vcs_solve.h:1526
void vcs_switch_elem_pos(size_t ipos, size_t jpos)
Swaps the indices for all of the global data for two elements, ipos and jpos.
std::vector< double > m_elemAbundancesGoal
Element abundances vector Goals.
Definition: vcs_solve.h:1695
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
void setElemGlobalIndex(const size_t eLocal, const size_t eGlobal)
sets a local phase element to a global index value