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