Cantera  2.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  */
9 #include "vcs_Exception.h"
10 #include "math.h"
11
12 namespace VCSnonideal
13 {
14
15
16 //! Computes the current elemental abundances vector
17 /*!
18  * Computes the elemental abundances vector, m_elemAbundances[], and stores it
19  * back into the global structure
20  */
22 {
23  for (size_t j = 0; j < m_numElemConstraints; ++j) {
24  m_elemAbundances[j] = 0.0;
25  for (size_t i = 0; i < m_numSpeciesTot; ++i) {
28  }
29  }
30  }
31 }
32
33
34 /*
35  *
36  * vcs_elabcheck:
37  *
38  * This function checks to see if the element abundances are in
39  * compliance. If they are, then TRUE is returned. If not,
40  * FALSE is returned. Note the number of constraints checked is
41  * usually equal to the number of components in the problem. This
42  * routine can check satisfaction of all of the constraints in the
43  * problem, which is equal to ne. However, the solver can't fix
44  * breakage of constraints above nc, because that nc is the
45  * range space by definition. Satisfaction of extra constraints would
46  * have had to occur in the problem specification.
47  *
48  * The constraints should be broken up into 2 sections. If
49  * a constraint involves a formula matrix with positive and
50  * negative signs, and eaSet = 0.0, then you can't expect that the
51  * sum will be zero. There may be roundoff that inhibits this.
52  * However, if the formula matrix is all of one sign, then
53  * this requires that all species with nonzero entries in the
54  * formula matrix be identically zero. We put this into
55  * the logic below.
56  *
57  * Input
58  * -------
59  * ibound = 1 : Checks constraints up to the number of elements
60  * 0 : Checks constraints up to the number of components.
61  *
62  */
63 bool VCS_SOLVE::vcs_elabcheck(int ibound)
64 {
65  size_t top = m_numComponents;
66  double eval, scale;
67  int numNonZero;
68  bool multisign = false;
69  if (ibound) {
71  }
72  /*
73  * Require 12 digits of accuracy on non-zero constraints.
74  */
75  for (size_t i = 0; i < top; ++i) {
76  if (m_elementActive[i]) {
77  if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > (fabs(m_elemAbundancesGoal[i]) * 1.0e-12)) {
78  /*
79  * This logic is for charge neutrality condition
80  */
82  AssertThrowVCS(m_elemAbundancesGoal[i] == 0.0, "vcs_elabcheck");
83  }
86  /*
87  * Find out if the constraint is a multisign constraint.
88  * If it is, then we have to worry about roundoff error
89  * in the addition of terms. We are limited to 13
90  * digits of finite arithmetic accuracy.
91  */
92  numNonZero = 0;
93  multisign = false;
94  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
95  eval = m_formulaMatrix[i][kspec];
96  if (eval < 0.0) {
97  multisign = true;
98  }
99  if (eval != 0.0) {
100  scale = std::max(scale, fabs(eval * m_molNumSpecies_old[kspec]));
101  numNonZero++;
102  }
103  }
104  if (multisign) {
105  if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > 1e-11 * scale) {
106  return false;
107  }
108  } else {
110  return false;
111  }
112  }
113  } else {
114  /*
115  * For normal element balances, we require absolute compliance
116  * even for rediculously small numbers.
117  */
118  if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) {
119  return false;
120  } else {
121  return false;
122  }
123  }
124  }
125  }
126  }
127  return true;
128 } /* vcs_elabcheck() *********************************************************/
129
130 /*****************************************************************************/
131 /*****************************************************************************/
132 /*****************************************************************************/
133
134 void VCS_SOLVE::vcs_elabPhase(size_t iphase, double* const elemAbundPhase)
135
136 /*************************************************************************
137  *
138  * vcs_elabPhase:
139  *
140  * Computes the elemental abundances vector for a single phase,
141  * elemAbundPhase[], and returns it through the argument list.
142  * The mole numbers of species are taken from the current value
143  * in m_molNumSpecies_old[].
144  *************************************************************************/
145 {
146  for (size_t j = 0; j < m_numElemConstraints; ++j) {
147  elemAbundPhase[j] = 0.0;
148  for (size_t i = 0; i < m_numSpeciesTot; ++i) {
150  if (m_phaseID[i] == iphase) {
151  elemAbundPhase[j] += m_formulaMatrix[j][i] * m_molNumSpecies_old[i];
152  }
153  }
154  }
155  }
156 }
157
158 /*****************************************************************************/
159 /*****************************************************************************/
160 /*****************************************************************************/
161
162 int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
163
164 /**************************************************************************
165  *
166  * vcs_elcorr:
167  *
168  * This subroutine corrects for element abundances. At the end of the
169  * surbroutine, the total moles in all phases are recalculated again,
170  * because we have changed the number of moles in this routine.
171  *
172  * Input
173  * -> temporary work vectors:
174  * aa[ne*ne]
175  * x[ne]
176  *
177  * Return Values:
178  * 0 = Nothing of significance happened,
179  * Element abundances were and still are good.
180  * 1 = The solution changed significantly;
181  * The element abundances are now good.
182  * 2 = The solution changed significantly,
183  * The element abundances are still bad.
184  * 3 = The solution changed significantly,
185  * The element abundances are still bad and a component
186  * species got zeroed out.
187  *
188  * Internal data to be worked on::
189  *
190  * ga Current element abundances
191  * m_elemAbundancesGoal Required elemental abundances
192  * m_molNumSpecies_old Current mole number of species.
193  * m_formulaMatrix[][] Formula matrix of the species
194  * ne Number of elements
195  * nc Number of components.
196  *
197  * NOTES:
198  * This routine is turning out to be very problematic. There are
199  * lots of special cases and problems with zeroing out species.
200  *
201  * Still need to check out when we do loops over nc vs. ne.
202  *
203  *************************************************************************/
204 {
205  int retn = 0, its;
206  bool goodSpec;
207  double xx, par, saveDir, dir;
208
209 #ifdef DEBUG_MODE
210  double l2before = 0.0, l2after = 0.0;
211  std::vector<double> ga_save(m_numElemConstraints, 0.0);
212  vcs_dcopy(VCS_DATA_PTR(ga_save), VCS_DATA_PTR(m_elemAbundances), m_numElemConstraints);
213  if (m_debug_print_lvl >= 2) {
214  plogf(" --- vcsc_elcorr: Element abundances correction routine");
215  if (m_numElemConstraints != m_numComponents) {
216  plogf(" (m_numComponents != m_numElemConstraints)");
217  }
218  plogf("\n");
219  }
220
221  for (size_t i = 0; i < m_numElemConstraints; ++i) {
222  x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i];
223  }
224  l2before = 0.0;
225  for (size_t i = 0; i < m_numElemConstraints; ++i) {
226  l2before += x[i] * x[i];
227  }
228  l2before = sqrt(l2before/m_numElemConstraints);
229 #endif
230
231  /*
232  * Special section to take out single species, single component,
233  * moles. These are species which have non-zero entries in the
234  * formula matrix, and no other species have zero values either.
235  *
236  */
237  int numNonZero = 0;
238  bool changed = false;
239  bool multisign = false;
240  for (size_t i = 0; i < m_numElemConstraints; ++i) {
241  numNonZero = 0;
242  multisign = false;
243  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
245  double eval = m_formulaMatrix[i][kspec];
246  if (eval < 0.0) {
247  multisign = true;
248  }
249  if (eval != 0.0) {
250  numNonZero++;
251  }
252  }
253  }
254  if (!multisign) {
255  if (numNonZero < 2) {
256  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
258  double eval = m_formulaMatrix[i][kspec];
259  if (eval > 0.0) {
260  m_molNumSpecies_old[kspec] = m_elemAbundancesGoal[i] / eval;
261  changed = true;
262  }
263  }
264  }
265  } else {
266  int numCompNonZero = 0;
267  size_t compID = npos;
268  for (size_t kspec = 0; kspec < m_numComponents; kspec++) {
270  double eval = m_formulaMatrix[i][kspec];
271  if (eval > 0.0) {
272  compID = kspec;
273  numCompNonZero++;
274  }
275  }
276  }
277  if (numCompNonZero == 1) {
278  double diff = m_elemAbundancesGoal[i];
279  for (size_t kspec = m_numComponents; kspec < m_numSpeciesTot; kspec++) {
281  double eval = m_formulaMatrix[i][kspec];
282  diff -= eval * m_molNumSpecies_old[kspec];
283  }
284  m_molNumSpecies_old[compID] = std::max(0.0,diff/m_formulaMatrix[i][compID]);
285  changed = true;
286  }
287  }
288  }
289  }
290  }
291  if (changed) {
292  vcs_elab();
293  }
294
295  /*
296  * Section to check for maximum bounds errors on all species
297  * due to elements.
298  * This may only be tried on element types which are VCS_ELEM_TYPE_ABSPOS.
299  * This is because no other species may have a negative number of these.
300  *
301  * Note, also we can do this over ne, the number of elements, not just
302  * the number of components.
303  */
304  changed = false;
305  for (size_t i = 0; i < m_numElemConstraints; ++i) {
306  int elType = m_elType[i];
307  if (elType == VCS_ELEM_TYPE_ABSPOS) {
308  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
310  double atomComp = m_formulaMatrix[i][kspec];
311  if (atomComp > 0.0) {
312  double maxPermissible = m_elemAbundancesGoal[i] / atomComp;
313  if (m_molNumSpecies_old[kspec] > maxPermissible) {
314
315 #ifdef DEBUG_MODE
316  if (m_debug_print_lvl >= 3) {
317  plogf(" --- vcs_elcorr: Reduced species %s from %g to %g "
318  "due to %s max bounds constraint\n",
319  m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec],
320  maxPermissible, m_elementName[i].c_str());
321  }
322 #endif
323  m_molNumSpecies_old[kspec] = maxPermissible;
324  changed = true;
326  m_molNumSpecies_old[kspec] = 0.0;
327  if (m_SSPhase[kspec]) {
329  } else {
331  }
332 #ifdef DEBUG_MODE
333  if (m_debug_print_lvl >= 2) {
334  plogf(" --- vcs_elcorr: Zeroed species %s and changed "
335  "status to %d due to max bounds constraint\n",
336  m_speciesName[kspec].c_str(), m_speciesStatus[kspec]);
337  }
338 #endif
339  }
340  }
341  }
342  }
343  }
344  }
345  }
346
347  // Recalculate the element abundances if something has changed.
348  if (changed) {
349  vcs_elab();
350  }
351
352  /*
353  * Ok, do the general case. Linear algebra problem is
354  * of length nc, not ne, as there may be degenerate rows when
355  * nc .ne. ne.
356  */
357  for (size_t i = 0; i < m_numComponents; ++i) {
358  x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i];
359  if (fabs(x[i]) > 1.0E-13) {
360  retn = 1;
361  }
362  for (size_t j = 0; j < m_numComponents; ++j) {
363  aa[j + i*m_numElemConstraints] = m_formulaMatrix[j][i];
364  }
365  }
366  int err = vcsUtil_mlequ(aa, m_numElemConstraints, m_numComponents, x, 1);
367  if (err == 1) {
368  plogf("vcs_elcorr ERROR: mlequ returned error condition\n");
369  return VCS_FAILED_CONVERGENCE;
370  }
371  /*
372  * Now apply the new direction without creating negative species.
373  */
374  par = 0.5;
375  for (size_t i = 0; i < m_numComponents; ++i) {
376  if (m_molNumSpecies_old[i] > 0.0) {
377  xx = -x[i] / m_molNumSpecies_old[i];
378  if (par < xx) {
379  par = xx;
380  }
381  }
382  }
383  if (par > 100.0) {
384  par = 100.0;
385  }
386  par = 1.0 / par;
387  if (par < 1.0 && par > 0.0) {
388  retn = 2;
389  par *= 0.9999;
390  for (size_t i = 0; i < m_numComponents; ++i) {
391  double tmp = m_molNumSpecies_old[i] + par * x[i];
392  if (tmp > 0.0) {
393  m_molNumSpecies_old[i] = tmp;
394  } else {
395  if (m_SSPhase[i]) {
396  m_molNumSpecies_old[i] = 0.0;
397  } else {
398  m_molNumSpecies_old[i] = m_molNumSpecies_old[i] * 0.0001;
399  }
400  }
401  }
402  } else {
403  for (size_t i = 0; i < m_numComponents; ++i) {
404  double tmp = m_molNumSpecies_old[i] + x[i];
405  if (tmp > 0.0) {
406  m_molNumSpecies_old[i] = tmp;
407  } else {
408  if (m_SSPhase[i]) {
409  m_molNumSpecies_old[i] = 0.0;
410  } else {
411  m_molNumSpecies_old[i] = m_molNumSpecies_old[i] * 0.0001;
412  }
413  }
414  }
415  }
416
417  /*
418  * We have changed the element abundances. Calculate them again
419  */
420  vcs_elab();
421  /*
422  * We have changed the total moles in each phase. Calculate them again
423  */
424  vcs_tmoles();
425
426  /*
427  * Try some ad hoc procedures for fixing the problem
428  */
429  if (retn >= 2) {
430  /*
431  * First find a species whose adjustment is a win-win
432  * situation.
433  */
434  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
436  continue;
437  }
438  saveDir = 0.0;
439  goodSpec = true;
440  for (size_t i = 0; i < m_numComponents; ++i) {
441  dir = m_formulaMatrix[i][kspec] * (m_elemAbundancesGoal[i] - m_elemAbundances[i]);
442  if (fabs(dir) > 1.0E-10) {
443  if (dir > 0.0) {
444  if (saveDir < 0.0) {
445  goodSpec = false;
446  break;
447  }
448  } else {
449  if (saveDir > 0.0) {
450  goodSpec = false;
451  break;
452  }
453  }
454  saveDir = dir;
455  } else {
456  if (m_formulaMatrix[i][kspec] != 0.) {
457  goodSpec = false;
458  break;
459  }
460  }
461  }
462  if (goodSpec) {
463  its = 0;
464  xx = 0.0;
465  for (size_t i = 0; i < m_numComponents; ++i) {
466  if (m_formulaMatrix[i][kspec] != 0.0) {
467  xx += (m_elemAbundancesGoal[i] - m_elemAbundances[i]) / m_formulaMatrix[i][kspec];
468  its++;
469  }
470  }
471  if (its > 0) {
472  xx /= its;
473  }
474  m_molNumSpecies_old[kspec] += xx;
475  m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 1.0E-10);
476  /*
477  * If we are dealing with a deleted species, then
478  * we need to reinsert it into the active list.
479  */
480  if (kspec >= m_numSpeciesRdc) {
481  vcs_reinsert_deleted(kspec);
483  vcs_elab();
484  goto L_CLEANUP;
485  }
486  vcs_elab();
487  }
488  }
489  }
490  if (vcs_elabcheck(0)) {
491  retn = 1;
492  goto L_CLEANUP;
493  }
494
495  for (size_t i = 0; i < m_numElemConstraints; ++i) {
497  (m_elType[i] == VCS_ELEM_TYPE_ABSPOS && m_elemAbundancesGoal[i] == 0.0)) {
498  for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
499  if (m_elemAbundances[i] > 0.0) {
500  if (m_formulaMatrix[i][kspec] < 0.0) {
501  m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix[i][kspec] ;
502  if (m_molNumSpecies_old[kspec] < 0.0) {
503  m_molNumSpecies_old[kspec] = 0.0;
504  }
505  vcs_elab();
506  break;
507  }
508  }
509  if (m_elemAbundances[i] < 0.0) {
510  if (m_formulaMatrix[i][kspec] > 0.0) {
511  m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / m_formulaMatrix[i][kspec];
512  if (m_molNumSpecies_old[kspec] < 0.0) {
513  m_molNumSpecies_old[kspec] = 0.0;
514  }
515  vcs_elab();
516  break;
517  }
518  }
519  }
520  }
521  }
522  if (vcs_elabcheck(1)) {
523  retn = 1;
524  goto L_CLEANUP;
525  }
526
527  /*
528  * For electron charges element types, we try positive deltas
529  * in the species concentrations to match the desired
530  * electron charge exactly.
531  */
532  for (size_t i = 0; i < m_numElemConstraints; ++i) {
533  double dev = m_elemAbundancesGoal[i] - m_elemAbundances[i];
534  if (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE && (fabs(dev) > 1.0E-300)) {
535  bool useZeroed = true;
536  for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
537  if (dev < 0.0) {
538  if (m_formulaMatrix[i][kspec] < 0.0) {
539  if (m_molNumSpecies_old[kspec] > 0.0) {
540  useZeroed = false;
541  }
542  }
543  } else {
544  if (m_formulaMatrix[i][kspec] > 0.0) {
545  if (m_molNumSpecies_old[kspec] > 0.0) {
546  useZeroed = false;
547  }
548  }
549  }
550  }
551  for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
552  if (m_molNumSpecies_old[kspec] > 0.0 || useZeroed) {
553  if (dev < 0.0) {
554  if (m_formulaMatrix[i][kspec] < 0.0) {
555  double delta = dev / m_formulaMatrix[i][kspec] ;
556  m_molNumSpecies_old[kspec] += delta;
557  if (m_molNumSpecies_old[kspec] < 0.0) {
558  m_molNumSpecies_old[kspec] = 0.0;
559  }
560  vcs_elab();
561  break;
562  }
563  }
564  if (dev > 0.0) {
565  if (m_formulaMatrix[i][kspec] > 0.0) {
566  double delta = dev / m_formulaMatrix[i][kspec] ;
567  m_molNumSpecies_old[kspec] += delta;
568  if (m_molNumSpecies_old[kspec] < 0.0) {
569  m_molNumSpecies_old[kspec] = 0.0;
570  }
571  vcs_elab();
572  break;
573  }
574  }
575  }
576  }
577  }
578  }
579  if (vcs_elabcheck(1)) {
580  retn = 1;
581  goto L_CLEANUP;
582  }
583
584 L_CLEANUP:
585  ;
586  vcs_tmoles();
587 #ifdef DEBUG_MODE
588  l2after = 0.0;
589  for (int i = 0; i < m_numElemConstraints; ++i) {
590  l2after += SQUARE(m_elemAbundances[i] - m_elemAbundancesGoal[i]);
591  }
592  l2after = sqrt(l2after/m_numElemConstraints);
593  if (m_debug_print_lvl >= 2) {
594  plogf(" --- Elem_Abund: Correct Initial "
595  " Final\n");
596  for (int i = 0; i < m_numElemConstraints; ++i) {
597  plogf(" --- ");
598  plogf("%-2.2s", m_elementName[i].c_str());
599  plogf(" %20.12E %20.12E %20.12E\n", m_elemAbundancesGoal[i], ga_save[i], m_elemAbundances[i]);
600  }
601  plogf(" --- Diff_Norm: %20.12E %20.12E\n",
602  l2before, l2after);
603  }
604 #endif
605  return retn;
606 }
607
608 }
609