Cantera  2.1.2
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  */
10 #include "math.h"
11 
12 namespace VCSnonideal
13 {
15 {
16  for (size_t j = 0; j < m_numElemConstraints; ++j) {
17  m_elemAbundances[j] = 0.0;
18  for (size_t i = 0; i < m_numSpeciesTot; ++i) {
21  }
22  }
23  }
24 }
25 
26 bool VCS_SOLVE::vcs_elabcheck(int ibound)
27 {
28  size_t top = m_numComponents;
29  double eval, scale;
30  int numNonZero;
31  bool multisign = false;
32  if (ibound) {
34  }
35  /*
36  * Require 12 digits of accuracy on non-zero constraints.
37  */
38  for (size_t i = 0; i < top; ++i) {
39  if (m_elementActive[i]) {
40  if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > (fabs(m_elemAbundancesGoal[i]) * 1.0e-12)) {
41  /*
42  * This logic is for charge neutrality condition
43  */
45  m_elemAbundancesGoal[i] != 0.0) {
46  throw Cantera::CanteraError("VCS_SOLVE::vcs_elabcheck",
47  "Problem with charge neutrality condition");
48  }
51  /*
52  * Find out if the constraint is a multisign constraint.
53  * If it is, then we have to worry about roundoff error
54  * in the addition of terms. We are limited to 13
55  * digits of finite arithmetic accuracy.
56  */
57  numNonZero = 0;
58  multisign = false;
59  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
60  eval = m_formulaMatrix[i][kspec];
61  if (eval < 0.0) {
62  multisign = true;
63  }
64  if (eval != 0.0) {
65  scale = std::max(scale, fabs(eval * m_molNumSpecies_old[kspec]));
66  numNonZero++;
67  }
68  }
69  if (multisign) {
70  if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > 1e-11 * scale) {
71  return false;
72  }
73  } else {
75  return false;
76  }
77  }
78  } else {
79  /*
80  * For normal element balances, we require absolute compliance
81  * even for ridiculously small numbers.
82  */
83  if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) {
84  return false;
85  } else {
86  return false;
87  }
88  }
89  }
90  }
91  }
92  return true;
93 }
94 
95 void VCS_SOLVE::vcs_elabPhase(size_t iphase, double* const elemAbundPhase)
96 {
97  for (size_t j = 0; j < m_numElemConstraints; ++j) {
98  elemAbundPhase[j] = 0.0;
99  for (size_t i = 0; i < m_numSpeciesTot; ++i) {
101  if (m_phaseID[i] == iphase) {
102  elemAbundPhase[j] += m_formulaMatrix[j][i] * m_molNumSpecies_old[i];
103  }
104  }
105  }
106  }
107 }
108 
109 int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
110 {
111  int retn = 0, its;
112  bool goodSpec;
113  double xx, par, saveDir, dir;
114 
115 #ifdef DEBUG_MODE
116  double l2before = 0.0, l2after = 0.0;
117  std::vector<double> ga_save(m_numElemConstraints, 0.0);
119  if (m_debug_print_lvl >= 2) {
120  plogf(" --- vcsc_elcorr: Element abundances correction routine");
122  plogf(" (m_numComponents != m_numElemConstraints)");
123  }
124  plogf("\n");
125  }
126 
127  for (size_t i = 0; i < m_numElemConstraints; ++i) {
128  x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i];
129  }
130  l2before = 0.0;
131  for (size_t i = 0; i < m_numElemConstraints; ++i) {
132  l2before += x[i] * x[i];
133  }
134  l2before = sqrt(l2before/m_numElemConstraints);
135 #endif
136 
137  /*
138  * Special section to take out single species, single component,
139  * moles. These are species which have non-zero entries in the
140  * formula matrix, and no other species have zero values either.
141  *
142  */
143  int numNonZero = 0;
144  bool changed = false;
145  bool multisign = false;
146  for (size_t i = 0; i < m_numElemConstraints; ++i) {
147  numNonZero = 0;
148  multisign = false;
149  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
151  double eval = m_formulaMatrix[i][kspec];
152  if (eval < 0.0) {
153  multisign = true;
154  }
155  if (eval != 0.0) {
156  numNonZero++;
157  }
158  }
159  }
160  if (!multisign) {
161  if (numNonZero < 2) {
162  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
164  double eval = m_formulaMatrix[i][kspec];
165  if (eval > 0.0) {
166  m_molNumSpecies_old[kspec] = m_elemAbundancesGoal[i] / eval;
167  changed = true;
168  }
169  }
170  }
171  } else {
172  int numCompNonZero = 0;
173  size_t compID = npos;
174  for (size_t kspec = 0; kspec < m_numComponents; kspec++) {
176  double eval = m_formulaMatrix[i][kspec];
177  if (eval > 0.0) {
178  compID = kspec;
179  numCompNonZero++;
180  }
181  }
182  }
183  if (numCompNonZero == 1) {
184  double diff = m_elemAbundancesGoal[i];
185  for (size_t kspec = m_numComponents; kspec < m_numSpeciesTot; kspec++) {
187  double eval = m_formulaMatrix[i][kspec];
188  diff -= eval * m_molNumSpecies_old[kspec];
189  }
190  m_molNumSpecies_old[compID] = std::max(0.0,diff/m_formulaMatrix[i][compID]);
191  changed = true;
192  }
193  }
194  }
195  }
196  }
197  if (changed) {
198  vcs_elab();
199  }
200 
201  /*
202  * Section to check for maximum bounds errors on all species
203  * due to elements.
204  * This may only be tried on element types which are VCS_ELEM_TYPE_ABSPOS.
205  * This is because no other species may have a negative number of these.
206  *
207  * Note, also we can do this over ne, the number of elements, not just
208  * the number of components.
209  */
210  changed = false;
211  for (size_t i = 0; i < m_numElemConstraints; ++i) {
212  int elType = m_elType[i];
213  if (elType == VCS_ELEM_TYPE_ABSPOS) {
214  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
216  double atomComp = m_formulaMatrix[i][kspec];
217  if (atomComp > 0.0) {
218  double maxPermissible = m_elemAbundancesGoal[i] / atomComp;
219  if (m_molNumSpecies_old[kspec] > maxPermissible) {
220 
221 #ifdef DEBUG_MODE
222  if (m_debug_print_lvl >= 3) {
223  plogf(" --- vcs_elcorr: Reduced species %s from %g to %g "
224  "due to %s max bounds constraint\n",
225  m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec],
226  maxPermissible, m_elementName[i].c_str());
227  }
228 #endif
229  m_molNumSpecies_old[kspec] = maxPermissible;
230  changed = true;
232  m_molNumSpecies_old[kspec] = 0.0;
233  if (m_SSPhase[kspec]) {
235  } else {
237  }
238 #ifdef DEBUG_MODE
239  if (m_debug_print_lvl >= 2) {
240  plogf(" --- vcs_elcorr: Zeroed species %s and changed "
241  "status to %d due to max bounds constraint\n",
242  m_speciesName[kspec].c_str(), m_speciesStatus[kspec]);
243  }
244 #endif
245  }
246  }
247  }
248  }
249  }
250  }
251  }
252 
253  // Recalculate the element abundances if something has changed.
254  if (changed) {
255  vcs_elab();
256  }
257 
258  /*
259  * Ok, do the general case. Linear algebra problem is
260  * of length nc, not ne, as there may be degenerate rows when
261  * nc .ne. ne.
262  */
263  for (size_t i = 0; i < m_numComponents; ++i) {
264  x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i];
265  if (fabs(x[i]) > 1.0E-13) {
266  retn = 1;
267  }
268  for (size_t j = 0; j < m_numComponents; ++j) {
269  aa[j + i*m_numElemConstraints] = m_formulaMatrix[j][i];
270  }
271  }
272  int err = vcsUtil_mlequ(aa, m_numElemConstraints, m_numComponents, x, 1);
273  if (err == 1) {
274  plogf("vcs_elcorr ERROR: mlequ returned error condition\n");
275  return VCS_FAILED_CONVERGENCE;
276  }
277  /*
278  * Now apply the new direction without creating negative species.
279  */
280  par = 0.5;
281  for (size_t i = 0; i < m_numComponents; ++i) {
282  if (m_molNumSpecies_old[i] > 0.0) {
283  xx = -x[i] / m_molNumSpecies_old[i];
284  if (par < xx) {
285  par = xx;
286  }
287  }
288  }
289  if (par > 100.0) {
290  par = 100.0;
291  }
292  par = 1.0 / par;
293  if (par < 1.0 && par > 0.0) {
294  retn = 2;
295  par *= 0.9999;
296  for (size_t i = 0; i < m_numComponents; ++i) {
297  double tmp = m_molNumSpecies_old[i] + par * x[i];
298  if (tmp > 0.0) {
299  m_molNumSpecies_old[i] = tmp;
300  } else {
301  if (m_SSPhase[i]) {
302  m_molNumSpecies_old[i] = 0.0;
303  } else {
304  m_molNumSpecies_old[i] = m_molNumSpecies_old[i] * 0.0001;
305  }
306  }
307  }
308  } else {
309  for (size_t i = 0; i < m_numComponents; ++i) {
310  double tmp = m_molNumSpecies_old[i] + x[i];
311  if (tmp > 0.0) {
312  m_molNumSpecies_old[i] = tmp;
313  } else {
314  if (m_SSPhase[i]) {
315  m_molNumSpecies_old[i] = 0.0;
316  } else {
317  m_molNumSpecies_old[i] = m_molNumSpecies_old[i] * 0.0001;
318  }
319  }
320  }
321  }
322 
323  /*
324  * We have changed the element abundances. Calculate them again
325  */
326  vcs_elab();
327  /*
328  * We have changed the total moles in each phase. Calculate them again
329  */
330  vcs_tmoles();
331 
332  /*
333  * Try some ad hoc procedures for fixing the problem
334  */
335  if (retn >= 2) {
336  /*
337  * First find a species whose adjustment is a win-win
338  * situation.
339  */
340  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
342  continue;
343  }
344  saveDir = 0.0;
345  goodSpec = true;
346  for (size_t i = 0; i < m_numComponents; ++i) {
347  dir = m_formulaMatrix[i][kspec] * (m_elemAbundancesGoal[i] - m_elemAbundances[i]);
348  if (fabs(dir) > 1.0E-10) {
349  if (dir > 0.0) {
350  if (saveDir < 0.0) {
351  goodSpec = false;
352  break;
353  }
354  } else {
355  if (saveDir > 0.0) {
356  goodSpec = false;
357  break;
358  }
359  }
360  saveDir = dir;
361  } else {
362  if (m_formulaMatrix[i][kspec] != 0.) {
363  goodSpec = false;
364  break;
365  }
366  }
367  }
368  if (goodSpec) {
369  its = 0;
370  xx = 0.0;
371  for (size_t i = 0; i < m_numComponents; ++i) {
372  if (m_formulaMatrix[i][kspec] != 0.0) {
373  xx += (m_elemAbundancesGoal[i] - m_elemAbundances[i]) / m_formulaMatrix[i][kspec];
374  its++;
375  }
376  }
377  if (its > 0) {
378  xx /= its;
379  }
380  m_molNumSpecies_old[kspec] += xx;
381  m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 1.0E-10);
382  /*
383  * If we are dealing with a deleted species, then
384  * we need to reinsert it into the active list.
385  */
386  if (kspec >= m_numSpeciesRdc) {
387  vcs_reinsert_deleted(kspec);
389  vcs_elab();
390  goto L_CLEANUP;
391  }
392  vcs_elab();
393  }
394  }
395  }
396  if (vcs_elabcheck(0)) {
397  retn = 1;
398  goto L_CLEANUP;
399  }
400 
401  for (size_t i = 0; i < m_numElemConstraints; ++i) {
403  (m_elType[i] == VCS_ELEM_TYPE_ABSPOS && m_elemAbundancesGoal[i] == 0.0)) {
404  for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
405  if (m_elemAbundances[i] > 0.0) {
406  if (m_formulaMatrix[i][kspec] < 0.0) {
407  m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix[i][kspec] ;
408  if (m_molNumSpecies_old[kspec] < 0.0) {
409  m_molNumSpecies_old[kspec] = 0.0;
410  }
411  vcs_elab();
412  break;
413  }
414  }
415  if (m_elemAbundances[i] < 0.0) {
416  if (m_formulaMatrix[i][kspec] > 0.0) {
417  m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix[i][kspec];
418  if (m_molNumSpecies_old[kspec] < 0.0) {
419  m_molNumSpecies_old[kspec] = 0.0;
420  }
421  vcs_elab();
422  break;
423  }
424  }
425  }
426  }
427  }
428  if (vcs_elabcheck(1)) {
429  retn = 1;
430  goto L_CLEANUP;
431  }
432 
433  /*
434  * For electron charges element types, we try positive deltas
435  * in the species concentrations to match the desired
436  * electron charge exactly.
437  */
438  for (size_t i = 0; i < m_numElemConstraints; ++i) {
439  double dev = m_elemAbundancesGoal[i] - m_elemAbundances[i];
440  if (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE && (fabs(dev) > 1.0E-300)) {
441  bool useZeroed = true;
442  for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
443  if (dev < 0.0) {
444  if (m_formulaMatrix[i][kspec] < 0.0) {
445  if (m_molNumSpecies_old[kspec] > 0.0) {
446  useZeroed = false;
447  }
448  }
449  } else {
450  if (m_formulaMatrix[i][kspec] > 0.0) {
451  if (m_molNumSpecies_old[kspec] > 0.0) {
452  useZeroed = false;
453  }
454  }
455  }
456  }
457  for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
458  if (m_molNumSpecies_old[kspec] > 0.0 || useZeroed) {
459  if (dev < 0.0) {
460  if (m_formulaMatrix[i][kspec] < 0.0) {
461  double delta = dev / m_formulaMatrix[i][kspec] ;
462  m_molNumSpecies_old[kspec] += delta;
463  if (m_molNumSpecies_old[kspec] < 0.0) {
464  m_molNumSpecies_old[kspec] = 0.0;
465  }
466  vcs_elab();
467  break;
468  }
469  }
470  if (dev > 0.0) {
471  if (m_formulaMatrix[i][kspec] > 0.0) {
472  double delta = dev / m_formulaMatrix[i][kspec] ;
473  m_molNumSpecies_old[kspec] += delta;
474  if (m_molNumSpecies_old[kspec] < 0.0) {
475  m_molNumSpecies_old[kspec] = 0.0;
476  }
477  vcs_elab();
478  break;
479  }
480  }
481  }
482  }
483  }
484  }
485  if (vcs_elabcheck(1)) {
486  retn = 1;
487  goto L_CLEANUP;
488  }
489 
490 L_CLEANUP:
491  ;
492  vcs_tmoles();
493 #ifdef DEBUG_MODE
494  l2after = 0.0;
495  for (int i = 0; i < m_numElemConstraints; ++i) {
496  l2after += SQUARE(m_elemAbundances[i] - m_elemAbundancesGoal[i]);
497  }
498  l2after = sqrt(l2after/m_numElemConstraints);
499  if (m_debug_print_lvl >= 2) {
500  plogf(" --- Elem_Abund: Correct Initial "
501  " Final\n");
502  for (int i = 0; i < m_numElemConstraints; ++i) {
503  plogf(" --- ");
504  plogf("%-2.2s", m_elementName[i].c_str());
505  plogf(" %20.12E %20.12E %20.12E\n", m_elemAbundancesGoal[i], ga_save[i], m_elemAbundances[i]);
506  }
507  plogf(" --- Diff_Norm: %20.12E %20.12E\n",
508  l2before, l2after);
509  }
510 #endif
511  return retn;
512 }
513 
514 }
#define VCS_ELEM_TYPE_ELECTRONCHARGE
This refers to conservation of electrons.
Definition: vcs_defs.h:312
void vcs_elab()
Computes the current elemental abundances vector.
Definition: vcs_elem.cpp:14
size_t m_numComponents
Number of components calculated for the problem.
Definition: vcs_solve.h:1506
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:173
#define VCS_DATA_PTR(vvv)
Points to the data in a std::vector<> object.
Definition: vcs_internal.h:24
void vcs_reinsert_deleted(size_t kspec)
int m_debug_print_lvl
Debug printing lvl.
Definition: vcs_solve.h:1999
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...
#define VCS_SPECIES_ZEROEDSS
Species is a SS phase, that is currently zeroed out.
Definition: vcs_defs.h:173
#define VCS_ELEM_TYPE_ABSPOS
Normal element constraint consisting of positive coefficients for the formula matrix.
Definition: vcs_defs.h:306
std::vector< double > m_elemAbundances
Element abundances vector.
Definition: vcs_solve.h:1690
int vcs_elcorr(double aa[], double x[])
Definition: vcs_elem.cpp:109
#define VCS_DELETE_MINORSPECIES_CUTOFF
Cutoff relative mole number value, below which species are deleted from the equilibrium problem...
Definition: vcs_defs.h:73
Internal declarations for the VCSnonideal package.
DoubleStarStar m_formulaMatrix
Formula matrix for the problem.
Definition: vcs_solve.h:1538
std::vector< double > m_molNumSpecies_old
Total moles of the species.
Definition: vcs_solve.h:1616
size_t m_numSpeciesTot
Total number of species in the problems.
Definition: vcs_solve.h:1497
int vcsUtil_mlequ(double *c, size_t idem, size_t n, double *b, size_t m)
Invert an n x n matrix and solve m rhs's.
Definition: vcs_util.cpp:297
#define VCS_ELEM_TYPE_CHARGENEUTRALITY
This refers to a charge neutrality of a single phase.
Definition: vcs_defs.h:318
std::vector< std::string > m_elementName
Vector of strings containing the element names.
Definition: vcs_solve.h:1849
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
#define VCS_SPECIES_ACTIVEBUTZERO
Species lies in a multicomponent phase that is active, but species concentration is zero...
Definition: vcs_defs.h:207
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< int > m_speciesStatus
Major -Minor status vector for the species in the problem.
Definition: vcs_solve.h:1833
std::vector< size_t > m_phaseID
Mapping from the species number to the phase number.
Definition: vcs_solve.h:1836
double vcs_tmoles()
Calculates the total number of moles of species in all phases.
void vcs_elabPhase(size_t iphase, double *const elemAbundPhase)
Definition: vcs_elem.cpp:95
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:154
std::vector< char > m_SSPhase
Boolean indicating whether a species belongs to a single-species phase.
Definition: vcs_solve.h:1840
bool vcs_elabcheck(int ibound)
Definition: vcs_elem.cpp:26
#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
#define VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
Unknown refers to the voltage level of a phase.
Definition: vcs_defs.h:368
std::vector< int > m_speciesUnknownType
Specifies the species unknown type.
Definition: vcs_solve.h:1629
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
std::vector< std::string > m_speciesName
Species string name for the kth species.
Definition: vcs_solve.h:1843
size_t m_numSpeciesRdc
Current number of species in the problems.
Definition: vcs_solve.h:1516