Cantera  2.4.0
ChemEquil.cpp
Go to the documentation of this file.
1 /**
2  * @file ChemEquil.cpp
3  * Chemical equilibrium. Implementation file for class
4  * ChemEquil.
5  */
6 
7 // This file is part of Cantera. See License.txt in the top-level directory or
8 // at http://www.cantera.org/license.txt for license and copyright information.
9 
13 
14 using namespace std;
15 
16 #include <cstdio>
17 
18 namespace Cantera
19 {
20 
21 int _equilflag(const char* xy)
22 {
23  string flag = string(xy);
24  if (flag == "TP") {
25  return TP;
26  } else if (flag == "TV") {
27  return TV;
28  } else if (flag == "HP") {
29  return HP;
30  } else if (flag == "UV") {
31  return UV;
32  } else if (flag == "SP") {
33  return SP;
34  } else if (flag == "SV") {
35  return SV;
36  } else if (flag == "UP") {
37  return UP;
38  } else {
39  throw CanteraError("_equilflag","unknown property pair "+flag);
40  }
41  return -1;
42 }
43 
44 ChemEquil::ChemEquil() : m_skip(npos), m_elementTotalSum(1.0),
45  m_p0(OneAtm), m_eloc(npos),
46  m_elemFracCutoff(1.0E-100),
47  m_doResPerturb(false)
48 {}
49 
50 ChemEquil::ChemEquil(thermo_t& s) :
51  m_skip(npos),
52  m_elementTotalSum(1.0),
53  m_p0(OneAtm), m_eloc(npos),
54  m_elemFracCutoff(1.0E-100),
55  m_doResPerturb(false)
56 {
57  initialize(s);
58 }
59 
60 ChemEquil::~ChemEquil()
61 {
62 }
63 
65 {
66  // store a pointer to s and some of its properties locally.
67  m_phase = &s;
68  m_p0 = s.refPressure();
69  m_kk = s.nSpecies();
70  m_mm = s.nElements();
72 
73  // allocate space in internal work arrays within the ChemEquil object
74  m_molefractions.resize(m_kk);
75  m_lambda.resize(m_mm, -100.0);
76  m_elementmolefracs.resize(m_mm);
77  m_comp.resize(m_mm * m_kk);
78  m_jwork1.resize(m_mm+2);
79  m_jwork2.resize(m_mm+2);
80  m_startSoln.resize(m_mm+1);
81  m_grt.resize(m_kk);
82  m_mu_RT.resize(m_kk);
83  m_muSS_RT.resize(m_kk);
84  m_component.resize(m_mm,npos);
85  m_orderVectorElements.resize(m_mm);
86 
87  for (size_t m = 0; m < m_mm; m++) {
88  m_orderVectorElements[m] = m;
89  }
90  m_orderVectorSpecies.resize(m_kk);
91  for (size_t k = 0; k < m_kk; k++) {
92  m_orderVectorSpecies[k] = k;
93  }
94 
95  // set up elemental composition matrix
96  size_t mneg = npos;
97  for (size_t m = 0; m < m_mm; m++) {
98  for (size_t k = 0; k < m_kk; k++) {
99  // handle the case of negative atom numbers (used to
100  // represent positive ions, where the 'element' is an
101  // electron
102  if (s.nAtoms(k,m) < 0.0) {
103  // if negative atom numbers have already been specified
104  // for some element other than this one, throw
105  // an exception
106  if (mneg != npos && mneg != m) {
107  throw CanteraError("ChemEquil::initialize",
108  "negative atom numbers allowed for only one element");
109  }
110  mneg = m;
111 
112  // the element should be an electron... if it isn't
113  // print a warning.
114  if (s.atomicWeight(m) > 1.0e-3) {
115  writelog("WARNING: species {} has {} atoms of element {},"
116  " but this element is not an electron.\n",
117  s.speciesName(k), s.nAtoms(k,m), s.elementName(m));
118  }
119  }
120  }
121  }
122  m_eloc = mneg;
123 
124  // set up the elemental composition matrix
125  for (size_t k = 0; k < m_kk; k++) {
126  for (size_t m = 0; m < m_mm; m++) {
127  m_comp[k*m_mm + m] = s.nAtoms(k,m);
128  }
129  }
130 }
131 
133  const vector_fp& lambda_RT, doublereal t)
134 {
135  // Construct the chemical potentials by summing element potentials
136  fill(m_mu_RT.begin(), m_mu_RT.end(), 0.0);
137  for (size_t k = 0; k < m_kk; k++) {
138  for (size_t m = 0; m < m_mm; m++) {
139  m_mu_RT[k] += lambda_RT[m]*nAtoms(k,m);
140  }
141  }
142 
143  // Set the temperature
144  s.setTemperature(t);
145 
146  // Call the phase-specific method to set the phase to the
147  // equilibrium state with the specified species chemical
148  // potentials.
149  s.setToEquilState(m_mu_RT.data());
150  update(s);
151 }
152 
154 {
155  // get the mole fractions, temperature, and density
157  m_temp = s.temperature();
158  m_dens = s.density();
159 
160  // compute the elemental mole fractions
161  double sum = 0.0;
162  for (size_t m = 0; m < m_mm; m++) {
163  m_elementmolefracs[m] = 0.0;
164  for (size_t k = 0; k < m_kk; k++) {
165  m_elementmolefracs[m] += nAtoms(k,m) * m_molefractions[k];
166  if (m_molefractions[k] < 0.0) {
167  throw CanteraError("ChemEquil::update",
168  "negative mole fraction for {}: {}",
169  s.speciesName(k), m_molefractions[k]);
170  }
171  }
172  sum += m_elementmolefracs[m];
173  }
174  // Store the sum for later use
175  m_elementTotalSum = sum;
176  // normalize the element mole fractions
177  for (size_t m = 0; m < m_mm; m++) {
178  m_elementmolefracs[m] /= sum;
179  }
180 }
181 
183  int loglevel)
184 {
185  MultiPhase mp;
186  mp.addPhase(&s, 1.0);
187  mp.init();
188  MultiPhaseEquil e(&mp, true, loglevel-1);
189  e.setInitialMixMoles(loglevel-1);
190 
191  // store component indices
192  m_nComponents = std::min(m_nComponents, m_kk);
193  for (size_t m = 0; m < m_nComponents; m++) {
194  m_component[m] = e.componentIndex(m);
195  }
196 
197  // Update the current values of the temp, density, and mole fraction,
198  // and element abundance vectors kept within the ChemEquil object.
199  update(s);
200 
201  if (m_loglevel > 0) {
202  writelog("setInitialMoles: Estimated Mole Fractions\n");
203  writelogf(" Temperature = %g\n", s.temperature());
204  writelogf(" Pressure = %g\n", s.pressure());
205  for (size_t k = 0; k < m_kk; k++) {
206  writelogf(" %-12s % -10.5g\n",
207  s.speciesName(k), s.moleFraction(k));
208  }
209  writelog(" Element_Name ElementGoal ElementMF\n");
210  for (size_t m = 0; m < m_mm; m++) {
211  writelogf(" %-12s % -10.5g% -10.5g\n",
212  s.elementName(m), elMoleGoal[m], m_elementmolefracs[m]);
213  }
214  }
215  return 0;
216 }
217 
219  vector_fp& elMolesGoal, int loglevel)
220 {
221  vector_fp b(m_mm, -999.0);
222  vector_fp mu_RT(m_kk, 0.0);
223  vector_fp xMF_est(m_kk, 0.0);
224 
225  s.getMoleFractions(xMF_est.data());
226  for (size_t n = 0; n < s.nSpecies(); n++) {
227  xMF_est[n] = std::max(xMF_est[n], 1e-20);
228  }
229  s.setMoleFractions(xMF_est.data());
230  s.getMoleFractions(xMF_est.data());
231 
232  MultiPhase mp;
233  mp.addPhase(&s, 1.0);
234  mp.init();
235  int usedZeroedSpecies = 0;
236  vector_fp formRxnMatrix;
237  m_nComponents = BasisOptimize(&usedZeroedSpecies, false,
238  &mp, m_orderVectorSpecies,
239  m_orderVectorElements, formRxnMatrix);
240 
241  for (size_t m = 0; m < m_nComponents; m++) {
242  size_t k = m_orderVectorSpecies[m];
243  m_component[m] = k;
244  xMF_est[k] = std::max(xMF_est[k], 1e-8);
245  }
246  s.setMoleFractions(xMF_est.data());
247  s.getMoleFractions(xMF_est.data());
248 
249  ElemRearrange(m_nComponents, elMolesGoal, &mp,
250  m_orderVectorSpecies, m_orderVectorElements);
251 
252  s.getChemPotentials(mu_RT.data());
253  scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(),
254  1.0/(GasConstant* s.temperature()));
255 
256  if (loglevel > 0) {
257  for (size_t m = 0; m < m_nComponents; m++) {
258  size_t isp = m_component[m];
259  writelogf("isp = %d, %s\n", isp, s.speciesName(isp));
260  }
261  writelogf("Pressure = %g\n", s.pressure());
262  writelogf("Temperature = %g\n", s.temperature());
263  writelog(" id Name MF mu/RT \n");
264  for (size_t n = 0; n < s.nSpecies(); n++) {
265  writelogf("%10d %15s %10.5g %10.5g\n",
266  n, s.speciesName(n), xMF_est[n], mu_RT[n]);
267  }
268  }
270  for (size_t m = 0; m < m_nComponents; m++) {
271  for (size_t n = 0; n < m_nComponents; n++) {
272  aa(m,n) = nAtoms(m_component[m], m_orderVectorElements[n]);
273  }
274  b[m] = mu_RT[m_component[m]];
275  }
276 
277  int info = solve(aa, b.data());
278  if (info) {
279  info = -2;
280  }
281  for (size_t m = 0; m < m_nComponents; m++) {
282  lambda_RT[m_orderVectorElements[m]] = b[m];
283  }
284  for (size_t m = m_nComponents; m < m_mm; m++) {
285  lambda_RT[m_orderVectorElements[m]] = 0.0;
286  }
287 
288  if (loglevel > 0) {
289  writelog(" id CompSpecies ChemPot EstChemPot Diff\n");
290  for (size_t m = 0; m < m_nComponents; m++) {
291  size_t isp = m_component[m];
292  double tmp = 0.0;
293  for (size_t n = 0; n < m_mm; n++) {
294  tmp += nAtoms(isp, n) * lambda_RT[n];
295  }
296  writelogf("%3d %16s %10.5g %10.5g %10.5g\n",
297  m, s.speciesName(isp), mu_RT[isp], tmp, tmp - mu_RT[isp]);
298  }
299 
300  writelog(" id ElName Lambda_RT\n");
301  for (size_t m = 0; m < m_mm; m++) {
302  writelogf(" %3d %6s %10.5g\n", m, s.elementName(m), lambda_RT[m]);
303  }
304  }
305  return info;
306 }
307 
308 int ChemEquil::equilibrate(thermo_t& s, const char* XY, int loglevel)
309 {
310  initialize(s);
311  update(s);
312  vector_fp elMolesGoal = m_elementmolefracs;
313  return equilibrate(s, XY, elMolesGoal, loglevel-1);
314 }
315 
316 int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
317  vector_fp& elMolesGoal, int loglevel)
318 {
319  int fail = 0;
320  bool tempFixed = true;
321  int XY = _equilflag(XYstr);
322  vector_fp state;
323  s.saveState(state);
324  m_loglevel = loglevel;
325 
326  // Check Compatibility
327  if (m_mm != s.nElements() || m_kk != s.nSpecies()) {
328  throw CanteraError("ChemEquil::equilibrate ERROR",
329  "Input ThermoPhase is incompatible with initialization");
330  }
331 
332  initialize(s);
333  update(s);
334  switch (XY) {
335  case TP:
336  case PT:
337  m_p1 = [](ThermoPhase& s) { return s.temperature(); };
338  m_p2 = [](ThermoPhase& s) { return s.pressure(); };
339  break;
340  case HP:
341  case PH:
342  tempFixed = false;
343  m_p1 = [](ThermoPhase& s) { return s.enthalpy_mass(); };
344  m_p2 = [](ThermoPhase& s) { return s.pressure(); };
345  break;
346  case SP:
347  case PS:
348  tempFixed = false;
349  m_p1 = [](ThermoPhase& s) { return s.entropy_mass(); };
350  m_p2 = [](ThermoPhase& s) { return s.pressure(); };
351  break;
352  case SV:
353  case VS:
354  tempFixed = false;
355  m_p1 = [](ThermoPhase& s) { return s.entropy_mass(); };
356  m_p2 = [](ThermoPhase& s) { return s.density(); };
357  break;
358  case TV:
359  case VT:
360  m_p1 = [](ThermoPhase& s) { return s.temperature(); };
361  m_p2 = [](ThermoPhase& s) { return s.density(); };
362  break;
363  case UV:
364  case VU:
365  tempFixed = false;
366  m_p1 = [](ThermoPhase& s) { return s.intEnergy_mass(); };
367  m_p2 = [](ThermoPhase& s) { return s.density(); };
368  break;
369  default:
370  throw CanteraError("equilibrate","illegal property pair.");
371  }
372  // If the temperature is one of the specified variables, and
373  // it is outside the valid range, throw an exception.
374  if (tempFixed) {
375  double tfixed = s.temperature();
376  if (tfixed > s.maxTemp() + 1.0 || tfixed < s.minTemp() - 1.0) {
377  throw CanteraError("ChemEquil::equilibrate", "Specified temperature"
378  " ({} K) outside valid range of {} K to {} K\n",
379  s.temperature(), s.minTemp(), s.maxTemp());
380  }
381  }
382 
383  // Before we do anything to change the ThermoPhase object, we calculate and
384  // store the two specified thermodynamic properties that we are after.
385  double xval = m_p1(s);
386  double yval = m_p2(s);
387 
388  size_t mm = m_mm;
389  size_t nvar = mm + 1;
390  DenseMatrix jac(nvar, nvar); // Jacobian
391  vector_fp x(nvar, -102.0); // solution vector
392  vector_fp res_trial(nvar, 0.0); // residual
393 
394  // Replace one of the element abundance fraction equations with the
395  // specified property calculation.
396  //
397  // We choose the equation of the element with the highest element abundance.
398  double tmp = -1.0;
399  for (size_t im = 0; im < m_nComponents; im++) {
400  size_t m = m_orderVectorElements[im];
401  if (elMolesGoal[m] > tmp) {
402  m_skip = m;
403  tmp = elMolesGoal[m];
404  }
405  }
406  if (tmp <= 0.0) {
407  throw CanteraError("ChemEquil",
408  "Element Abundance Vector is zeroed");
409  }
410 
411  // start with a composition with everything non-zero. Note that since we
412  // have already save the target element moles, changing the composition at
413  // this point only affects the starting point, not the final solution.
414  vector_fp xmm(m_kk, 0.0);
415  for (size_t k = 0; k < m_kk; k++) {
416  xmm[k] = s.moleFraction(k) + 1.0E-32;
417  }
418  s.setMoleFractions(xmm.data());
419 
420  // Update the internally stored values of m_temp, m_dens, and the element
421  // mole fractions.
422  update(s);
423 
424  doublereal tmaxPhase = s.maxTemp();
425  doublereal tminPhase = s.minTemp();
426  // loop to estimate T
427  if (!tempFixed) {
428  doublereal tmin = std::max(s.temperature(), tminPhase);
429  if (tmin > tmaxPhase) {
430  tmin = tmaxPhase - 20;
431  }
432  doublereal tmax = std::min(tmin + 10., tmaxPhase);
433  if (tmax < tminPhase) {
434  tmax = tminPhase + 20;
435  }
436 
437  doublereal slope, phigh, plow, pval, dt;
438 
439  // first get the property values at the upper and lower temperature
440  // limits. Since p1 (h, s, or u) is monotonic in T, these values
441  // determine the upper and lower bounds (phigh, plow) for p1.
442 
443  s.setTemperature(tmax);
444  setInitialMoles(s, elMolesGoal, loglevel - 1);
445  phigh = m_p1(s);
446 
447  s.setTemperature(tmin);
448  setInitialMoles(s, elMolesGoal, loglevel - 1);
449  plow = m_p1(s);
450 
451  // start with T at the midpoint of the range
452  doublereal t0 = 0.5*(tmin + tmax);
453  s.setTemperature(t0);
454 
455  // loop up to 5 times
456  for (int it = 0; it < 10; it++) {
457  // set the composition and get p1
458  setInitialMoles(s, elMolesGoal, loglevel - 1);
459  pval = m_p1(s);
460 
461  // If this value of p1 is greater than the specified property value,
462  // then the current temperature is too high. Use it as the new upper
463  // bound. Otherwise, it is too low, so use it as the new lower
464  // bound.
465  if (pval > xval) {
466  tmax = t0;
467  phigh = pval;
468  } else {
469  tmin = t0;
470  plow = pval;
471  }
472 
473  // Determine the new T estimate by linearly interpolating
474  // between the upper and lower bounds
475  slope = (phigh - plow)/(tmax - tmin);
476  dt = (xval - pval)/slope;
477 
478  // If within 50 K, terminate the search
479  if (fabs(dt) < 50.0) {
480  break;
481  }
482  dt = clip(dt, -200.0, 200.0);
483  if ((t0 + dt) < tminPhase) {
484  dt = 0.5*((t0) + tminPhase) - t0;
485  }
486  if ((t0 + dt) > tmaxPhase) {
487  dt = 0.5*((t0) + tmaxPhase) - t0;
488  }
489  // update the T estimate
490  t0 += dt;
491  if (t0 <= tminPhase || t0 >= tmaxPhase || t0 < 100.0) {
492  throw CanteraError("ChemEquil::equilibrate", "T out of bounds");
493  }
494  s.setTemperature(t0);
495  }
496  }
497 
498  setInitialMoles(s, elMolesGoal,loglevel);
499 
500  // Calculate initial estimates of the element potentials. This algorithm
501  // uses the MultiPhaseEquil object's initialization capabilities to
502  // calculate an initial estimate of the mole fractions for a set of linearly
503  // independent component species. Then, the element potentials are solved
504  // for based on the chemical potentials of the component species.
505  estimateElementPotentials(s, x, elMolesGoal);
506 
507  // Do a better estimate of the element potentials. We have found that the
508  // current estimate may not be good enough to avoid drastic numerical issues
509  // associated with the use of a numerically generated Jacobian.
510  //
511  // The Brinkley algorithm assumes a constant T, P system and uses a
512  // linearized analytical Jacobian that turns out to be very stable.
513  int info = estimateEP_Brinkley(s, x, elMolesGoal);
514  if (info == 0) {
515  setToEquilState(s, x, s.temperature());
516  }
517 
518  // Install the log(temp) into the last solution unknown slot.
519  x[m_mm] = log(s.temperature());
520 
521  // Setting the max and min values for x[]. Also, if element abundance vector
522  // is zero, setting x[] to -1000. This effectively zeroes out all species
523  // containing that element.
524  vector_fp above(nvar);
525  vector_fp below(nvar);
526  for (size_t m = 0; m < mm; m++) {
527  above[m] = 200.0;
528  below[m] = -2000.0;
529  if (elMolesGoal[m] < m_elemFracCutoff && m != m_eloc) {
530  x[m] = -1000.0;
531  }
532  }
533 
534  // Set the temperature bounds to be 25 degrees different than the max and
535  // min temperatures.
536  above[mm] = log(s.maxTemp() + 25.0);
537  below[mm] = log(s.minTemp() - 25.0);
538 
539  vector_fp grad(nvar, 0.0); // gradient of f = F*F/2
540  vector_fp oldx(nvar, 0.0); // old solution
541  vector_fp oldresid(nvar, 0.0);
542 
543  for (int iter = 0; iter < options.maxIterations; iter++) {
544  // check for convergence.
545  equilResidual(s, x, elMolesGoal, res_trial, xval, yval);
546  double f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin());
547  double xx = m_p1(s);
548  double yy = m_p2(s);
549  double deltax = (xx - xval)/xval;
550  double deltay = (yy - yval)/yval;
551  bool passThis = true;
552  for (size_t m = 0; m < nvar; m++) {
553  double tval = options.relTolerance;
554  if (m < mm) {
555  // Special case convergence requirements for electron element.
556  // This is a special case because the element coefficients may
557  // be both positive and negative. And, typically they sum to
558  // 0.0. Therefore, there is no natural absolute value for this
559  // quantity. We supply the absolute value tolerance here. Note,
560  // this is made easier since the element abundances are
561  // normalized to one within this routine.
562  //
563  // Note, the 1.0E-13 value was recently relaxed from 1.0E-15,
564  // because convergence failures were found to occur for the
565  // lower value at small pressure (0.01 pascal).
566  if (m == m_eloc) {
567  tval = elMolesGoal[m] * options.relTolerance + options.absElemTol
568  + 1.0E-13;
569  } else {
570  tval = elMolesGoal[m] * options.relTolerance + options.absElemTol;
571  }
572  }
573  if (fabs(res_trial[m]) > tval) {
574  passThis = false;
575  }
576  }
577  if (iter > 0 && passThis && fabs(deltax) < options.relTolerance
578  && fabs(deltay) < options.relTolerance) {
579  options.iterations = iter;
580  for (size_t m = 0; m < m_mm; m++) {
581  m_lambda[m] = x[m]* s.RT();
582  }
583 
584  if (m_eloc != npos) {
585  adjustEloc(s, elMolesGoal);
586  }
587 
588  if (s.temperature() > s.maxTemp() + 1.0 ||
589  s.temperature() < s.minTemp() - 1.0) {
590  writelog("Warning: Temperature ({} K) outside valid range of "
591  "{} K to {} K\n",
592  s.temperature(), s.minTemp(), s.maxTemp());
593  }
594  return 0;
595  }
596  // compute the residual and the Jacobian using the current
597  // solution vector
598  equilResidual(s, x, elMolesGoal, res_trial, xval, yval);
599  f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin());
600 
601  // Compute the Jacobian matrix
602  equilJacobian(s, x, elMolesGoal, jac, xval, yval);
603 
604  if (m_loglevel > 0) {
605  writelogf("Jacobian matrix %d:\n", iter);
606  for (size_t m = 0; m <= m_mm; m++) {
607  writelog(" [ ");
608  for (size_t n = 0; n <= m_mm; n++) {
609  writelog("{:10.5g} ", jac(m,n));
610  }
611  writelog(" ]");
612  if (m < m_mm) {
613  writelog("x_{:10s}", s.elementName(m));
614  } else if (m_eloc == m) {
615  writelog("x_ELOC");
616  } else if (m == m_skip) {
617  writelog("x_YY");
618  } else {
619  writelog("x_XX");
620  }
621  writelog(" = - ({:10.5g})\n", res_trial[m]);
622  }
623  }
624 
625  oldx = x;
626  double oldf = f;
627  scale(res_trial.begin(), res_trial.end(), res_trial.begin(), -1.0);
628 
629  // Solve the system
630  try {
631  info = solve(jac, res_trial.data());
632  } catch (CanteraError& err) {
633  s.restoreState(state);
634  throw CanteraError("equilibrate",
635  "Jacobian is singular. \nTry adding more species, "
636  "changing the elemental composition slightly, \nor removing "
637  "unused elements.\n\n" + err.getMessage());
638  }
639 
640  // find the factor by which the Newton step can be multiplied
641  // to keep the solution within bounds.
642  double fctr = 1.0;
643  for (size_t m = 0; m < nvar; m++) {
644  double newval = x[m] + res_trial[m];
645  if (newval > above[m]) {
646  fctr = std::max(0.0,
647  std::min(fctr,0.8*(above[m] - x[m])/(newval - x[m])));
648  } else if (newval < below[m]) {
649  if (m < m_mm && (m != m_skip)) {
650  res_trial[m] = -50;
651  if (x[m] < below[m] + 50.) {
652  res_trial[m] = below[m] - x[m];
653  }
654  } else {
655  fctr = std::min(fctr, 0.8*(x[m] - below[m])/(x[m] - newval));
656  }
657  }
658  // Delta Damping
659  if (m == mm && fabs(res_trial[mm]) > 0.2) {
660  fctr = std::min(fctr, 0.2/fabs(res_trial[mm]));
661  }
662  }
663  if (fctr != 1.0 && loglevel > 0) {
664  writelogf("WARNING Soln Damping because of bounds: %g\n", fctr);
665  }
666 
667  // multiply the step by the scaling factor
668  scale(res_trial.begin(), res_trial.end(), res_trial.begin(), fctr);
669 
670  if (!dampStep(s, oldx, oldf, grad, res_trial,
671  x, f, elMolesGoal , xval, yval)) {
672  fail++;
673  if (fail > 3) {
674  s.restoreState(state);
675  throw CanteraError("equilibrate",
676  "Cannot find an acceptable Newton damping coefficient.");
677  }
678  } else {
679  fail = 0;
680  }
681  }
682 
683  // no convergence
684  s.restoreState(state);
685  throw CanteraError("ChemEquil::equilibrate",
686  "no convergence in {} iterations.", options.maxIterations);
687 }
688 
689 
691  double oldf, vector_fp& grad, vector_fp& step, vector_fp& x,
692  double& f, vector_fp& elmols, double xval, double yval)
693 {
694  // Carry out a delta damping approach on the dimensionless element
695  // potentials.
696  double damp = 1.0;
697  for (size_t m = 0; m < m_mm; m++) {
698  if (m == m_eloc) {
699  if (step[m] > 1.25) {
700  damp = std::min(damp, 1.25 /step[m]);
701  }
702  if (step[m] < -1.25) {
703  damp = std::min(damp, -1.25 / step[m]);
704  }
705  } else {
706  if (step[m] > 0.75) {
707  damp = std::min(damp, 0.75 /step[m]);
708  }
709  if (step[m] < -0.75) {
710  damp = std::min(damp, -0.75 / step[m]);
711  }
712  }
713  }
714 
715  // Update the solution unknown
716  for (size_t m = 0; m < x.size(); m++) {
717  x[m] = oldx[m] + damp * step[m];
718  }
719  if (m_loglevel > 0) {
720  writelogf("Solution Unknowns: damp = %g\n", damp);
721  writelog(" X_new X_old Step\n");
722  for (size_t m = 0; m < m_mm; m++) {
723  writelogf(" % -10.5g % -10.5g % -10.5g\n", x[m], oldx[m], step[m]);
724  }
725  }
726  return 1;
727 }
728 
730  const vector_fp& elmFracGoal, vector_fp& resid,
731  doublereal xval, doublereal yval, int loglevel)
732 {
733  setToEquilState(s, x, exp(x[m_mm]));
734 
735  // residuals are the total element moles
736  vector_fp& elmFrac = m_elementmolefracs;
737  for (size_t n = 0; n < m_mm; n++) {
738  size_t m = m_orderVectorElements[n];
739  // drive element potential for absent elements to -1000
740  if (elmFracGoal[m] < m_elemFracCutoff && m != m_eloc) {
741  resid[m] = x[m] + 1000.0;
742  } else if (n >= m_nComponents) {
743  resid[m] = x[m];
744  } else {
745  // Change the calculation for small element number, using
746  // L'Hopital's rule. The log formulation is unstable.
747  if (elmFracGoal[m] < 1.0E-10 || elmFrac[m] < 1.0E-10 || m == m_eloc) {
748  resid[m] = elmFracGoal[m] - elmFrac[m];
749  } else {
750  resid[m] = log((1.0 + elmFracGoal[m]) / (1.0 + elmFrac[m]));
751  }
752  }
753  }
754 
755  if (loglevel > 0 && !m_doResPerturb) {
756  writelog("Residual: ElFracGoal ElFracCurrent Resid\n");
757  for (size_t n = 0; n < m_mm; n++) {
758  writelogf(" % -14.7E % -14.7E % -10.5E\n",
759  elmFracGoal[n], elmFrac[n], resid[n]);
760  }
761  }
762 
763  double xx = m_p1(s);
764  double yy = m_p2(s);
765  resid[m_mm] = xx/xval - 1.0;
766  resid[m_skip] = yy/yval - 1.0;
767 
768  if (loglevel > 0 && !m_doResPerturb) {
769  writelog(" Goal Xvalue Resid\n");
770  writelogf(" XX : % -14.7E % -14.7E % -10.5E\n", xval, xx, resid[m_mm]);
771  writelogf(" YY(%1d): % -14.7E % -14.7E % -10.5E\n", m_skip, yval, yy, resid[m_skip]);
772  }
773 }
774 
775 void ChemEquil::equilJacobian(thermo_t& s, vector_fp& x,
776  const vector_fp& elmols, DenseMatrix& jac,
777  doublereal xval, doublereal yval, int loglevel)
778 {
779  vector_fp& r0 = m_jwork1;
780  vector_fp& r1 = m_jwork2;
781  size_t len = x.size();
782  r0.resize(len);
783  r1.resize(len);
784  doublereal atol = 1.e-10;
785 
786  equilResidual(s, x, elmols, r0, xval, yval, loglevel-1);
787 
788  m_doResPerturb = false;
789  for (size_t n = 0; n < len; n++) {
790  double xsave = x[n];
791  double dx = std::max(atol, fabs(xsave) * 1.0E-7);
792  x[n] = xsave + dx;
793  dx = x[n] - xsave;
794  double rdx = 1.0/dx;
795 
796  // calculate perturbed residual
797  equilResidual(s, x, elmols, r1, xval, yval, loglevel-1);
798 
799  // compute nth column of Jacobian
800  for (size_t m = 0; m < x.size(); m++) {
801  jac(m, n) = (r1[m] - r0[m])*rdx;
802  }
803  x[n] = xsave;
804  }
805  m_doResPerturb = false;
806 }
807 
808 double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double& n_t,
809  const vector_fp& Xmol_i_calc,
810  vector_fp& eMolesCalc, vector_fp& n_i_calc,
811  double pressureConst)
812 {
813  double n_t_calc = 0.0;
814 
815  // Calculate the activity coefficients of the solution, at the previous
816  // solution state.
817  vector_fp actCoeff(m_kk, 1.0);
818  s.setMoleFractions(Xmol_i_calc.data());
819  s.setPressure(pressureConst);
820  s.getActivityCoefficients(actCoeff.data());
821 
822  for (size_t k = 0; k < m_kk; k++) {
823  double tmp = - (m_muSS_RT[k] + log(actCoeff[k]));
824  for (size_t m = 0; m < m_mm; m++) {
825  tmp += nAtoms(k,m) * x[m];
826  }
827  tmp = std::min(tmp, 100.0);
828  if (tmp < -300.) {
829  n_i_calc[k] = 0.0;
830  } else {
831  n_i_calc[k] = n_t * exp(tmp);
832  }
833  n_t_calc += n_i_calc[k];
834  }
835  for (size_t m = 0; m < m_mm; m++) {
836  eMolesCalc[m] = 0.0;
837  for (size_t k = 0; k < m_kk; k++) {
838  eMolesCalc[m] += nAtoms(k,m) * n_i_calc[k];
839  }
840  }
841  return n_t_calc;
842 }
843 
845  vector_fp& elMoles)
846 {
847  // Before we do anything, we will save the state of the solution. Then, if
848  // things go drastically wrong, we will restore the saved state.
849  vector_fp state;
850  s.saveState(state);
851  bool modifiedMatrix = false;
852  size_t neq = m_mm+1;
853  int retn = 1;
854  DenseMatrix a1(neq, neq, 0.0);
855  vector_fp b(neq, 0.0);
856  vector_fp n_i(m_kk,0.0);
857  vector_fp n_i_calc(m_kk,0.0);
858  vector_fp actCoeff(m_kk, 1.0);
859  double beta = 1.0;
860 
861  s.getMoleFractions(n_i.data());
862  double pressureConst = s.pressure();
863  vector_fp Xmol_i_calc = n_i;
864 
865  vector_fp x_old(m_mm+1, 0.0);
866  vector_fp resid(m_mm+1, 0.0);
867  vector_int lumpSum(m_mm+1, 0);
868 
869  // Get the nondimensional Gibbs functions for the species at their standard
870  // states of solution at the current T and P of the solution.
871  s.getGibbs_RT(m_muSS_RT.data());
872 
873  vector_fp eMolesCalc(m_mm, 0.0);
874  vector_fp eMolesFix(m_mm, 0.0);
875  double elMolesTotal = 0.0;
876  for (size_t m = 0; m < m_mm; m++) {
877  elMolesTotal += elMoles[m];
878  for (size_t k = 0; k < m_kk; k++) {
879  eMolesFix[m] += nAtoms(k,m) * n_i[k];
880  }
881  }
882 
883  for (size_t m = 0; m < m_mm; m++) {
884  if (elMoles[m] > 1.0E-70) {
885  x[m] = clip(x[m], -100.0, 50.0);
886  } else {
887  x[m] = clip(x[m], -1000.0, 50.0);
888  }
889  }
890 
891  double n_t = 0.0;
892  double nAtomsMax = 1.0;
893  s.setMoleFractions(Xmol_i_calc.data());
894  s.setPressure(pressureConst);
895  s.getActivityCoefficients(actCoeff.data());
896  for (size_t k = 0; k < m_kk; k++) {
897  double tmp = - (m_muSS_RT[k] + log(actCoeff[k]));
898  double sum2 = 0.0;
899  for (size_t m = 0; m < m_mm; m++) {
900  double sum = nAtoms(k,m);
901  tmp += sum * x[m];
902  sum2 += sum;
903  nAtomsMax = std::max(nAtomsMax, sum2);
904  }
905  if (tmp > 100.) {
906  n_t += 2.8E43;
907  } else {
908  n_t += exp(tmp);
909  }
910  }
911 
912  if (m_loglevel > 0) {
913  writelog("estimateEP_Brinkley::\n\n");
914  writelogf("temp = %g\n", s.temperature());
915  writelogf("pres = %g\n", s.pressure());
916  writelog("Initial mole numbers and mu_SS:\n");
917  writelog(" Name MoleNum mu_SS actCoeff\n");
918  for (size_t k = 0; k < m_kk; k++) {
919  writelogf("%15s %13.5g %13.5g %13.5g\n",
920  s.speciesName(k), n_i[k], m_muSS_RT[k], actCoeff[k]);
921  }
922  writelogf("Initial n_t = %10.5g\n", n_t);
923  writelog("Comparison of Goal Element Abundance with Initial Guess:\n");
924  writelog(" eName eCurrent eGoal\n");
925  for (size_t m = 0; m < m_mm; m++) {
926  writelogf("%5s %13.5g %13.5g\n",
927  s.elementName(m), eMolesFix[m], elMoles[m]);
928  }
929  }
930  for (size_t m = 0; m < m_mm; m++) {
931  if (m != m_eloc && elMoles[m] <= options.absElemTol) {
932  x[m] = -200.;
933  }
934  }
935 
936  // Main Loop.
937  for (int iter = 0; iter < 20* options.maxIterations; iter++) {
938  // Save the old solution
939  for (size_t m = 0; m < m_mm; m++) {
940  x_old[m] = x[m];
941  }
942  x_old[m_mm] = n_t;
943  // Calculate the mole numbers of species
944  if (m_loglevel > 0) {
945  writelogf("START ITERATION %d:\n", iter);
946  }
947  // Calculate the mole numbers of species and elements.
948  double n_t_calc = calcEmoles(s, x, n_t, Xmol_i_calc, eMolesCalc, n_i_calc,
949  pressureConst);
950 
951  for (size_t k = 0; k < m_kk; k++) {
952  Xmol_i_calc[k] = n_i_calc[k]/n_t_calc;
953  }
954 
955  if (m_loglevel > 0) {
956  writelog(" Species: Calculated_Moles Calculated_Mole_Fraction\n");
957  for (size_t k = 0; k < m_kk; k++) {
958  writelogf("%15s: %10.5g %10.5g\n",
959  s.speciesName(k), n_i_calc[k], Xmol_i_calc[k]);
960  }
961  writelogf("%15s: %10.5g\n", "Total Molar Sum", n_t_calc);
962  writelogf("(iter %d) element moles bal: Goal Calculated\n", iter);
963  for (size_t m = 0; m < m_mm; m++) {
964  writelogf(" %8s: %10.5g %10.5g \n",
965  s.elementName(m), elMoles[m], eMolesCalc[m]);
966  }
967  }
968 
969  bool normalStep = true;
970  // Decide if we are to do a normal step or a modified step
971  size_t iM = npos;
972  for (size_t m = 0; m < m_mm; m++) {
973  if (elMoles[m] > 0.001 * elMolesTotal) {
974  if (eMolesCalc[m] > 1000. * elMoles[m]) {
975  normalStep = false;
976  iM = m;
977  }
978  if (1000 * eMolesCalc[m] < elMoles[m]) {
979  normalStep = false;
980  iM = m;
981  }
982  }
983  }
984  if (m_loglevel > 0 && !normalStep) {
985  writelogf(" NOTE: iter(%d) Doing an abnormal step due to row %d\n", iter, iM);
986  }
987  if (!normalStep) {
988  beta = 1.0;
989  resid[m_mm] = 0.0;
990  for (size_t im = 0; im < m_mm; im++) {
991  size_t m = m_orderVectorElements[im];
992  resid[m] = 0.0;
993  if (im < m_nComponents && elMoles[m] > 0.001 * elMolesTotal) {
994  if (eMolesCalc[m] > 1000. * elMoles[m]) {
995  resid[m] = -0.5;
996  resid[m_mm] -= 0.5;
997  }
998  if (1000 * eMolesCalc[m] < elMoles[m]) {
999  resid[m] = 0.5;
1000  resid[m_mm] += 0.5;
1001  }
1002  }
1003  }
1004  if (n_t < (elMolesTotal / nAtomsMax)) {
1005  if (resid[m_mm] < 0.0) {
1006  resid[m_mm] = 0.1;
1007  }
1008  } else if (n_t > elMolesTotal) {
1009  resid[m_mm] = std::min(resid[m_mm], 0.0);
1010  }
1011  } else {
1012  // Determine whether the matrix should be dumbed down because the
1013  // coefficient matrix of species (with significant concentrations)
1014  // is rank deficient.
1015  //
1016  // The basic idea is that at any time during the calculation only a
1017  // small subset of species with sufficient concentration matters. If
1018  // the rank of the element coefficient matrix for that subset of
1019  // species is less than the number of elements, then the matrix
1020  // created by the Brinkley method below may become singular.
1021  //
1022  // The logic below looks for obvious cases where the current element
1023  // coefficient matrix is rank deficient.
1024  //
1025  // The way around rank-deficiency is to lump-sum the corresponding
1026  // row of the matrix. Note, lump-summing seems to work very well in
1027  // terms of its stability properties, i.e., it heads in the right
1028  // direction, albeit with lousy convergence rates.
1029  //
1030  // NOTE: This probably should be extended to a full blown Gauss-
1031  // Jordan factorization scheme in the future. For Example the scheme
1032  // below would fail for the set: HCl NH4Cl, NH3. Hopefully, it's
1033  // caught by the equal rows logic below.
1034  for (size_t m = 0; m < m_mm; m++) {
1035  lumpSum[m] = 1;
1036  }
1037 
1038  double nCutoff = 1.0E-9 * n_t_calc;
1039  if (m_loglevel > 0) {
1040  writelog(" Lump Sum Elements Calculation: \n");
1041  }
1042  for (size_t m = 0; m < m_mm; m++) {
1043  size_t kMSp = npos;
1044  size_t kMSp2 = npos;
1045  int nSpeciesWithElem = 0;
1046  for (size_t k = 0; k < m_kk; k++) {
1047  if (n_i_calc[k] > nCutoff && fabs(nAtoms(k,m)) > 0.001) {
1048  nSpeciesWithElem++;
1049  if (kMSp != npos) {
1050  kMSp2 = k;
1051  double factor = fabs(nAtoms(kMSp,m) / nAtoms(kMSp2,m));
1052  for (size_t n = 0; n < m_mm; n++) {
1053  if (fabs(factor * nAtoms(kMSp2,n) - nAtoms(kMSp,n)) > 1.0E-8) {
1054  lumpSum[m] = 0;
1055  break;
1056  }
1057  }
1058  } else {
1059  kMSp = k;
1060  }
1061  }
1062  }
1063  if (m_loglevel > 0) {
1064  writelogf(" %5s %3d : %5d %5d\n",
1065  s.elementName(m), lumpSum[m], kMSp, kMSp2);
1066  }
1067  }
1068 
1069  // Formulate the matrix.
1070  for (size_t im = 0; im < m_mm; im++) {
1071  size_t m = m_orderVectorElements[im];
1072  if (im < m_nComponents) {
1073  for (size_t n = 0; n < m_mm; n++) {
1074  a1(m,n) = 0.0;
1075  for (size_t k = 0; k < m_kk; k++) {
1076  a1(m,n) += nAtoms(k,m) * nAtoms(k,n) * n_i_calc[k];
1077  }
1078  }
1079  a1(m,m_mm) = eMolesCalc[m];
1080  a1(m_mm, m) = eMolesCalc[m];
1081  } else {
1082  for (size_t n = 0; n <= m_mm; n++) {
1083  a1(m,n) = 0.0;
1084  }
1085  a1(m,m) = 1.0;
1086  }
1087  }
1088  a1(m_mm, m_mm) = 0.0;
1089 
1090  // Formulate the residual, resid, and the estimate for the
1091  // convergence criteria, sum
1092  double sum = 0.0;
1093  for (size_t im = 0; im < m_mm; im++) {
1094  size_t m = m_orderVectorElements[im];
1095  if (im < m_nComponents) {
1096  resid[m] = elMoles[m] - eMolesCalc[m];
1097  } else {
1098  resid[m] = 0.0;
1099  }
1100 
1101  // For equations with positive and negative coefficients,
1102  // (electronic charge), we must mitigate the convergence
1103  // criteria by a condition limited by finite precision of
1104  // inverting a matrix. Other equations with just positive
1105  // coefficients aren't limited by this.
1106  double tmp;
1107  if (m == m_eloc) {
1108  tmp = resid[m] / (elMoles[m] + elMolesTotal*1.0E-6 + options.absElemTol);
1109  } else {
1110  tmp = resid[m] / (elMoles[m] + options.absElemTol);
1111  }
1112  sum += tmp * tmp;
1113  }
1114 
1115  for (size_t m = 0; m < m_mm; m++) {
1116  if (a1(m,m) < 1.0E-50) {
1117  if (m_loglevel > 0) {
1118  writelogf(" NOTE: Diagonalizing the analytical Jac row %d\n", m);
1119  }
1120  for (size_t n = 0; n < m_mm; n++) {
1121  a1(m,n) = 0.0;
1122  }
1123  a1(m,m) = 1.0;
1124  if (resid[m] > 0.0) {
1125  resid[m] = 1.0;
1126  } else if (resid[m] < 0.0) {
1127  resid[m] = -1.0;
1128  } else {
1129  resid[m] = 0.0;
1130  }
1131  }
1132  }
1133 
1134  resid[m_mm] = n_t - n_t_calc;
1135 
1136  if (m_loglevel > 0) {
1137  writelog("Matrix:\n");
1138  for (size_t m = 0; m <= m_mm; m++) {
1139  writelog(" [");
1140  for (size_t n = 0; n <= m_mm; n++) {
1141  writelogf(" %10.5g", a1(m,n));
1142  }
1143  writelogf("] = %10.5g\n", resid[m]);
1144  }
1145  }
1146 
1147  sum += pow(resid[m_mm] /(n_t + 1.0E-15), 2);
1148  if (m_loglevel > 0) {
1149  writelogf("(it %d) Convergence = %g\n", iter, sum);
1150  }
1151 
1152  // Insist on 20x accuracy compared to the top routine. There are
1153  // instances, for ill-conditioned or singular matrices where this is
1154  // needed to move the system to a point where the matrices aren't
1155  // singular.
1156  if (sum < 0.05 * options.relTolerance) {
1157  retn = 0;
1158  break;
1159  }
1160 
1161  // Row Sum scaling
1162  for (size_t m = 0; m <= m_mm; m++) {
1163  double tmp = 0.0;
1164  for (size_t n = 0; n <= m_mm; n++) {
1165  tmp += fabs(a1(m,n));
1166  }
1167  if (m < m_mm && tmp < 1.0E-30) {
1168  if (m_loglevel > 0) {
1169  writelogf(" NOTE: Diagonalizing row %d\n", m);
1170  }
1171  for (size_t n = 0; n <= m_mm; n++) {
1172  if (n != m) {
1173  a1(m,n) = 0.0;
1174  a1(n,m) = 0.0;
1175  }
1176  }
1177  }
1178  tmp = 1.0/tmp;
1179  for (size_t n = 0; n <= m_mm; n++) {
1180  a1(m,n) *= tmp;
1181  }
1182  resid[m] *= tmp;
1183  }
1184 
1185  if (m_loglevel > 0) {
1186  writelog("Row Summed Matrix:\n");
1187  for (size_t m = 0; m <= m_mm; m++) {
1188  writelog(" [");
1189  for (size_t n = 0; n <= m_mm; n++) {
1190  writelogf(" %10.5g", a1(m,n));
1191  }
1192  writelogf("] = %10.5g\n", resid[m]);
1193  }
1194  }
1195 
1196  // Next Step: We have row-summed the equations. However, there are
1197  // some degenerate cases where two rows will be multiplies of each
1198  // other in terms of 0 < m, 0 < m part of the matrix. This occurs on
1199  // a case by case basis, and depends upon the current state of the
1200  // element potential values, which affect the concentrations of
1201  // species.
1202  //
1203  // So, the way we have found to eliminate this problem is to lump-
1204  // sum one of the rows of the matrix, except for the last column,
1205  // and stick it all on the diagonal. Then, we at least have a non-
1206  // singular matrix, and the modified equation moves the
1207  // corresponding unknown in the correct direction.
1208  //
1209  // The previous row-sum operation has made the identification of
1210  // identical rows much simpler.
1211  //
1212  // Note at least 6E-4 is necessary for the comparison. I'm guessing
1213  // 1.0E-3. If two rows are anywhere close to being equivalent, the
1214  // algorithm can get stuck in an oscillatory mode.
1215  modifiedMatrix = false;
1216  for (size_t m = 0; m < m_mm; m++) {
1217  size_t sameAsRow = npos;
1218  for (size_t im = 0; im < m; im++) {
1219  bool theSame = true;
1220  for (size_t n = 0; n < m_mm; n++) {
1221  if (fabs(a1(m,n) - a1(im,n)) > 1.0E-7) {
1222  theSame = false;
1223  break;
1224  }
1225  }
1226  if (theSame) {
1227  sameAsRow = im;
1228  }
1229  }
1230  if (sameAsRow != npos || lumpSum[m]) {
1231  if (m_loglevel > 0) {
1232  if (lumpSum[m]) {
1233  writelogf("Lump summing row %d, due to rank deficiency analysis\n", m);
1234  } else if (sameAsRow != npos) {
1235  writelogf("Identified that rows %d and %d are the same\n", m, sameAsRow);
1236  }
1237  }
1238  modifiedMatrix = true;
1239  for (size_t n = 0; n < m_mm; n++) {
1240  if (n != m) {
1241  a1(m,m) += fabs(a1(m,n));
1242  a1(m,n) = 0.0;
1243  }
1244  }
1245  }
1246  }
1247 
1248  if (m_loglevel > 0 && modifiedMatrix) {
1249  writelog("Row Summed, MODIFIED Matrix:\n");
1250  for (size_t m = 0; m <= m_mm; m++) {
1251  writelog(" [");
1252  for (size_t n = 0; n <= m_mm; n++) {
1253  writelogf(" %10.5g", a1(m,n));
1254  }
1255  writelogf("] = %10.5g\n", resid[m]);
1256  }
1257  }
1258 
1259  try {
1260  solve(a1, resid.data());
1261  } catch (CanteraError& err) {
1262  s.restoreState(state);
1263  throw CanteraError("equilibrate:estimateEP_Brinkley()",
1264  "Jacobian is singular. \nTry adding more species, "
1265  "changing the elemental composition slightly, \nor removing "
1266  "unused elements.\n\n" + err.getMessage());
1267  }
1268 
1269  // Figure out the damping coefficient: Use a delta damping
1270  // coefficient formulation: magnitude of change is capped to exp(1).
1271  beta = 1.0;
1272  for (size_t m = 0; m < m_mm; m++) {
1273  if (resid[m] > 1.0) {
1274  beta = std::min(beta, 1.0 / resid[m]);
1275  }
1276  if (resid[m] < -1.0) {
1277  beta = std::min(beta, -1.0 / resid[m]);
1278  }
1279  }
1280  if (m_loglevel > 0 && beta != 1.0) {
1281  writelogf("(it %d) Beta = %g\n", iter, beta);
1282  }
1283  }
1284  // Update the solution vector
1285  for (size_t m = 0; m < m_mm; m++) {
1286  x[m] += beta * resid[m];
1287  }
1288  n_t *= exp(beta * resid[m_mm]);
1289 
1290  if (m_loglevel > 0) {
1291  writelogf("(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter);
1292  for (size_t m = 0; m < m_mm; m++) {
1293  writelogf(" %5s %10.5g %10.5g %10.5g\n",
1294  s.elementName(m), x_old[m], x[m], resid[m]);
1295  }
1296  writelogf(" n_t %10.5g %10.5g %10.5g \n", x_old[m_mm], n_t, exp(resid[m_mm]));
1297  }
1298  }
1299  if (m_loglevel > 0) {
1300  double temp = s.temperature();
1301  double pres = s.pressure();
1302 
1303  if (retn == 0) {
1304  writelogf(" ChemEquil::estimateEP_Brinkley() SUCCESS: equilibrium found at T = %g, Pres = %g\n",
1305  temp, pres);
1306  } else {
1307  writelogf(" ChemEquil::estimateEP_Brinkley() FAILURE: equilibrium not found at T = %g, Pres = %g\n",
1308  temp, pres);
1309  }
1310  }
1311  return retn;
1312 }
1313 
1314 
1315 void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal)
1316 {
1317  if (m_eloc == npos) {
1318  return;
1319  }
1320  if (fabs(elMolesGoal[m_eloc]) > 1.0E-20) {
1321  return;
1322  }
1324  size_t maxPosEloc = npos;
1325  size_t maxNegEloc = npos;
1326  double maxPosVal = -1.0;
1327  double maxNegVal = -1.0;
1328  if (m_loglevel > 0) {
1329  for (size_t k = 0; k < m_kk; k++) {
1330  if (nAtoms(k,m_eloc) > 0.0 && m_molefractions[k] > maxPosVal && m_molefractions[k] > 0.0) {
1331  maxPosVal = m_molefractions[k];
1332  maxPosEloc = k;
1333  }
1334  if (nAtoms(k,m_eloc) < 0.0 && m_molefractions[k] > maxNegVal && m_molefractions[k] > 0.0) {
1335  maxNegVal = m_molefractions[k];
1336  maxNegEloc = k;
1337  }
1338  }
1339  }
1340 
1341  double sumPos = 0.0;
1342  double sumNeg = 0.0;
1343  for (size_t k = 0; k < m_kk; k++) {
1344  if (nAtoms(k,m_eloc) > 0.0) {
1345  sumPos += nAtoms(k,m_eloc) * m_molefractions[k];
1346  }
1347  if (nAtoms(k,m_eloc) < 0.0) {
1348  sumNeg += nAtoms(k,m_eloc) * m_molefractions[k];
1349  }
1350  }
1351  sumNeg = - sumNeg;
1352 
1353  if (sumPos >= sumNeg) {
1354  if (sumPos <= 0.0) {
1355  return;
1356  }
1357  double factor = (elMolesGoal[m_eloc] + sumNeg) / sumPos;
1358  if (m_loglevel > 0 && factor < 0.9999999999) {
1359  writelogf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n",
1360  s.speciesName(maxPosEloc),
1361  m_molefractions[maxPosEloc], m_molefractions[maxPosEloc]*factor);
1362  }
1363  for (size_t k = 0; k < m_kk; k++) {
1364  if (nAtoms(k,m_eloc) > 0.0) {
1365  m_molefractions[k] *= factor;
1366  }
1367  }
1368  } else {
1369  double factor = (-elMolesGoal[m_eloc] + sumPos) / sumNeg;
1370  if (m_loglevel > 0 && factor < 0.9999999999) {
1371  writelogf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n",
1372  s.speciesName(maxNegEloc),
1373  m_molefractions[maxNegEloc], m_molefractions[maxNegEloc]*factor);
1374  }
1375  for (size_t k = 0; k < m_kk; k++) {
1376  if (nAtoms(k,m_eloc) < 0.0) {
1377  m_molefractions[k] *= factor;
1378  }
1379  }
1380  }
1381 
1384 }
1385 
1386 } // namespace
EquilOpt options
Options controlling how the calculation is carried out.
Definition: ChemEquil.h:136
size_t nElements() const
Number of elements.
Definition: Phase.cpp:88
void restoreState(const vector_fp &state)
Restore a state saved on a previous call to saveState.
Definition: Phase.cpp:233
vector_fp m_muSS_RT
Dimensionless values of the Gibbs free energy for the standard state of each species, at the temperature and pressure of the solution (the star standard state).
Definition: ChemEquil.h:302
size_t BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, MultiPhase *mphase, std::vector< size_t > &orderVectorSpecies, std::vector< size_t > &orderVectorElements, vector_fp &formRxnMatrix)
Choose the optimum basis of species for the equilibrium calculations.
thermo_t * m_phase
Pointer to the ThermoPhase object used to initialize this object.
Definition: ChemEquil.h:145
doublereal relTolerance
Relative tolerance.
Definition: ChemEquil.h:34
const doublereal OneAtm
One atmosphere [Pa].
Definition: ct_defs.h:69
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
int setInitialMoles(thermo_t &s, vector_fp &elMoleGoal, int loglevel=0)
Estimate the initial mole numbers.
Definition: ChemEquil.cpp:182
doublereal nAtoms(size_t k, size_t m) const
number of atoms of element m in species k.
Definition: ChemEquil.h:148
void saveState(vector_fp &state) const
Save the current internal state of the phase.
Definition: Phase.cpp:221
doublereal moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition: Phase.cpp:471
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
Definition: ThermoPhase.h:461
virtual void getGibbs_RT(doublereal *grt) const
Get the nondimensional Gibbs functions for the species in their standard states at the current T and ...
Definition: ThermoPhase.h:576
size_t m_kk
number of species in the phase
Definition: ChemEquil.h:258
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:153
void addPhase(ThermoPhase *p, doublereal moles)
Add a phase to the mixture.
Definition: MultiPhase.cpp:49
doublereal m_elementTotalSum
Current value of the sum of the element abundances given the current element potentials.
Definition: ChemEquil.h:276
T clip(const T &value, const T &lower, const T &upper)
Clip value such that lower <= value <= upper.
Definition: global.h:269
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:266
STL namespace.
virtual doublereal density() const
Density (kg/m^3).
Definition: Phase.h:607
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
void initialize(thermo_t &s)
Definition: ChemEquil.cpp:64
virtual doublereal minTemp(size_t k=npos) const
Minimum temperature for which the thermodynamic data for the species or phase are valid...
Definition: ThermoPhase.h:131
doublereal enthalpy_mass() const
Specific enthalpy. Units: J/kg.
Definition: ThermoPhase.h:714
size_t m_nComponents
This is equal to the rank of the stoichiometric coefficient matrix when it is computed.
Definition: ChemEquil.h:263
vector_fp m_elementmolefracs
Current value of the element mole fractions.
Definition: ChemEquil.h:280
void ElemRearrange(size_t nComponents, const vector_fp &elementAbundances, MultiPhase *mphase, std::vector< size_t > &orderVectorSpecies, std::vector< size_t > &orderVectorElements)
Handles the potential rearrangement of the constraint equations represented by the Formula Matrix...
void equilResidual(thermo_t &s, const vector_fp &x, const vector_fp &elmtotal, vector_fp &resid, double xval, double yval, int loglevel=0)
Evaluates the residual vector F, of length m_mm.
Definition: ChemEquil.cpp:729
void update(const thermo_t &s)
Update internally stored state information.
Definition: ChemEquil.cpp:153
int iterations
Iteration counter.
Definition: ChemEquil.h:37
doublereal RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:748
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:93
virtual std::string getMessage() const
Method overridden by derived classes to format the error message.
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:159
doublereal atomicWeight(size_t m) const
Atomic weight of element m.
Definition: Phase.cpp:128
void setToEquilState(thermo_t &s, const vector_fp &x, doublereal t)
Definition: ChemEquil.cpp:132
A class for multiphase mixtures.
Definition: MultiPhase.h:57
size_t m_eloc
Index of the element id corresponding to the electric charge of each species.
Definition: ChemEquil.h:292
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:191
vector_fp m_comp
Storage of the element compositions. natom(k,m) = m_comp[k*m_mm+ m];.
Definition: ChemEquil.h:286
doublereal absElemTol
Abs Tol in element number.
Definition: ChemEquil.h:35
double m_elemFracCutoff
element fractional cutoff, below which the element will be zeroed.
Definition: ChemEquil.h:306
doublereal entropy_mass() const
Specific entropy. Units: J/kg/K.
Definition: ThermoPhase.h:724
int estimateEP_Brinkley(thermo_t &s, vector_fp &lambda, vector_fp &elMoles)
Definition: ChemEquil.cpp:844
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
virtual void setMoleFractions(const doublereal *const x)
Set the mole fractions to the specified values.
Definition: Phase.cpp:251
vector_fp & data()
Return a reference to the data vector.
Definition: Array.h:277
int estimateElementPotentials(thermo_t &s, vector_fp &lambda, vector_fp &elMolesGoal, int loglevel=0)
Generate a starting estimate for the element potentials.
Definition: ChemEquil.cpp:218
size_t m_mm
number of elements in the phase
Definition: ChemEquil.h:257
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:466
virtual void getActivityCoefficients(doublereal *ac) const
Get the array of non-dimensional molar-based activity coefficients at the current solution temperatur...
Definition: ThermoPhase.h:420
doublereal dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
Definition: utilities.h:107
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
Definition: ThermoPhase.h:245
int equilibrate(thermo_t &s, const char *XY, int loglevel=0)
Definition: ChemEquil.cpp:308
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:174
int m_loglevel
Verbosity of printed output.
Definition: ChemEquil.h:314
virtual void setPressure(doublereal p)
Set the internally stored pressure (Pa) at constant temperature and composition.
Definition: ThermoPhase.h:770
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
virtual void setTemperature(const doublereal temp)
Set the internally stored temperature of the phase (K).
Definition: Phase.h:637
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:130
vector_fp m_lambda
Current value of the dimensional element potentials.
Definition: ChemEquil.h:272
virtual doublereal refPressure() const
Returns the reference pressure in Pa.
Definition: ThermoPhase.h:116
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
void init()
Process phases and build atomic composition array.
Definition: MultiPhase.cpp:116
Contains declarations for string manipulation functions within Cantera.
vector_fp m_molefractions
Current value of the mole fractions in the single phase. length = m_kk.
Definition: ChemEquil.h:268
double calcEmoles(thermo_t &s, vector_fp &x, const double &n_t, const vector_fp &Xmol_i_calc, vector_fp &eMolesCalc, vector_fp &n_i_calc, double pressureConst)
Given a vector of dimensionless element abundances, this routine calculates the moles of the elements...
Definition: ChemEquil.cpp:808
int maxIterations
Maximum number of iterations.
Definition: ChemEquil.h:36
int _equilflag(const char *xy)
map property strings to integers
Definition: ChemEquil.cpp:21
virtual void setToEquilState(const doublereal *lambda_RT)
This method is used by the ChemEquil equilibrium solver.
Definition: ThermoPhase.h:1227
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
doublereal intEnergy_mass() const
Specific internal energy. Units: J/kg.
Definition: ThermoPhase.h:719
virtual doublereal maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
Definition: ThermoPhase.h:184
doublereal nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
Definition: Phase.cpp:161
Chemical equilibrium.
std::string elementName(size_t m) const
Name of the element with index m.
Definition: Phase.cpp:107
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition: DenseMatrix.h:50
int dampStep(thermo_t &s, vector_fp &oldx, double oldf, vector_fp &grad, vector_fp &step, vector_fp &x, double &f, vector_fp &elmols, double xval, double yval)
Find an acceptable step size and take it.
Definition: ChemEquil.cpp:690