Cantera  2.4.0
vcs_elem.cpp
Go to the documentation of this file.
1 /**
2  * @file vcs_elem.cpp
3  * This file contains the algorithm for checking the satisfaction of the
4  * element abundances constraints and the algorithm for fixing violations
5  * of the element abundances constraints.
6  */
7 
8 // This file is part of Cantera. See License.txt in the top-level directory or
9 // at http://www.cantera.org/license.txt for license and copyright information.
10 
14 
15 namespace Cantera
16 {
18 {
19  for (size_t j = 0; j < m_nelem; ++j) {
20  m_elemAbundances[j] = 0.0;
21  for (size_t i = 0; i < m_nsp; ++i) {
24  }
25  }
26  }
27 }
28 
29 bool VCS_SOLVE::vcs_elabcheck(int ibound)
30 {
31  size_t top = m_numComponents;
32  if (ibound) {
33  top = m_nelem;
34  }
35 
36  for (size_t i = 0; i < top; ++i) {
37  // Require 12 digits of accuracy on non-zero constraints.
38  if (m_elementActive[i] && fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > fabs(m_elemAbundancesGoal[i]) * 1.0e-12) {
39  // This logic is for charge neutrality condition
41  m_elemAbundancesGoal[i] != 0.0) {
42  throw CanteraError("VCS_SOLVE::vcs_elabcheck",
43  "Problem with charge neutrality condition");
44  }
47 
48  // Find out if the constraint is a multisign constraint. If it
49  // is, then we have to worry about roundoff error in the
50  // addition of terms. We are limited to 13 digits of finite
51  // arithmetic accuracy.
52  bool multisign = false;
53  for (size_t kspec = 0; kspec < m_nsp; kspec++) {
54  double eval = m_formulaMatrix(kspec,i);
55  if (eval < 0.0) {
56  multisign = true;
57  }
58  if (eval != 0.0) {
59  scale = std::max(scale, fabs(eval * m_molNumSpecies_old[kspec]));
60  }
61  }
62  if (multisign) {
63  if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > 1e-11 * scale) {
64  return false;
65  }
66  } else {
68  return false;
69  }
70  }
71  } else {
72  // For normal element balances, we require absolute compliance
73  // even for ridiculously small numbers.
74  if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) {
75  return false;
76  } else {
77  return false;
78  }
79  }
80  }
81  }
82  return true;
83 }
84 
85 void VCS_SOLVE::vcs_elabPhase(size_t iphase, double* const elemAbundPhase)
86 {
87  for (size_t j = 0; j < m_nelem; ++j) {
88  elemAbundPhase[j] = 0.0;
89  for (size_t i = 0; i < m_nsp; ++i) {
91  elemAbundPhase[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i];
92  }
93  }
94  }
95 }
96 
97 int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
98 {
99  int retn = 0;
100 
101  vector_fp ga_save(m_elemAbundances);
102  if (m_debug_print_lvl >= 2) {
103  plogf(" --- vcsc_elcorr: Element abundances correction routine");
104  if (m_nelem != m_numComponents) {
105  plogf(" (m_numComponents != m_nelem)");
106  }
107  plogf("\n");
108  }
109 
110  for (size_t i = 0; i < m_nelem; ++i) {
111  x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i];
112  }
113  double l2before = 0.0;
114  for (size_t i = 0; i < m_nelem; ++i) {
115  l2before += x[i] * x[i];
116  }
117  l2before = sqrt(l2before/m_nelem);
118 
119  // Special section to take out single species, single component,
120  // moles. These are species which have non-zero entries in the
121  // formula matrix, and no other species have zero values either.
122  bool changed = false;
123  for (size_t i = 0; i < m_nelem; ++i) {
124  int numNonZero = 0;
125  bool multisign = false;
126  for (size_t kspec = 0; kspec < m_nsp; kspec++) {
128  double eval = m_formulaMatrix(kspec,i);
129  if (eval < 0.0) {
130  multisign = true;
131  }
132  if (eval != 0.0) {
133  numNonZero++;
134  }
135  }
136  }
137  if (!multisign) {
138  if (numNonZero < 2) {
139  for (size_t kspec = 0; kspec < m_nsp; kspec++) {
141  double eval = m_formulaMatrix(kspec,i);
142  if (eval > 0.0) {
143  m_molNumSpecies_old[kspec] = m_elemAbundancesGoal[i] / eval;
144  changed = true;
145  }
146  }
147  }
148  } else {
149  int numCompNonZero = 0;
150  size_t compID = npos;
151  for (size_t kspec = 0; kspec < m_numComponents; kspec++) {
153  double eval = m_formulaMatrix(kspec,i);
154  if (eval > 0.0) {
155  compID = kspec;
156  numCompNonZero++;
157  }
158  }
159  }
160  if (numCompNonZero == 1) {
161  double diff = m_elemAbundancesGoal[i];
162  for (size_t kspec = m_numComponents; kspec < m_nsp; kspec++) {
164  double eval = m_formulaMatrix(kspec,i);
165  diff -= eval * m_molNumSpecies_old[kspec];
166  }
167  m_molNumSpecies_old[compID] = std::max(0.0,diff/m_formulaMatrix(compID,i));
168  changed = true;
169  }
170  }
171  }
172  }
173  }
174  if (changed) {
175  vcs_elab();
176  }
177 
178  // Section to check for maximum bounds errors on all species due to
179  // elements. This may only be tried on element types which are
180  // VCS_ELEM_TYPE_ABSPOS. This is because no other species may have a
181  // negative number of these.
182  //
183  // Note, also we can do this over ne, the number of elements, not just the
184  // number of components.
185  changed = false;
186  for (size_t i = 0; i < m_nelem; ++i) {
187  int elType = m_elType[i];
188  if (elType == VCS_ELEM_TYPE_ABSPOS) {
189  for (size_t kspec = 0; kspec < m_nsp; kspec++) {
191  double atomComp = m_formulaMatrix(kspec,i);
192  if (atomComp > 0.0) {
193  double maxPermissible = m_elemAbundancesGoal[i] / atomComp;
194  if (m_molNumSpecies_old[kspec] > maxPermissible) {
195  if (m_debug_print_lvl >= 3) {
196  plogf(" --- vcs_elcorr: Reduced species %s from %g to %g "
197  "due to %s max bounds constraint\n",
198  m_speciesName[kspec], m_molNumSpecies_old[kspec],
199  maxPermissible, m_elementName[i]);
200  }
201  m_molNumSpecies_old[kspec] = maxPermissible;
202  changed = true;
204  m_molNumSpecies_old[kspec] = 0.0;
205  if (m_SSPhase[kspec]) {
207  } else {
209  }
210  if (m_debug_print_lvl >= 2) {
211  plogf(" --- vcs_elcorr: Zeroed species %s and changed "
212  "status to %d due to max bounds constraint\n",
213  m_speciesName[kspec], m_speciesStatus[kspec]);
214  }
215  }
216  }
217  }
218  }
219  }
220  }
221  }
222 
223  // Recalculate the element abundances if something has changed.
224  if (changed) {
225  vcs_elab();
226  }
227 
228  // Ok, do the general case. Linear algebra problem is of length nc, not ne,
229  // as there may be degenerate rows when nc .ne. ne.
231  for (size_t i = 0; i < m_numComponents; ++i) {
232  x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i];
233  if (fabs(x[i]) > 1.0E-13) {
234  retn = 1;
235  }
236  for (size_t j = 0; j < m_numComponents; ++j) {
237  A(j, i) = - m_formulaMatrix(i,j);
238  }
239  }
240 
241  solve(A, x, 1, m_nelem);
242 
243  // Now apply the new direction without creating negative species.
244  double par = 0.5;
245  for (size_t i = 0; i < m_numComponents; ++i) {
246  if (m_molNumSpecies_old[i] > 0.0) {
247  par = std::max(par, -x[i] / m_molNumSpecies_old[i]);
248  }
249  }
250  par = std::min(par, 100.0);
251  par = 1.0 / par;
252  if (par < 1.0 && par > 0.0) {
253  retn = 2;
254  par *= 0.9999;
255  for (size_t i = 0; i < m_numComponents; ++i) {
256  double tmp = m_molNumSpecies_old[i] + par * x[i];
257  if (tmp > 0.0) {
258  m_molNumSpecies_old[i] = tmp;
259  } else {
260  if (m_SSPhase[i]) {
261  m_molNumSpecies_old[i] = 0.0;
262  } else {
263  m_molNumSpecies_old[i] = m_molNumSpecies_old[i] * 0.0001;
264  }
265  }
266  }
267  } else {
268  for (size_t i = 0; i < m_numComponents; ++i) {
269  double tmp = m_molNumSpecies_old[i] + x[i];
270  if (tmp > 0.0) {
271  m_molNumSpecies_old[i] = tmp;
272  } else {
273  if (m_SSPhase[i]) {
274  m_molNumSpecies_old[i] = 0.0;
275  } else {
276  m_molNumSpecies_old[i] = m_molNumSpecies_old[i] * 0.0001;
277  }
278  }
279  }
280  }
281 
282  // We have changed the element abundances. Calculate them again
283  vcs_elab();
284 
285  // We have changed the total moles in each phase. Calculate them again
286  vcs_tmoles();
287 
288  // Try some ad hoc procedures for fixing the problem
289  if (retn >= 2) {
290  // First find a species whose adjustment is a win-win situation.
291  for (size_t kspec = 0; kspec < m_nsp; kspec++) {
293  continue;
294  }
295  double saveDir = 0.0;
296  bool goodSpec = true;
297  for (size_t i = 0; i < m_numComponents; ++i) {
298  double dir = m_formulaMatrix(kspec,i) * (m_elemAbundancesGoal[i] - m_elemAbundances[i]);
299  if (fabs(dir) > 1.0E-10) {
300  if (dir > 0.0) {
301  if (saveDir < 0.0) {
302  goodSpec = false;
303  break;
304  }
305  } else {
306  if (saveDir > 0.0) {
307  goodSpec = false;
308  break;
309  }
310  }
311  saveDir = dir;
312  } else {
313  if (m_formulaMatrix(kspec,i) != 0.) {
314  goodSpec = false;
315  break;
316  }
317  }
318  }
319  if (goodSpec) {
320  int its = 0;
321  double xx = 0.0;
322  for (size_t i = 0; i < m_numComponents; ++i) {
323  if (m_formulaMatrix(kspec,i) != 0.0) {
324  xx += (m_elemAbundancesGoal[i] - m_elemAbundances[i]) / m_formulaMatrix(kspec,i);
325  its++;
326  }
327  }
328  if (its > 0) {
329  xx /= its;
330  }
331  m_molNumSpecies_old[kspec] += xx;
332  m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 1.0E-10);
333 
334  // If we are dealing with a deleted species, then we need to
335  // reinsert it into the active list.
336  if (kspec >= m_numSpeciesRdc) {
337  vcs_reinsert_deleted(kspec);
339  vcs_elab();
340  goto L_CLEANUP;
341  }
342  vcs_elab();
343  }
344  }
345  }
346  if (vcs_elabcheck(0)) {
347  retn = 1;
348  goto L_CLEANUP;
349  }
350 
351  for (size_t i = 0; i < m_nelem; ++i) {
353  (m_elType[i] == VCS_ELEM_TYPE_ABSPOS && m_elemAbundancesGoal[i] == 0.0)) {
354  for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
355  if (m_elemAbundances[i] > 0.0 && m_formulaMatrix(kspec,i) < 0.0) {
356  m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix(kspec,i);
357  m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0);
358  vcs_elab();
359  break;
360  }
361  if (m_elemAbundances[i] < 0.0 && m_formulaMatrix(kspec,i) > 0.0) {
362  m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix(kspec,i);
363  m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0);
364  vcs_elab();
365  break;
366  }
367  }
368  }
369  }
370  if (vcs_elabcheck(1)) {
371  retn = 1;
372  goto L_CLEANUP;
373  }
374 
375  // For electron charges element types, we try positive deltas in the species
376  // concentrations to match the desired electron charge exactly.
377  for (size_t i = 0; i < m_nelem; ++i) {
378  double dev = m_elemAbundancesGoal[i] - m_elemAbundances[i];
379  if (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE && (fabs(dev) > 1.0E-300)) {
380  bool useZeroed = true;
381  for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
382  if (dev < 0.0) {
383  if (m_formulaMatrix(kspec,i) < 0.0 && m_molNumSpecies_old[kspec] > 0.0) {
384  useZeroed = false;
385  }
386  } else {
387  if (m_formulaMatrix(kspec,i) > 0.0 && m_molNumSpecies_old[kspec] > 0.0) {
388  useZeroed = false;
389  }
390  }
391  }
392  for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
393  if (m_molNumSpecies_old[kspec] > 0.0 || useZeroed) {
394  if (dev < 0.0 && m_formulaMatrix(kspec,i) < 0.0) {
395  double delta = dev / m_formulaMatrix(kspec,i);
396  m_molNumSpecies_old[kspec] += delta;
397  m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0);
398  vcs_elab();
399  break;
400  }
401  if (dev > 0.0 && m_formulaMatrix(kspec,i) > 0.0) {
402  double delta = dev / m_formulaMatrix(kspec,i);
403  m_molNumSpecies_old[kspec] += delta;
404  m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0);
405  vcs_elab();
406  break;
407  }
408  }
409  }
410  }
411  }
412  if (vcs_elabcheck(1)) {
413  retn = 1;
414  goto L_CLEANUP;
415  }
416 
417 L_CLEANUP:
418  ;
419  vcs_tmoles();
420  double l2after = 0.0;
421  for (size_t i = 0; i < m_nelem; ++i) {
422  l2after += pow(m_elemAbundances[i] - m_elemAbundancesGoal[i], 2);
423  }
424  l2after = sqrt(l2after/m_nelem);
425  if (m_debug_print_lvl >= 2) {
426  plogf(" --- Elem_Abund: Correct Initial "
427  " Final\n");
428  for (size_t i = 0; i < m_nelem; ++i) {
429  plogf(" --- ");
430  plogf("%-2.2s", m_elementName[i]);
431  plogf(" %20.12E %20.12E %20.12E\n", m_elemAbundancesGoal[i], ga_save[i], m_elemAbundances[i]);
432  }
433  plogf(" --- Diff_Norm: %20.12E %20.12E\n",
434  l2before, l2after);
435  }
436  return retn;
437 }
438 
439 }
#define VCS_ELEM_TYPE_ELECTRONCHARGE
This refers to conservation of electrons.
Definition: vcs_defs.h:234
std::vector< size_t > m_phaseID
Mapping from the species number to the phase number.
Definition: vcs_solve.h:1357
vector_int m_speciesUnknownType
Specifies the species unknown type.
Definition: vcs_solve.h:1165
vector_fp m_elemAbundances
Element abundances vector.
Definition: vcs_solve.h:1222
vector_fp m_molNumSpecies_old
Total moles of the species.
Definition: vcs_solve.h:1152
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
int vcs_elcorr(double aa[], double x[])
Definition: vcs_elem.cpp:97
void vcs_elabPhase(size_t iphase, double *const elemAbundPhase)
Definition: vcs_elem.cpp:85
double vcs_tmoles()
Calculates the total number of moles of species in all phases.
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
size_t m_numSpeciesRdc
Current number of species in the problems.
Definition: vcs_solve.h:1057
vector_int m_elementActive
Specifies whether an element constraint is active.
Definition: vcs_solve.h:1391
std::vector< std::string > m_elementName
Vector of strings containing the element names.
Definition: vcs_solve.h:1370
Header file for the internal object that holds the vcs equilibrium problem (see Class VCS_SOLVE and E...
#define VCS_SPECIES_ZEROEDSS
Species is a SS phase, that is currently zeroed out.
Definition: vcs_defs.h:129
#define VCS_ELEM_TYPE_ABSPOS
Normal element constraint consisting of positive coefficients for the formula matrix.
Definition: vcs_defs.h:228
#define VCS_DELETE_MINORSPECIES_CUTOFF
Cutoff relative mole number value, below which species are deleted from the equilibrium problem...
Definition: vcs_defs.h:44
#define VCS_ELEM_TYPE_CHARGENEUTRALITY
This refers to a charge neutrality of a single phase.
Definition: vcs_defs.h:240
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
#define VCS_SPECIES_ACTIVEBUTZERO
Species lies in a multicomponent phase that is active, but species concentration is zero...
Definition: vcs_defs.h:163
int m_debug_print_lvl
Debug printing lvl.
Definition: vcs_solve.h:1498
size_t m_numComponents
Number of components calculated for the problem.
Definition: vcs_solve.h:1050
vector_fp m_elemAbundancesGoal
Element abundances vector Goals.
Definition: vcs_solve.h:1231
void vcs_reinsert_deleted(size_t kspec)
std::vector< char > m_SSPhase
Boolean indicating whether a species belongs to a single-species phase.
Definition: vcs_solve.h:1361
std::vector< std::string > m_speciesName
Species string name for the kth species.
Definition: vcs_solve.h:1364
Array2D m_formulaMatrix
Formula matrix for the problem.
Definition: vcs_solve.h:1076
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:157
void vcs_elab()
Computes the current elemental abundances vector.
Definition: vcs_elem.cpp:17
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:130
size_t m_nsp
Total number of species in the problems.
Definition: vcs_solve.h:1044
Headers for the DenseMatrix object, which deals with dense rectangular matrices and description of th...
vector_int m_speciesStatus
Major -Minor status vector for the species in the problem.
Definition: vcs_solve.h:1354
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:18
#define VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
Unknown refers to the voltage level of a phase.
Definition: vcs_defs.h:289
vector_int m_elType
Type of the element constraint.
Definition: vcs_solve.h:1385
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
size_t m_nelem
Number of element constraints in the problem.
Definition: vcs_solve.h:1047
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition: DenseMatrix.h:50
bool vcs_elabcheck(int ibound)
Definition: vcs_elem.cpp:29