Cantera  3.0.0
Loading...
Searching...
No Matches
MultiPhase.cpp
Go to the documentation of this file.
1/**
2 * @file MultiPhase.cpp
3 * Definitions for the @link Cantera::MultiPhase MultiPhase@endlink
4 * object that is used to set up multiphase equilibrium problems (see @ref equilGroup).
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
17
18using namespace std;
19
20namespace Cantera
21{
22
24{
25 for (size_t n = 0; n < mix.nPhases(); n++) {
26 addPhase(mix.m_phase[n], mix.m_moles[n]);
27 }
28}
29
30void MultiPhase::addPhases(vector<ThermoPhase*>& phases,
31 const vector<double>& phaseMoles)
32{
33 for (size_t n = 0; n < phases.size(); n++) {
34 addPhase(phases[n], phaseMoles[n]);
35 }
36 init();
37}
38
39void MultiPhase::addPhase(ThermoPhase* p, double moles)
40{
41 if (m_init) {
42 throw CanteraError("MultiPhase::addPhase",
43 "phases cannot be added after init() has been called.");
44 }
45
46 if (!p->compatibleWithMultiPhase()) {
47 throw CanteraError("MultiPhase::addPhase", "Phase '{}'' is not "
48 "compatible with MultiPhase equilibrium solver", p->name());
49 }
50
51 // save the pointer to the phase object
52 m_phase.push_back(p);
53
54 // store its number of moles
55 m_moles.push_back(moles);
56 m_temp_OK.push_back(true);
57
58 // update the total number of species
59 m_nsp += p->nSpecies();
60
61 // determine if this phase has new elements for each new element, add an
62 // entry in the map from names to index number + 1:
63
64 // iterate over the elements in this phase
65 for (size_t m = 0; m < p->nElements(); m++) {
66 string ename = p->elementName(m);
67
68 // if no entry is found for this element name, then it is a new element.
69 // In this case, add the name to the list of names, increment the
70 // element count, and add an entry to the name->(index+1) map.
71 if (m_enamemap.find(ename) == m_enamemap.end()) {
72 m_enamemap[ename] = m_nel + 1;
73 m_enames.push_back(ename);
74 m_atomicNumber.push_back(p->atomicNumber(m));
75
76 // Element 'E' (or 'e') is special. Note its location.
77 if (ename == "E" || ename == "e") {
78 m_eloc = m_nel;
79 }
80
81 m_nel++;
82 }
83 }
84
85 // If the mixture temperature hasn't been set, then set the temperature and
86 // pressure to the values for the phase being added. There is no good way to
87 // do this. However, this will be overridden later.
88 if (m_temp == 298.15 && p->temperature() > 2.0E-3) {
89 m_temp = p->temperature();
90 m_press = p->pressure();
91 }
92
93 // If this is a solution phase, update the minimum and maximum mixture
94 // temperatures. Stoichiometric phases are excluded, since a mixture may
95 // define multiple stoichiometric phases, each of which has thermo data
96 // valid only over a limited range. For example, a mixture might be defined
97 // to contain a phase representing water ice and one representing liquid
98 // water, only one of which should be present if the mixture represents an
99 // equilibrium state.
100 if (p->nSpecies() > 1) {
101 m_Tmin = std::max(p->minTemp(), m_Tmin);
102 m_Tmax = std::min(p->maxTemp(), m_Tmax);
103 }
104}
105
107{
108 if (m_init) {
109 return;
110 }
111
112 // allocate space for the atomic composition matrix
113 m_atoms.resize(m_nel, m_nsp, 0.0);
114 m_moleFractions.resize(m_nsp, 0.0);
115 m_elemAbundances.resize(m_nel, 0.0);
116
117 // iterate over the elements
118 // -> fill in m_atoms(m,k), m_snames(k), m_spphase(k), m_spstart(ip)
119 for (size_t m = 0; m < m_nel; m++) {
120 size_t k = 0;
121 // iterate over the phases
122 for (size_t ip = 0; ip < nPhases(); ip++) {
123 ThermoPhase* p = m_phase[ip];
124 size_t nsp = p->nSpecies();
125 size_t mlocal = p->elementIndex(m_enames[m]);
126 for (size_t kp = 0; kp < nsp; kp++) {
127 if (mlocal != npos) {
128 m_atoms(m, k) = p->nAtoms(kp, mlocal);
129 }
130 if (m == 0) {
131 m_snames.push_back(p->speciesName(kp));
132 if (kp == 0) {
133 m_spstart.push_back(m_spphase.size());
134 }
135 m_spphase.push_back(ip);
136 }
137 k++;
138 }
139 }
140 }
141
142 // set the initial composition within each phase to the
143 // mole fractions stored in the phase objects
144 m_init = true;
146 updatePhases();
147}
148
150{
151 if (!m_init) {
152 init();
153 }
154 m_phase[n]->setTemperature(m_temp);
155 m_phase[n]->setMoleFractions_NoNorm(&m_moleFractions[m_spstart[n]]);
156 m_phase[n]->setPressure(m_press);
157 return *m_phase[n];
158}
159
160void MultiPhase::checkPhaseIndex(size_t m) const
161{
162 if (m >= nPhases()) {
163 throw IndexError("MultiPhase::checkPhaseIndex", "phase", m, nPhases()-1);
164 }
165}
166
168{
169 if (nPhases() > mm) {
170 throw ArraySizeError("MultiPhase::checkPhaseIndex", mm, nPhases());
171 }
172}
173
174double MultiPhase::speciesMoles(size_t k) const
175{
176 size_t ip = m_spphase[k];
177 return m_moles[ip]*m_moleFractions[k];
178}
179
180double MultiPhase::elementMoles(size_t m) const
181{
182 double sum = 0.0;
183 for (size_t i = 0; i < nPhases(); i++) {
184 double phasesum = 0.0;
185 size_t nsp = m_phase[i]->nSpecies();
186 for (size_t ik = 0; ik < nsp; ik++) {
187 size_t k = speciesIndex(ik, i);
188 phasesum += m_atoms(m,k)*m_moleFractions[k];
189 }
190 sum += phasesum * m_moles[i];
191 }
192 return sum;
193}
194
195double MultiPhase::charge() const
196{
197 double sum = 0.0;
198 for (size_t i = 0; i < nPhases(); i++) {
199 sum += phaseCharge(i);
200 }
201 return sum;
202}
203
204size_t MultiPhase::speciesIndex(const string& speciesName, const string& phaseName)
205{
206 if (!m_init) {
207 init();
208 }
209 size_t p = phaseIndex(phaseName);
210 if (p == npos) {
211 throw CanteraError("MultiPhase::speciesIndex", "phase not found: " + phaseName);
212 }
213 size_t k = m_phase[p]->speciesIndex(speciesName);
214 if (k == npos) {
215 throw CanteraError("MultiPhase::speciesIndex", "species not found: " + speciesName);
216 }
217 return m_spstart[p] + k;
218}
219
220double MultiPhase::phaseCharge(size_t p) const
221{
222 double phasesum = 0.0;
223 size_t nsp = m_phase[p]->nSpecies();
224 for (size_t ik = 0; ik < nsp; ik++) {
225 size_t k = speciesIndex(ik, p);
226 phasesum += m_phase[p]->charge(ik)*m_moleFractions[k];
227 }
228 return Faraday*phasesum*m_moles[p];
229}
230
231void MultiPhase::getChemPotentials(double* mu) const
232{
233 updatePhases();
234 size_t loc = 0;
235 for (size_t i = 0; i < nPhases(); i++) {
236 m_phase[i]->getChemPotentials(mu + loc);
237 loc += m_phase[i]->nSpecies();
238 }
239}
240
241void MultiPhase::getValidChemPotentials(double not_mu, double* mu, bool standard) const
242{
243 updatePhases();
244 // iterate over the phases
245 size_t loc = 0;
246 for (size_t i = 0; i < nPhases(); i++) {
247 if (tempOK(i) || m_phase[i]->nSpecies() > 1) {
248 if (!standard) {
249 m_phase[i]->getChemPotentials(mu + loc);
250 } else {
251 m_phase[i]->getStandardChemPotentials(mu + loc);
252 }
253 } else {
254 fill(mu + loc, mu + loc + m_phase[i]->nSpecies(), not_mu);
255 }
256 loc += m_phase[i]->nSpecies();
257 }
258}
259
260bool MultiPhase::solutionSpecies(size_t k) const
261{
262 if (m_phase[m_spphase[k]]->nSpecies() > 1) {
263 return true;
264 } else {
265 return false;
266 }
267}
268
269double MultiPhase::gibbs() const
270{
271 double sum = 0.0;
272 updatePhases();
273 for (size_t i = 0; i < nPhases(); i++) {
274 if (m_moles[i] > 0.0) {
275 sum += m_phase[i]->gibbs_mole() * m_moles[i];
276 }
277 }
278 return sum;
279}
280
282{
283 double sum = 0.0;
284 updatePhases();
285 for (size_t i = 0; i < nPhases(); i++) {
286 if (m_moles[i] > 0.0) {
287 sum += m_phase[i]->enthalpy_mole() * m_moles[i];
288 }
289 }
290 return sum;
291}
292
294{
295 double sum = 0.0;
296 updatePhases();
297 for (size_t i = 0; i < nPhases(); i++) {
298 if (m_moles[i] > 0.0) {
299 sum += m_phase[i]->intEnergy_mole() * m_moles[i];
300 }
301 }
302 return sum;
303}
304
306{
307 double sum = 0.0;
308 updatePhases();
309 for (size_t i = 0; i < nPhases(); i++) {
310 if (m_moles[i] > 0.0) {
311 sum += m_phase[i]->entropy_mole() * m_moles[i];
312 }
313 }
314 return sum;
315}
316
317double MultiPhase::cp() const
318{
319 double sum = 0.0;
320 updatePhases();
321 for (size_t i = 0; i < nPhases(); i++) {
322 if (m_moles[i] > 0.0) {
323 sum += m_phase[i]->cp_mole() * m_moles[i];
324 }
325 }
326 return sum;
327}
328
329void MultiPhase::setPhaseMoleFractions(const size_t n, const double* const x)
330{
331 if (!m_init) {
332 init();
333 }
334 ThermoPhase* p = m_phase[n];
336 size_t istart = m_spstart[n];
337 for (size_t k = 0; k < p->nSpecies(); k++) {
338 m_moleFractions[istart+k] = x[k];
339 }
340}
341
343{
344 size_t kk = nSpecies();
345 vector<double> moles(kk, 0.0);
346 for (size_t k = 0; k < kk; k++) {
347 moles[k] = std::max(getValue(xMap, speciesName(k), 0.0), 0.0);
348 }
349 setMoles(moles.data());
350}
351
352void MultiPhase::setMolesByName(const string& x)
353{
354 // build the composition map from the string, and then set the moles.
356 setMolesByName(xx);
357}
358
359void MultiPhase::getMoles(double* molNum) const
360{
361 // First copy in the mole fractions
362 copy(m_moleFractions.begin(), m_moleFractions.end(), molNum);
363 double* dtmp = molNum;
364 for (size_t ip = 0; ip < nPhases(); ip++) {
365 double phasemoles = m_moles[ip];
366 ThermoPhase* p = m_phase[ip];
367 size_t nsp = p->nSpecies();
368 for (size_t ik = 0; ik < nsp; ik++) {
369 *(dtmp++) *= phasemoles;
370 }
371 }
372}
373
374void MultiPhase::setMoles(const double* n)
375{
376 if (!m_init) {
377 init();
378 }
379 size_t loc = 0;
380 size_t k = 0;
381 for (size_t ip = 0; ip < nPhases(); ip++) {
382 ThermoPhase* p = m_phase[ip];
383 size_t nsp = p->nSpecies();
384 double phasemoles = 0.0;
385 for (size_t ik = 0; ik < nsp; ik++) {
386 phasemoles += n[k];
387 k++;
388 }
389 m_moles[ip] = phasemoles;
390 if (nsp > 1) {
391 if (phasemoles > 0.0) {
392 p->setState_TPX(m_temp, m_press, n + loc);
394 } else {
396 }
397 } else {
398 m_moleFractions[loc] = 1.0;
399 }
400 loc += nsp;
401 }
402}
403
404void MultiPhase::addSpeciesMoles(const int indexS, const double addedMoles)
405{
406 vector<double> tmpMoles(m_nsp, 0.0);
407 getMoles(tmpMoles.data());
408 tmpMoles[indexS] += addedMoles;
409 tmpMoles[indexS] = std::max(tmpMoles[indexS], 0.0);
410 setMoles(tmpMoles.data());
411}
412
413void MultiPhase::setState_TP(const double T, const double Pres)
414{
415 if (!m_init) {
416 init();
417 }
418 m_temp = T;
419 m_press = Pres;
420 updatePhases();
421}
422
423void MultiPhase::setState_TPMoles(const double T, const double Pres, const double* n)
424{
425 m_temp = T;
426 m_press = Pres;
427 setMoles(n);
428}
429
430void MultiPhase::getElemAbundances(double* elemAbundances) const
431{
433 for (size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) {
434 elemAbundances[eGlobal] = m_elemAbundances[eGlobal];
435 }
436}
437
439{
440 size_t loc = 0;
441 double spMoles;
442 for (size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) {
443 m_elemAbundances[eGlobal] = 0.0;
444 }
445 for (size_t ip = 0; ip < nPhases(); ip++) {
446 ThermoPhase* p = m_phase[ip];
447 size_t nspPhase = p->nSpecies();
448 double phasemoles = m_moles[ip];
449 for (size_t ik = 0; ik < nspPhase; ik++) {
450 size_t kGlobal = loc + ik;
451 spMoles = m_moleFractions[kGlobal] * phasemoles;
452 for (size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) {
453 m_elemAbundances[eGlobal] += m_atoms(eGlobal, kGlobal) * spMoles;
454 }
455 }
456 loc += nspPhase;
457 }
458}
459
460double MultiPhase::volume() const
461{
462 double sum = 0;
463 for (size_t i = 0; i < nPhases(); i++) {
464 double vol = 1.0/m_phase[i]->molarDensity();
465 sum += m_moles[i] * vol;
466 }
467 return sum;
468}
469
470double MultiPhase::equilibrate_MultiPhaseEquil(int XY, double err, int maxsteps,
471 int maxiter, int loglevel)
472{
473 bool strt = false;
474 double dta = 0.0;
475 if (!m_init) {
476 init();
477 }
478
479 if (XY == TP) {
480 // create an equilibrium manager
481 MultiPhaseEquil e(this);
482 return e.equilibrate(XY, err, maxsteps, loglevel);
483 } else if (XY == HP) {
484 double h0 = enthalpy();
485 double Tlow = 0.5*m_Tmin; // lower bound on T
486 double Thigh = 2.0*m_Tmax; // upper bound on T
487 double Hlow = Undef, Hhigh = Undef;
488 for (int n = 0; n < maxiter; n++) {
489 // if 'strt' is false, the current composition will be used as
490 // the starting estimate; otherwise it will be estimated
491 MultiPhaseEquil e(this, strt);
492 // start with a loose error tolerance, but tighten it as we get
493 // close to the final temperature
494
495 try {
496 e.equilibrate(TP, err, maxsteps, loglevel);
497 double hnow = enthalpy();
498 // the equilibrium enthalpy monotonically increases with T;
499 // if the current value is below the target, the we know the
500 // current temperature is too low. Set
501 if (hnow < h0) {
502 if (m_temp > Tlow) {
503 Tlow = m_temp;
504 Hlow = hnow;
505 }
506 } else {
507 // the current enthalpy is greater than the target;
508 // therefore the current temperature is too high.
509 if (m_temp < Thigh) {
510 Thigh = m_temp;
511 Hhigh = hnow;
512 }
513 }
514 double dt;
515 if (Hlow != Undef && Hhigh != Undef) {
516 double cpb = (Hhigh - Hlow)/(Thigh - Tlow);
517 dt = (h0 - hnow)/cpb;
518 dta = fabs(dt);
519 double dtmax = 0.5*fabs(Thigh - Tlow);
520 if (dta > dtmax) {
521 dt *= dtmax/dta;
522 }
523 } else {
524 double tnew = sqrt(Tlow*Thigh);
525 dt = tnew - m_temp;
526 }
527
528 double herr = fabs((h0 - hnow)/h0);
529
530 if (herr < err) {
531 return err;
532 }
533 double tnew = m_temp + dt;
534 if (tnew < 0.0) {
535 tnew = 0.5*m_temp;
536 }
537 setTemperature(tnew);
538
539 // if the size of Delta T is not too large, use
540 // the current composition as the starting estimate
541 if (dta < 100.0) {
542 strt = false;
543 }
544
545 } catch (CanteraError&) {
546 if (!strt) {
547 strt = true;
548 } else {
549 double tnew = 0.5*(m_temp + Thigh);
550 if (fabs(tnew - m_temp) < 1.0) {
551 tnew = m_temp + 1.0;
552 }
553 setTemperature(tnew);
554 }
555 }
556 }
557 throw CanteraError("MultiPhase::equilibrate_MultiPhaseEquil",
558 "No convergence for T");
559 } else if (XY == SP) {
560 double s0 = entropy();
561 double Tlow = 1.0; // lower bound on T
562 double Thigh = 1.0e6; // upper bound on T
563 for (int n = 0; n < maxiter; n++) {
564 MultiPhaseEquil e(this, strt);
565
566 try {
567 e.equilibrate(TP, err, maxsteps, loglevel);
568 double snow = entropy();
569 if (snow < s0) {
570 Tlow = std::max(Tlow, m_temp);
571 } else {
572 Thigh = std::min(Thigh, m_temp);
573 }
574 double dt = (s0 - snow)*m_temp/cp();
575 double dtmax = 0.5*fabs(Thigh - Tlow);
576 dtmax = (dtmax > 500.0 ? 500.0 : dtmax);
577 dta = fabs(dt);
578 if (dta > dtmax) {
579 dt *= dtmax/dta;
580 }
581 if (dta < 1.0e-4) {
582 return err;
583 }
584 double tnew = m_temp + dt;
585 setTemperature(tnew);
586
587 // if the size of Delta T is not too large, use
588 // the current composition as the starting estimate
589 if (dta < 100.0) {
590 strt = false;
591 }
592 } catch (CanteraError&) {
593 if (!strt) {
594 strt = true;
595 } else {
596 double tnew = 0.5*(m_temp + Thigh);
597 setTemperature(tnew);
598 }
599 }
600 }
601 throw CanteraError("MultiPhase::equilibrate_MultiPhaseEquil",
602 "No convergence for T");
603 } else if (XY == TV) {
604 double v0 = volume();
605 bool start = true;
606 for (int n = 0; n < maxiter; n++) {
607 double pnow = pressure();
608 MultiPhaseEquil e(this, start);
609 start = false;
610
611 e.equilibrate(TP, err, maxsteps, loglevel);
612 double vnow = volume();
613 double verr = fabs((v0 - vnow)/v0);
614
615 if (verr < err) {
616 return err;
617 }
618 // find dV/dP
619 setPressure(pnow*1.01);
620 double dVdP = (volume() - vnow)/(0.01*pnow);
621 setPressure(pnow + 0.5*(v0 - vnow)/dVdP);
622 }
623 } else {
624 throw CanteraError("MultiPhase::equilibrate_MultiPhaseEquil",
625 "unknown option");
626 }
627 return -1.0;
628}
629
630void MultiPhase::equilibrate(const string& XY, const string& solver,
631 double rtol, int max_steps, int max_iter,
632 int estimate_equil, int log_level)
633{
634 // Save the initial state so that it can be restored in case one of the
635 // solvers fails
636 vector<double> initial_moleFractions = m_moleFractions;
637 vector<double> initial_moles = m_moles;
638 double initial_T = m_temp;
639 double initial_P = m_press;
640 int ixy = _equilflag(XY.c_str());
641 if (solver == "auto" || solver == "vcs") {
642 try {
643 debuglog("Trying VCS equilibrium solver\n", log_level);
644 vcs_MultiPhaseEquil eqsolve(this, log_level-1);
645 int ret = eqsolve.equilibrate(ixy, estimate_equil, log_level-1,
646 rtol, max_steps);
647 if (ret) {
648 throw CanteraError("MultiPhase::equilibrate",
649 "VCS solver failed. Return code: {}", ret);
650 }
651 debuglog("VCS solver succeeded\n", log_level);
652 return;
653 } catch (std::exception& err) {
654 debuglog("VCS solver failed.\n", log_level);
655 debuglog(err.what(), log_level);
656 m_moleFractions = initial_moleFractions;
657 m_moles = initial_moles;
658 m_temp = initial_T;
659 m_press = initial_P;
660 updatePhases();
661 if (solver == "auto") {
662 } else {
663 throw;
664 }
665 }
666 }
667
668 if (solver == "auto" || solver == "gibbs") {
669 try {
670 debuglog("Trying MultiPhaseEquil (Gibbs) equilibrium solver\n",
671 log_level);
672 equilibrate_MultiPhaseEquil(ixy, rtol, max_steps, max_iter,
673 log_level-1);
674 debuglog("MultiPhaseEquil solver succeeded\n", log_level);
675 return;
676 } catch (std::exception& err) {
677 debuglog("MultiPhaseEquil solver failed.\n", log_level);
678 debuglog(err.what(), log_level);
679 m_moleFractions = initial_moleFractions;
680 m_moles = initial_moles;
681 m_temp = initial_T;
682 m_press = initial_P;
683 updatePhases();
684 throw;
685 }
686 }
687
688 if (solver != "auto") {
689 throw CanteraError("MultiPhase::equilibrate",
690 "Invalid solver specified: '" + solver + "'");
691 }
692}
693
694void MultiPhase::setTemperature(const double T)
695{
696 if (!m_init) {
697 init();
698 }
699 m_temp = T;
700 updatePhases();
701}
702
704{
705 if (m >= m_nel) {
706 throw IndexError("MultiPhase::checkElementIndex", "elements", m, m_nel-1);
707 }
708}
709
711{
712 if (m_nel > mm) {
713 throw ArraySizeError("MultiPhase::checkElementArraySize", mm, m_nel);
714 }
715}
716
717string MultiPhase::elementName(size_t m) const
718{
719 return m_enames[m];
720}
721
722size_t MultiPhase::elementIndex(const string& name) const
723{
724 for (size_t e = 0; e < m_nel; e++) {
725 if (m_enames[e] == name) {
726 return e;
727 }
728 }
729 return npos;
730}
731
733{
734 if (k >= m_nsp) {
735 throw IndexError("MultiPhase::checkSpeciesIndex", "species", k, m_nsp-1);
736 }
737}
738
740{
741 if (m_nsp > kk) {
742 throw ArraySizeError("MultiPhase::checkSpeciesArraySize", kk, m_nsp);
743 }
744}
745
746string MultiPhase::speciesName(const size_t k) const
747{
748 return m_snames[k];
749}
750
751double MultiPhase::nAtoms(const size_t kGlob, const size_t mGlob) const
752{
753 return m_atoms(mGlob, kGlob);
754}
755
756void MultiPhase::getMoleFractions(double* const x) const
757{
758 std::copy(m_moleFractions.begin(), m_moleFractions.end(), x);
759}
760
761string MultiPhase::phaseName(const size_t iph) const
762{
763 return m_phase[iph]->name();
764}
765
766int MultiPhase::phaseIndex(const string& pName) const
767{
768 for (int iph = 0; iph < (int) nPhases(); iph++) {
769 if (m_phase[iph]->name() == pName) {
770 return iph;
771 }
772 }
773 return -1;
774}
775
776double MultiPhase::phaseMoles(const size_t n) const
777{
778 return m_moles[n];
779}
780
781void MultiPhase::setPhaseMoles(const size_t n, const double moles)
782{
783 m_moles[n] = moles;
784}
785
786size_t MultiPhase::speciesPhaseIndex(const size_t kGlob) const
787{
788 return m_spphase[kGlob];
789}
790
791double MultiPhase::moleFraction(const size_t kGlob) const
792{
793 return m_moleFractions[kGlob];
794}
795
796bool MultiPhase::tempOK(const size_t p) const
797{
798 return m_temp_OK[p];
799}
800
802{
803 size_t loc = 0;
804 for (size_t ip = 0; ip < nPhases(); ip++) {
805 ThermoPhase* p = m_phase[ip];
806 p->getMoleFractions(m_moleFractions.data() + loc);
807 loc += p->nSpecies();
808 }
810}
811
813{
814 size_t loc = 0;
815 for (size_t p = 0; p < nPhases(); p++) {
816 m_phase[p]->setState_TPX(m_temp, m_press, m_moleFractions.data() + loc);
817 loc += m_phase[p]->nSpecies();
818 m_temp_OK[p] = true;
819 if (m_temp < m_phase[p]->minTemp() || m_temp > m_phase[p]->maxTemp()) {
820 m_temp_OK[p] = false;
821 }
822 }
823}
824
825std::ostream& operator<<(std::ostream& s, MultiPhase& x)
826{
827 x.updatePhases();
828 for (size_t ip = 0; ip < x.nPhases(); ip++) {
829 if (x.phase(ip).name() != "") {
830 s << "*************** " << x.phase(ip).name() << " *****************" << std::endl;
831 } else {
832 s << "*************** Phase " << ip << " *****************" << std::endl;
833 }
834 s << "Moles: " << x.phaseMoles(ip) << std::endl;
835
836 s << x.phase(ip).report() << std::endl;
837 }
838 return s;
839}
840
841}
Chemical equilibrium.
Headers for the MultiPhase object that is used to set up multiphase equilibrium problems (see Chemica...
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
Array size error.
Base class for exceptions thrown by Cantera classes.
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
An array index is out of range.
Multiphase chemical equilibrium solver.
A class for multiphase mixtures.
Definition MultiPhase.h:57
void init()
Process phases and build atomic composition array.
void checkPhaseIndex(size_t m) const
Check that the specified phase index is in range Throws an exception if m is greater than nPhases()
size_t speciesIndex(size_t k, size_t p) const
Return the global index of the species belonging to phase number p with local index k within the phas...
Definition MultiPhase.h:226
bool solutionSpecies(size_t kGlob) const
Return true is species kGlob is a species in a multicomponent solution phase.
vector< ThermoPhase * > m_phase
Vector of the ThermoPhase pointers.
Definition MultiPhase.h:571
double nAtoms(const size_t kGlob, const size_t mGlob) const
Returns the Number of atoms of global element mGlob in global species kGlob.
void setMolesByName(const Composition &xMap)
Set the number of moles of species in the mixture.
void setMoles(const double *n)
Sets all of the global species mole numbers.
DenseMatrix m_atoms
Global Stoichiometric Coefficient array.
Definition MultiPhase.h:579
double gibbs() const
The Gibbs function of the mixture [J].
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range.
size_t m_nel
Number of distinct elements in all of the phases.
Definition MultiPhase.h:624
double speciesMoles(size_t kGlob) const
Returns the moles of global species k. units = kmol.
void getValidChemPotentials(double not_mu, double *mu, bool standard=false) const
Returns a vector of Valid chemical potentials.
double m_temp
Current value of the temperature (kelvin)
Definition MultiPhase.h:618
void calcElemAbundances() const
Calculate the element abundance vector.
size_t nSpecies() const
Number of species, summed over all phases.
Definition MultiPhase.h:123
double enthalpy() const
The enthalpy of the mixture [J].
double pressure() const
Pressure [Pa].
Definition MultiPhase.h:388
vector< size_t > m_spstart
Vector of ints containing of first species index in the global list of species for each phase.
Definition MultiPhase.h:598
vector< size_t > m_spphase
Mapping between the global species number and the phase ID.
Definition MultiPhase.h:590
void getMoles(double *molNum) const
Get the mole numbers of all species in the multiphase object.
double minTemp() const
Minimum temperature for which all solution phases have valid thermo data.
Definition MultiPhase.h:246
vector< double > m_moleFractions
Locally stored vector of mole fractions of all species comprising the MultiPhase object.
Definition MultiPhase.h:583
vector< double > m_elemAbundances
Vector of element abundances.
Definition MultiPhase.h:656
double equilibrate_MultiPhaseEquil(int XY, double err, int maxsteps, int maxiter, int loglevel)
Set the mixture to a state of chemical equilibrium using the MultiPhaseEquil solver.
void checkPhaseArraySize(size_t mm) const
Check that an array size is at least nPhases() Throws an exception if mm is less than nPhases().
vector< bool > m_temp_OK
Vector of bools indicating whether temperatures are ok for phases.
Definition MultiPhase.h:641
double phaseCharge(size_t p) const
Charge (Coulombs) of phase with index p.
int phaseIndex(const string &pName) const
Returns the index, given the phase name.
size_t nPhases() const
Number of phases.
Definition MultiPhase.h:425
size_t m_eloc
Global ID of the element corresponding to the electronic charge.
Definition MultiPhase.h:634
double entropy() const
The entropy of the mixture [J/K].
void getChemPotentials(double *mu) const
Returns a vector of Chemical potentials.
double moleFraction(const size_t kGlob) const
Returns the mole fraction of global species k.
double m_press
Current value of the pressure (Pa)
Definition MultiPhase.h:621
bool tempOK(size_t p) const
Return true if the phase p has valid thermo data for the current temperature.
void addPhases(vector< ThermoPhase * > &phases, const vector< double > &phaseMoles)
Add a vector of phases to the mixture.
void addPhase(ThermoPhase *p, double moles)
Add a phase to the mixture.
map< string, size_t > m_enamemap
Returns the global element index, given the element string name.
Definition MultiPhase.h:615
void addSpeciesMoles(const int indexS, const double addedMoles)
Adds moles of a certain species to the mixture.
size_t elementIndex(const string &name) const
Returns the index of the element with name name.
void checkElementArraySize(size_t mm) const
Check that an array size is at least nElements().
size_t speciesPhaseIndex(const size_t kGlob) const
Returns the phase index of the Kth "global" species.
void setPressure(double P)
Set the pressure [Pa].
Definition MultiPhase.h:403
void setState_TPMoles(const double T, const double Pres, const double *Moles)
Set the state of the underlying ThermoPhase objects in one call.
vector< string > m_enames
String names of the global elements.
Definition MultiPhase.h:602
string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
double phaseMoles(const size_t n) const
Return the number of moles in phase n.
void getMoleFractions(double *const x) const
Returns the global Species mole fractions.
size_t m_nsp
Number of distinct species in all of the phases.
Definition MultiPhase.h:627
void setPhaseMoleFractions(const size_t n, const double *const x)
Set the Mole fractions of the nth phase.
vector< int > m_atomicNumber
Atomic number of each global element.
Definition MultiPhase.h:605
double volume() const
The total mixture volume [m^3].
void uploadMoleFractionsFromPhases()
Update the locally-stored composition within this object to match the current compositions of the pha...
bool m_init
True if the init() routine has been called, and the MultiPhase frozen.
Definition MultiPhase.h:630
vector< string > m_snames
Vector of species names in the problem.
Definition MultiPhase.h:609
double m_Tmin
Minimum temperature for which thermo parameterizations are valid.
Definition MultiPhase.h:645
void updatePhases() const
Set the states of the phase objects to the locally-stored state within this MultiPhase object.
void getElemAbundances(double *elemAbundances) const
Retrieves a vector of element abundances.
void checkSpeciesArraySize(size_t kk) const
Check that an array size is at least nSpecies().
double charge() const
Total charge summed over all phases (Coulombs).
string phaseName(const size_t iph) const
Returns the name of the n'th phase.
double cp() const
Heat capacity at constant pressure [J/K].
void setTemperature(const double T)
Set the temperature [K].
void setState_TP(const double T, const double Pres)
Set the state of the underlying ThermoPhase objects in one call.
void checkElementIndex(size_t m) const
Check that the specified element index is in range.
ThermoPhase & phase(size_t n)
Return a reference to phase n.
double IntEnergy() const
The internal energy of the mixture [J].
string elementName(size_t m) const
Returns the name of the global element m.
vector< double > m_moles
Vector of the number of moles in each phase.
Definition MultiPhase.h:568
void setPhaseMoles(const size_t n, const double moles)
Set the number of moles of phase with index n.
double elementMoles(size_t m) const
Total moles of global element m, summed over all phases.
double m_Tmax
Minimum temperature for which thermo parameterizations are valid.
Definition MultiPhase.h:649
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:245
int atomicNumber(size_t m) const
Atomic number of element m.
Definition Phase.cpp:86
double temperature() const
Temperature (K).
Definition Phase.h:662
size_t elementIndex(const string &name) const
Return the index of element named 'name'.
Definition Phase.cpp:55
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 nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
Definition Phase.cpp:103
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
string name() const
Return the name of the phase.
Definition Phase.cpp:20
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 void setState_TPX(double t, double p, const double *x)
Set the temperature (K), pressure (Pa), and mole fractions.
virtual string report(bool show_thermo=true, double threshold=-1e-14) const
returns a summary of the state of the phase as a string
virtual double maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
Cantera's Interface to the Multiphase chemical equilibrium solver.
int equilibrate(int XY, int estimateEquil=0, int printLvl=0, double err=1.0e-6, int maxsteps=VCS_MAXSTEPS, int loglevel=-99)
Equilibrate the solution using the current element abundances stored in the MultiPhase object.
Composition parseCompString(const string &ss, const vector< string > &names)
Parse a composition string into a map consisting of individual key:composition pairs.
void equilibrate(const string &XY, const string &solver="auto", double rtol=1e-9, int max_steps=50000, int max_iter=100, int estimate_equil=0, int log_level=0)
Equilibrate a MultiPhase object.
virtual bool compatibleWithMultiPhase() const
Indicates whether this phase type can be used with class MultiPhase for equilibrium calculations.
void debuglog(const string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Definition global.h:158
const double Faraday
Faraday constant [C/kmol].
Definition ct_defs.h:131
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
const double Undef
Fairly random number to be used to initialize variables against to see if they are subsequently defin...
Definition ct_defs.h:164
int _equilflag(const char *xy)
map property strings to integers
Definition ChemEquil.cpp:20
std::ostream & operator<<(std::ostream &s, const Array2D &m)
Output the current contents of the Array2D object.
Definition Array.cpp:120
const U & getValue(const map< T, U > &m, const T &key, const U &default_val)
Const accessor for a value in a map.
Definition utilities.h:190
map< string, double > Composition
Map from string names to doubles.
Definition ct_defs.h:184
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...
Interface class for the vcsnonlinear solver.