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