Cantera  2.3.0
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 // This file is part of Cantera. See License.txt in the top-level directory or
9
14
15 namespace Cantera
16 {
17
18 int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
19  double* const sm, double* const ss)
20 {
21  size_t ncomponents = m_numComponents;
22  if (m_debug_print_lvl >= 2) {
23  plogf(" ");
24  for (size_t i=0; i<77; i++) {
25  plogf("-");
26  }
27  plogf("\n");
28  plogf(" --- Subroutine elem_rearrange() called to ");
29  plogf("check stoich. coefficient matrix\n");
30  plogf(" --- and to rearrange the element ordering once");
31  plogendl();
32  }
33
34  // Use a temporary work array for the element numbers
35  // Also make sure the value of test is unique.
36  bool lindep = true;
37  double test = -1.0E10;
38  while (lindep) {
39  lindep = false;
40  for (size_t i = 0; i < m_numElemConstraints; ++i) {
41  test -= 1.0;
42  aw[i] = m_elemAbundancesGoal[i];
43  if (test == aw[i]) {
44  lindep = true;
45  }
46  }
47  }
48
49  // Top of a loop of some sort based on the index JR. JR is the current
50  // number independent elements found.
51  size_t jr = 0;
52  while (jr < ncomponents) {
53  size_t k;
54
55  // Top of another loop point based on finding a linearly independent
56  // species
57  while (true) {
58  // Search the remaining part of the mole fraction vector, AW, for
59  // the largest remaining species. Return its identity in K.
61  for (size_t ielem = jr; ielem < m_numElemConstraints; ielem++) {
62  if (m_elementActive[ielem] && aw[ielem] != test) {
63  k = ielem;
64  break;
65  }
66  }
67  if (k == m_numElemConstraints) {
68  throw CanteraError("vcs_elem_rearrange",
69  "Shouldn't be here. Algorithm misfired.");
70  }
71
72  // Assign a large negative number to the element that we have just
73  // found, in order to take it out of further consideration.
74  aw[k] = test;
75
76  // CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX LINE WITH
77  // PREVIOUS LINES OF THE FORMULA MATRIX
78  //
79  // Modified Gram-Schmidt Method, p. 202 Dalquist QR factorization of
80  // a matrix without row pivoting.
81  size_t jl = jr;
82
83  // Fill in the row for the current element, k, under consideration
84  // The row will contain the Formula matrix value for that element
85  // from the current component.
86  for (size_t j = 0; j < ncomponents; ++j) {
87  sm[j + jr*ncomponents] = m_formulaMatrix(j,k);
88  }
89  if (jl > 0) {
90  // Compute the coefficients of JA column of the the upper
91  // triangular R matrix, SS(J) = R_J_JR (this is slightly
92  // different than Dalquist) R_JA_JA = 1
93  for (size_t j = 0; j < jl; ++j) {
94  ss[j] = 0.0;
95  for (size_t i = 0; i < ncomponents; ++i) {
96  ss[j] += sm[i + jr*ncomponents] * sm[i + j*ncomponents];
97  }
98  ss[j] /= sa[j];
99  }
100
101  // Now make the new column, (*,JR), orthogonal to the previous
102  // columns
103  for (size_t j = 0; j < jl; ++j) {
104  for (size_t i = 0; i < ncomponents; ++i) {
105  sm[i + jr*ncomponents] -= ss[j] * sm[i + j*ncomponents];
106  }
107  }
108  }
109
110  // Find the new length of the new column in Q. It will be used in
111  // the denominator in future row calcs.
112  sa[jr] = 0.0;
113  for (size_t ml = 0; ml < ncomponents; ++ml) {
114  sa[jr] += pow(sm[ml + jr*ncomponents], 2);
115  }
116  // IF NORM OF NEW ROW .LT. 1E-6 REJECT
117  if (sa[jr] > 1.0e-6) {
118  break;
119  }
120  }
121  // REARRANGE THE DATA
122  if (jr != k) {
123  if (m_debug_print_lvl >= 2) {
124  plogf(" --- ");
125  plogf("%-2.2s", m_elementName[k]);
126  plogf("(%9.2g) replaces ", m_elemAbundancesGoal[k]);
127  plogf("%-2.2s", m_elementName[jr]);
128  plogf("(%9.2g) as element %3d", m_elemAbundancesGoal[jr], jr);
129  plogendl();
130  }
131  vcs_switch_elem_pos(jr, k);
132  std::swap(aw[jr], aw[k]);
133  }
134
135  // If we haven't found enough components, go back and find some more.
136  jr++;
137  }
138  return VCS_SUCCESS;
139 }
140
141 void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos)
142 {
143  if (ipos == jpos) {
144  return;
145  }
147  "vcs_switch_elem_pos",
148  "inappropriate args: {} {}", ipos, jpos);
149
150  // Change the element Global Index list in each vcs_VolPhase object
151  // to reflect the switch in the element positions.
152  for (size_t iph = 0; iph < m_numPhases; iph++) {
153  vcs_VolPhase* volPhase = m_VolPhaseList[iph];
154  for (size_t e = 0; e < volPhase->nElemConstraints(); e++) {
155  if (volPhase->elemGlobalIndex(e) == ipos) {
156  volPhase->setElemGlobalIndex(e, jpos);
157  }
158  if (volPhase->elemGlobalIndex(e) == jpos) {
159  volPhase->setElemGlobalIndex(e, ipos);
160  }
161  }
162  }
163  std::swap(m_elemAbundancesGoal[ipos], m_elemAbundancesGoal[jpos]);
164  std::swap(m_elemAbundances[ipos], m_elemAbundances[jpos]);
165  std::swap(m_elementMapIndex[ipos], m_elementMapIndex[jpos]);
166  std::swap(m_elType[ipos], m_elType[jpos]);
167  std::swap(m_elementActive[ipos], m_elementActive[jpos]);
168  for (size_t j = 0; j < m_numSpeciesTot; ++j) {
169  std::swap(m_formulaMatrix(j,ipos), m_formulaMatrix(j,jpos));
170  }
171  std::swap(m_elementName[ipos], m_elementName[jpos]);
172 }
173 }
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...
vector_fp m_elemAbundances
Element abundances vector.
Definition: vcs_solve.h:1568
size_t nElemConstraints() const
Returns the number of element constraints.
std::vector< vcs_VolPhase * > m_VolPhaseList
Array of Phase Structures. Length = number of phases.
Definition: vcs_solve.h:1740
std::vector< size_t > m_elementMapIndex
Index vector that keeps track of the rearrangement of the elements.
Definition: vcs_solve.h:1679
size_t elemGlobalIndex(const size_t e) const
Returns the global index of the local element index for the phase.
vector_int m_elementActive
Specifies whether an element constraint is active.
Definition: vcs_solve.h:1737
std::vector< std::string > m_elementName
Vector of strings containing the element names.
Definition: vcs_solve.h:1716
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:18
size_t m_numElemConstraints
Number of element constraints in the problem.
Definition: vcs_solve.h:1393
size_t m_numSpeciesTot
Total number of species in the problems.
Definition: vcs_solve.h:1387
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
int m_debug_print_lvl
Debug printing lvl.
Definition: vcs_solve.h:1863
size_t m_numComponents
Number of components calculated for the problem.
Definition: vcs_solve.h:1396
vector_fp m_elemAbundancesGoal
Element abundances vector Goals.
Definition: vcs_solve.h:1577
Phase information and Phase calculations for vcs.
Definition: vcs_VolPhase.h:81
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:270
Array2D m_formulaMatrix
Formula matrix for the problem.
Definition: vcs_solve.h:1422
#define plogendl()
define this Cantera function to replace cout << endl;
Definition: vcs_internal.h:25
Contains declarations for string manipulation functions within Cantera.
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:18
vector_int m_elType
Type of the element constraint.
Definition: vcs_solve.h:1731
Namespace for the Cantera kernel.
Definition: application.cpp:29
size_t m_numPhases
Number of Phases in the problem.
Definition: vcs_solve.h:1414
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.
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