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