Cantera  3.1.0a1
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 
16 #include "cantera/base/utilities.h"
17 
18 using namespace std;
19 
20 namespace Cantera
21 {
22 
23 void MultiPhase::addPhases(MultiPhase& mix)
24 {
25  for (size_t n = 0; n < mix.nPhases(); n++) {
26  addPhase(mix.m_phase[n], mix.m_moles[n]);
27  }
28 }
29 
30 void 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 
39 void 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 
106 void MultiPhase::init()
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;
145  uploadMoleFractionsFromPhases();
146  updatePhases();
147 }
148 
149 ThermoPhase& MultiPhase::phase(size_t n)
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 
160 void MultiPhase::checkPhaseIndex(size_t m) const
161 {
162  if (m >= nPhases()) {
163  throw IndexError("MultiPhase::checkPhaseIndex", "phase", m, nPhases()-1);
164  }
165 }
166 
167 void MultiPhase::checkPhaseArraySize(size_t mm) const
168 {
169  if (nPhases() > mm) {
170  throw ArraySizeError("MultiPhase::checkPhaseIndex", mm, nPhases());
171  }
172 }
173 
174 double MultiPhase::speciesMoles(size_t k) const
175 {
176  size_t ip = m_spphase[k];
177  return m_moles[ip]*m_moleFractions[k];
178 }
179 
180 double 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 
195 double 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 
204 size_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 
220 double 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 
231 void 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 
241 void 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 
260 bool 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 
269 double 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 
281 double MultiPhase::enthalpy() const
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 
293 double MultiPhase::IntEnergy() const
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 
305 double MultiPhase::entropy() const
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 
317 double 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 
329 void MultiPhase::setPhaseMoleFractions(const size_t n, const double* const x)
330 {
331  if (!m_init) {
332  init();
333  }
334  ThermoPhase* p = m_phase[n];
335  p->setState_TPX(m_temp, m_press, x);
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 
342 void MultiPhase::setMolesByName(const Composition& xMap)
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 
352 void MultiPhase::setMolesByName(const string& x)
353 {
354  // build the composition map from the string, and then set the moles.
355  Composition xx = parseCompString(x, m_snames);
356  setMolesByName(xx);
357 }
358 
359 void 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 
374 void 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);
393  p->getMoleFractions(&m_moleFractions[loc]);
394  } else {
395  p->getMoleFractions(&m_moleFractions[loc]);
396  }
397  } else {
398  m_moleFractions[loc] = 1.0;
399  }
400  loc += nsp;
401  }
402 }
403 
404 void 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 
413 void 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 
423 void 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 
430 void MultiPhase::getElemAbundances(double* elemAbundances) const
431 {
432  calcElemAbundances();
433  for (size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) {
434  elemAbundances[eGlobal] = m_elemAbundances[eGlobal];
435  }
436 }
437 
438 void MultiPhase::calcElemAbundances() const
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 
460 double 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 
470 double 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 
630 void 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 
694 void MultiPhase::setTemperature(const double T)
695 {
696  if (!m_init) {
697  init();
698  }
699  m_temp = T;
700  updatePhases();
701 }
702 
703 void MultiPhase::checkElementIndex(size_t m) const
704 {
705  if (m >= m_nel) {
706  throw IndexError("MultiPhase::checkElementIndex", "elements", m, m_nel-1);
707  }
708 }
709 
710 void MultiPhase::checkElementArraySize(size_t mm) const
711 {
712  if (m_nel > mm) {
713  throw ArraySizeError("MultiPhase::checkElementArraySize", mm, m_nel);
714  }
715 }
716 
717 string MultiPhase::elementName(size_t m) const
718 {
719  return m_enames[m];
720 }
721 
722 size_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 
732 void MultiPhase::checkSpeciesIndex(size_t k) const
733 {
734  if (k >= m_nsp) {
735  throw IndexError("MultiPhase::checkSpeciesIndex", "species", k, m_nsp-1);
736  }
737 }
738 
739 void MultiPhase::checkSpeciesArraySize(size_t kk) const
740 {
741  if (m_nsp > kk) {
742  throw ArraySizeError("MultiPhase::checkSpeciesArraySize", kk, m_nsp);
743  }
744 }
745 
746 string MultiPhase::speciesName(const size_t k) const
747 {
748  return m_snames[k];
749 }
750 
751 double MultiPhase::nAtoms(const size_t kGlob, const size_t mGlob) const
752 {
753  return m_atoms(mGlob, kGlob);
754 }
755 
756 void MultiPhase::getMoleFractions(double* const x) const
757 {
758  std::copy(m_moleFractions.begin(), m_moleFractions.end(), x);
759 }
760 
761 string MultiPhase::phaseName(const size_t iph) const
762 {
763  return m_phase[iph]->name();
764 }
765 
766 int 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 
776 double MultiPhase::phaseMoles(const size_t n) const
777 {
778  return m_moles[n];
779 }
780 
781 void MultiPhase::setPhaseMoles(const size_t n, const double moles)
782 {
783  m_moles[n] = moles;
784 }
785 
786 size_t MultiPhase::speciesPhaseIndex(const size_t kGlob) const
787 {
788  return m_spphase[kGlob];
789 }
790 
791 double MultiPhase::moleFraction(const size_t kGlob) const
792 {
793  return m_moleFractions[kGlob];
794 }
795 
796 bool MultiPhase::tempOK(const size_t p) const
797 {
798  return m_temp_OK[p];
799 }
800 
801 void MultiPhase::uploadMoleFractionsFromPhases()
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  }
809  calcElemAbundances();
810 }
811 
812 void MultiPhase::updatePhases() const
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 
825 std::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.
Definition: ctexceptions.h:135
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
An array index is out of range.
Definition: ctexceptions.h:165
Multiphase chemical equilibrium solver.
A class for multiphase mixtures.
Definition: MultiPhase.h:57
vector< ThermoPhase * > m_phase
Vector of the ThermoPhase pointers.
Definition: MultiPhase.h:571
size_t nPhases() const
Number of phases.
Definition: MultiPhase.h:425
double phaseMoles(const size_t n) const
Return the number of moles in phase n.
Definition: MultiPhase.cpp:776
void updatePhases() const
Set the states of the phase objects to the locally-stored state within this MultiPhase object.
Definition: MultiPhase.cpp:812
ThermoPhase & phase(size_t n)
Return a reference to phase n.
Definition: MultiPhase.cpp:149
vector< double > m_moles
Vector of the number of moles in each phase.
Definition: MultiPhase.h:568
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:231
int atomicNumber(size_t m) const
Atomic number of element m.
Definition: Phase.cpp:86
double temperature() const
Temperature (K).
Definition: Phase.h:562
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:142
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:434
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:580
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.
Definition: ThermoPhase.h:390
virtual double minTemp(size_t k=npos) const
Minimum temperature for which the thermodynamic data for the species or phase are valid.
Definition: ThermoPhase.h:451
virtual void setState_TPX(double t, double p, const double *x)
Set the temperature (K), pressure (Pa), and mole fractions.
Definition: ThermoPhase.cpp:85
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.
Definition: ThermoPhase.h:504
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.
Definition: stringUtils.cpp:58
virtual bool compatibleWithMultiPhase() const
Indicates whether this phase type can be used with class MultiPhase for equilibrium calculations.
Definition: ThermoPhase.h:1661
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:180
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
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
std::ostream & operator<<(std::ostream &s, MultiPhase &x)
Function to output a MultiPhase description to a stream.
Definition: MultiPhase.cpp:825
int _equilflag(const char *xy)
map property strings to integers
Definition: ChemEquil.cpp:20
map< string, double > Composition
Map from string names to doubles.
Definition: ct_defs.h:177
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.