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