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