Cantera  2.5.1
vcs_setMolesLinProg.cpp
Go to the documentation of this file.
1 /*!
2  * @file vcs_setMolesLinProg.cpp
3  */
4 
5 // This file is part of Cantera. See License.txt in the top-level directory or
6 // at https://cantera.org/license.txt for license and copyright information.
7 
9 
10 using namespace std;
11 
12 namespace Cantera
13 {
14 
15 static void printProgress(const vector<string> &spName,
16  const vector_fp &soln,
17  const vector_fp &ff)
18 {
19  double sum = 0.0;
20  plogf(" --- Summary of current progress:\n");
21  plogf(" --- Name Moles - SSGibbs \n");
22  plogf(" -------------------------------------------------------------------------------------\n");
23  for (size_t k = 0; k < soln.size(); k++) {
24  plogf(" --- %20s %12.4g - %12.4g\n", spName[k], soln[k], ff[k]);
25  sum += soln[k] * ff[k];
26  }
27  plogf(" --- Total sum to be minimized = %g\n", sum);
28 }
29 
30 int VCS_SOLVE::vcs_setMolesLinProg()
31 {
32  double test = -1.0E-10;
33 
34  if (m_debug_print_lvl >= 2) {
35  plogf(" --- call setInitialMoles\n");
36  }
37 
38  double dxi_min = 1.0e10;
39  int retn;
40  int iter = 0;
41  bool abundancesOK = true;
42  bool usedZeroedSpecies;
43  vector_fp sm(m_nelem * m_nelem, 0.0);
44  vector_fp ss(m_nelem, 0.0);
45  vector_fp sa(m_nelem, 0.0);
46  vector_fp wx(m_nelem, 0.0);
47  vector_fp aw(m_nsp, 0.0);
48 
49  for (size_t ik = 0; ik < m_nsp; ik++) {
50  if (m_speciesUnknownType[ik] != VCS_SPECIES_INTERFACIALVOLTAGE) {
51  m_molNumSpecies_old[ik] = max(0.0, m_molNumSpecies_old[ik]);
52  }
53  }
54 
55  if (m_debug_print_lvl >= 2) {
56  printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies);
57  }
58 
59  bool redo = true;
60  while (redo) {
61  if (!vcs_elabcheck(0)) {
62  if (m_debug_print_lvl >= 2) {
63  plogf(" --- seMolesLinProg Mole numbers failing element abundances\n");
64  plogf(" --- seMolesLinProg Call vcs_elcorr to attempt fix\n");
65  }
66  retn = vcs_elcorr(&sm[0], &wx[0]);
67  if (retn >= 2) {
68  abundancesOK = false;
69  } else {
70  abundancesOK = true;
71  }
72  } else {
73  abundancesOK = true;
74  }
75 
76  // Now find the optimized basis that spans the stoichiometric
77  // coefficient matrix, based on the current composition,
78  // m_molNumSpecies_old[] We also calculate sc[][], the reaction matrix.
79  retn = vcs_basopt(false, &aw[0], &sa[0], &sm[0], &ss[0],
80  test, &usedZeroedSpecies);
81  if (retn != VCS_SUCCESS) {
82  return retn;
83  }
84 
85  if (m_debug_print_lvl >= 2) {
86  plogf("iteration %d\n", iter);
87  }
88  redo = false;
89  iter++;
90  if (iter > 15) {
91  break;
92  }
93 
94  // loop over all reactions
95  for (size_t irxn = 0; irxn < m_numRxnTot; irxn++) {
96  // dg_rt is the Delta_G / RT value for the reaction
97  size_t ik = m_numComponents + irxn;
98  double dg_rt = m_SSfeSpecies[ik];
99  dxi_min = 1.0e10;
100  const double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
101  for (size_t jcomp = 0; jcomp < m_nelem; jcomp++) {
102  dg_rt += m_SSfeSpecies[jcomp] * sc_irxn[jcomp];
103  }
104  // fwd or rev direction.
105  // idir > 0 implies increasing the current species
106  // idir < 0 implies decreasing the current species
107  int idir = (dg_rt < 0.0 ? 1 : -1);
108  if (idir < 0) {
109  dxi_min = m_molNumSpecies_old[ik];
110  }
111 
112  for (size_t jcomp = 0; jcomp < m_numComponents; jcomp++) {
113  double nu = sc_irxn[jcomp];
114  // set max change in progress variable by
115  // non-negativity requirement
116  if (nu*idir < 0) {
117  double delta_xi = fabs(m_molNumSpecies_old[jcomp]/nu);
118  // if a component has nearly zero moles, redo
119  // with a new set of components
120  if (!redo && delta_xi < 1.0e-10 && (m_molNumSpecies_old[ik] >= 1.0E-10)) {
121  if (m_debug_print_lvl >= 2) {
122  plogf(" --- Component too small: %s\n", m_speciesName[jcomp]);
123  }
124  redo = true;
125  }
126  dxi_min = std::min(dxi_min, delta_xi);
127  }
128  }
129 
130  // step the composition by dxi_min, check against zero, since
131  // we are zeroing components and species on every step.
132  // Redo the iteration, if a component went from positive to zero on this step.
133  double dsLocal = idir*dxi_min;
134  m_molNumSpecies_old[ik] += dsLocal;
135  m_molNumSpecies_old[ik] = max(0.0, m_molNumSpecies_old[ik]);
136  for (size_t jcomp = 0; jcomp < m_numComponents; jcomp++) {
137  bool full = false;
138  if (m_molNumSpecies_old[jcomp] > 1.0E-15) {
139  full = true;
140  }
141  m_molNumSpecies_old[jcomp] += sc_irxn[jcomp] * dsLocal;
142  m_molNumSpecies_old[jcomp] = max(0.0, m_molNumSpecies_old[jcomp]);
143  if (full && m_molNumSpecies_old[jcomp] < 1.0E-60) {
144  redo = true;
145  }
146  }
147  }
148 
149  if (m_debug_print_lvl >= 2) {
150  printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies);
151  }
152  }
153 
154  if (m_debug_print_lvl == 1) {
155  printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies);
156  plogf(" --- setInitialMoles end\n");
157  }
158  retn = 0;
159  if (!abundancesOK) {
160  retn = -1;
161  } else if (iter > 15) {
162  retn = 1;
163  }
164  return retn;
165 }
166 
167 }
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:180
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264
#define VCS_SPECIES_INTERFACIALVOLTAGE
Species refers to an electron in the metal.
Definition: vcs_defs.h:145
#define VCS_SUCCESS
Definition: vcs_defs.h:18
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:18
Header file for the internal object that holds the vcs equilibrium problem (see Class VCS_SOLVE and E...