Cantera  2.1.2
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  */
8 
12 #include "cantera/base/global.h"
13 
14 using namespace std;
15 
16 namespace Cantera
17 {
18 
19 MultiPhase::MultiPhase() :
20  m_np(0),
21  m_temp(298.15),
22  m_press(OneBar),
23  m_nel(0),
24  m_nsp(0),
25  m_init(false),
26  m_eloc(npos),
27  m_Tmin(1.0),
28  m_Tmax(100000.0)
29 {
30 }
31 
32 MultiPhase::MultiPhase(const MultiPhase& right) :
33  m_np(0),
34  m_temp(298.15),
35  m_press(OneBar),
36  m_nel(0),
37  m_nsp(0),
38  m_init(false),
39  m_eloc(npos),
40  m_Tmin(1.0),
41  m_Tmax(100000.0)
42 {
43  operator=(right);
44 }
45 
47 {
48 }
49 
50 MultiPhase& MultiPhase::operator=(const MultiPhase& right)
51 {
52  if (&right != this) {
53  m_moles = right.m_moles;
54  // shallow copy of phase pointers
55  m_phase = right.m_phase;
56  m_atoms = right.m_atoms;
58  m_spphase = right.m_spphase;
59  m_spstart = right.m_spstart;
60  m_enames = right.m_enames;
61  m_enamemap = right.m_enamemap;
62  m_np = right.m_np;
63  m_temp = right.m_temp;
64  m_press = right.m_press;
65  m_nel = right.m_nel;
66  m_nsp = right.m_nsp;
67  m_init = right.m_init;
68  m_eloc = right.m_eloc;
69  m_temp_OK = right.m_temp_OK;
70  m_Tmin = right.m_Tmin;
71  m_Tmax = right.m_Tmax;
73  }
74  return *this;
75 }
76 
77 void MultiPhase::
79 {
80  size_t n;
81  for (n = 0; n < mix.m_np; n++) {
82  addPhase(mix.m_phase[n], mix.m_moles[n]);
83  }
84 }
85 
86 void MultiPhase::
87 addPhases(std::vector<ThermoPhase*>& phases, const vector_fp& phaseMoles)
88 {
89  size_t np = phases.size();
90  size_t n;
91  for (n = 0; n < np; n++) {
92  addPhase(phases[n], phaseMoles[n]);
93  }
94  init();
95 }
96 
97 void MultiPhase::
98 addPhase(ThermoPhase* p, doublereal moles)
99 {
100  if (m_init) {
101  throw CanteraError("addPhase",
102  "phases cannot be added after init() has been called.");
103  }
104 
105  // save the pointer to the phase object
106  m_phase.push_back(p);
107 
108  // store its number of moles
109  m_moles.push_back(moles);
110  m_temp_OK.push_back(true);
111 
112  // update the number of phases and the total number of
113  // species
114  m_np = m_phase.size();
115  m_nsp += p->nSpecies();
116 
117  // determine if this phase has new elements
118  // for each new element, add an entry in the map
119  // from names to index number + 1:
120 
121  string ename;
122  // iterate over the elements in this phase
123  size_t m, nel = p->nElements();
124  for (m = 0; m < nel; m++) {
125  ename = p->elementName(m);
126 
127  // if no entry is found for this element name, then
128  // it is a new element. In this case, add the name
129  // to the list of names, increment the element count,
130  // and add an entry to the name->(index+1) map.
131  if (m_enamemap.find(ename) == m_enamemap.end()) {
132  m_enamemap[ename] = m_nel + 1;
133  m_enames.push_back(ename);
134  m_atomicNumber.push_back(p->atomicNumber(m));
135 
136  // Element 'E' (or 'e') is special. Note its location.
137  if (ename == "E" || ename == "e") {
138  m_eloc = m_nel;
139  }
140 
141  m_nel++;
142  }
143  }
144 
145  // If the mixture temperature hasn't been set, then set the
146  // temperature and pressure to the values for the phase being
147  // added. There is no good way to do this. However, this will be overridden later.
148  if (m_temp == 298.15 && p->temperature() > 2.0E-3) {
149  m_temp = p->temperature();
150  m_press = p->pressure();
151  }
152 
153  // If this is a solution phase, update the minimum and maximum
154  // mixture temperatures. Stoichiometric phases are excluded,
155  // since a mixture may define multiple stoichiometric phases,
156  // each of which has thermo data valid only over a limited
157  // range. For example, a mixture might be defined to contain a
158  // phase representing water ice and one representing liquid
159  // water, only one of which should be present if the mixture
160  // represents an equilibrium state.
161  if (p->nSpecies() > 1) {
162  double t = p->minTemp();
163  if (t > m_Tmin) {
164  m_Tmin = t;
165  }
166  t = p->maxTemp();
167  if (t < m_Tmax) {
168  m_Tmax = t;
169  }
170  }
171 }
172 
174 {
175  if (m_init) {
176  return;
177  }
178  size_t ip, kp, k = 0, nsp, m;
179  size_t mlocal;
180  string sym;
181 
182  // allocate space for the atomic composition matrix
183  m_atoms.resize(m_nel, m_nsp, 0.0);
184  m_moleFractions.resize(m_nsp, 0.0);
185  m_elemAbundances.resize(m_nel, 0.0);
186 
187  // iterate over the elements
188  // -> fill in m_atoms(m,k), m_snames(k), m_spphase(k),
189  // m_sptart(ip)
190  for (m = 0; m < m_nel; m++) {
191  sym = m_enames[m];
192  k = 0;
193  // iterate over the phases
194  for (ip = 0; ip < m_np; ip++) {
195  ThermoPhase* p = m_phase[ip];
196  nsp = p->nSpecies();
197  mlocal = p->elementIndex(sym);
198  for (kp = 0; kp < nsp; kp++) {
199  if (mlocal != npos) {
200  m_atoms(m, k) = p->nAtoms(kp, mlocal);
201  }
202  if (m == 0) {
203  m_snames.push_back(p->speciesName(kp));
204  if (kp == 0) {
205  m_spstart.push_back(m_spphase.size());
206  }
207  m_spphase.push_back(ip);
208  }
209  k++;
210  }
211  }
212  }
213 
214  if (m_eloc != npos) {
215  doublereal esum;
216  for (k = 0; k < m_nsp; k++) {
217  esum = 0.0;
218  for (m = 0; m < m_nel; m++) {
219  if (m != m_eloc) {
220  esum += m_atoms(m,k) * m_atomicNumber[m];
221  }
222  }
223  //m_atoms(m_eloc, k) += esum;
224  }
225  }
226 
227  /// set the initial composition within each phase to the
228  /// mole fractions stored in the phase objects
229  m_init = true;
230 
232 
233  updatePhases();
234 }
235 
237 {
238  if (!m_init) {
239  init();
240  }
241  m_phase[n]->setTemperature(m_temp);
242  m_phase[n]->setMoleFractions_NoNorm(DATA_PTR(m_moleFractions) + m_spstart[n]);
243  m_phase[n]->setPressure(m_press);
244  return *m_phase[n];
245 }
246 
247 void MultiPhase::checkPhaseIndex(size_t m) const
248 {
249  if (m >= nPhases()) {
250  throw IndexError("checkPhaseIndex", "phase", m, nPhases()-1);
251  }
252 }
253 
254 void MultiPhase::checkPhaseArraySize(size_t mm) const
255 {
256  if (nPhases() > mm) {
257  throw ArraySizeError("checkPhaseIndex", mm, nPhases());
258  }
259 }
260 
261 doublereal MultiPhase::speciesMoles(size_t k) const
262 {
263  size_t ip = m_spphase[k];
264  return m_moles[ip]*m_moleFractions[k];
265 }
266 
267 doublereal MultiPhase::elementMoles(size_t m) const
268 {
269  doublereal sum = 0.0, phasesum;
270  size_t i, k = 0, ik, nsp;
271  for (i = 0; i < m_np; i++) {
272  phasesum = 0.0;
273  nsp = m_phase[i]->nSpecies();
274  for (ik = 0; ik < nsp; ik++) {
275  k = speciesIndex(ik, i);
276  phasesum += m_atoms(m,k)*m_moleFractions[k];
277  }
278  sum += phasesum * m_moles[i];
279  }
280  return sum;
281 }
282 
283 doublereal MultiPhase::charge() const
284 {
285  doublereal sum = 0.0;
286  size_t i;
287  for (i = 0; i < m_np; i++) {
288  sum += phaseCharge(i);
289  }
290  return sum;
291 }
292 
293 size_t MultiPhase::speciesIndex(const std::string& speciesName, const std::string& phaseName)
294 {
295  if (!m_init) {
296  init();
297  }
298  size_t p = phaseIndex(phaseName);
299  if (p == npos) {
300  throw CanteraError("MultiPhase::speciesIndex", "phase not found: " + phaseName);
301  }
302  size_t k = m_phase[p]->speciesIndex(speciesName);
303  if (k == npos) {
304  throw CanteraError("MultiPhase::speciesIndex", "species not found: " + speciesName);
305  }
306  return m_spstart[p] + k;
307 }
308 
309 doublereal MultiPhase::phaseCharge(size_t p) const
310 {
311  doublereal phasesum = 0.0;
312  size_t ik, k, nsp = m_phase[p]->nSpecies();
313  for (ik = 0; ik < nsp; ik++) {
314  k = speciesIndex(ik, p);
315  phasesum += m_phase[p]->charge(ik)*m_moleFractions[k];
316  }
317  return Faraday*phasesum*m_moles[p];
318 }
319 
320 void MultiPhase::getChemPotentials(doublereal* mu) const
321 {
322  size_t i, loc = 0;
323  updatePhases();
324  for (i = 0; i < m_np; i++) {
325  m_phase[i]->getChemPotentials(mu + loc);
326  loc += m_phase[i]->nSpecies();
327  }
328 }
329 
330 void MultiPhase::getValidChemPotentials(doublereal not_mu,
331  doublereal* mu, bool standard) const
332 {
333  size_t i, loc = 0;
334 
335  updatePhases();
336  // iterate over the phases
337  for (i = 0; i < m_np; i++) {
338  if (tempOK(i) || m_phase[i]->nSpecies() > 1) {
339  if (!standard) {
340  m_phase[i]->getChemPotentials(mu + loc);
341  } else {
342  m_phase[i]->getStandardChemPotentials(mu + loc);
343  }
344  } else {
345  fill(mu + loc, mu + loc + m_phase[i]->nSpecies(), not_mu);
346  }
347  loc += m_phase[i]->nSpecies();
348  }
349 }
350 
351 bool MultiPhase::solutionSpecies(size_t k) const
352 {
353  if (m_phase[m_spphase[k]]->nSpecies() > 1) {
354  return true;
355  } else {
356  return false;
357  }
358 }
359 
360 doublereal MultiPhase::gibbs() const
361 {
362  size_t i;
363  doublereal sum = 0.0;
364  updatePhases();
365  for (i = 0; i < m_np; i++) {
366  if (m_moles[i] > 0.0) {
367  sum += m_phase[i]->gibbs_mole() * m_moles[i];
368  }
369  }
370  return sum;
371 }
372 
373 doublereal MultiPhase::enthalpy() const
374 {
375  size_t i;
376  doublereal sum = 0.0;
377  updatePhases();
378  for (i = 0; i < m_np; i++) {
379  if (m_moles[i] > 0.0) {
380  sum += m_phase[i]->enthalpy_mole() * m_moles[i];
381  }
382  }
383  return sum;
384 }
385 
386 doublereal MultiPhase::IntEnergy() const
387 {
388  size_t i;
389  doublereal sum = 0.0;
390  updatePhases();
391  for (i = 0; i < m_np; i++) {
392  if (m_moles[i] > 0.0) {
393  sum += m_phase[i]->intEnergy_mole() * m_moles[i];
394  }
395  }
396  return sum;
397 }
398 
399 doublereal MultiPhase::entropy() const
400 {
401  size_t i;
402  doublereal sum = 0.0;
403  updatePhases();
404  for (i = 0; i < m_np; i++) {
405  if (m_moles[i] > 0.0) {
406  sum += m_phase[i]->entropy_mole() * m_moles[i];
407  }
408  }
409  return sum;
410 }
411 
412 doublereal MultiPhase::cp() const
413 {
414  size_t i;
415  doublereal sum = 0.0;
416  updatePhases();
417  for (i = 0; i < m_np; i++) {
418  if (m_moles[i] > 0.0) {
419  sum += m_phase[i]->cp_mole() * m_moles[i];
420  }
421  }
422  return sum;
423 }
424 
425 void MultiPhase::setPhaseMoleFractions(const size_t n, const doublereal* const x)
426 {
427  if (!m_init) {
428  init();
429  }
430  ThermoPhase* p = m_phase[n];
431  p->setState_TPX(m_temp, m_press, x);
432  size_t istart = m_spstart[n];
433  for (size_t k = 0; k < p->nSpecies(); k++) {
434  m_moleFractions[istart+k] = x[k];
435  }
436 }
437 
439 {
440  size_t kk = nSpecies();
441  doublereal x;
442  vector_fp moles(kk, 0.0);
443  for (size_t k = 0; k < kk; k++) {
444  x = xMap[speciesName(k)];
445  if (x > 0.0) {
446  moles[k] = x;
447  }
448  }
449  setMoles(DATA_PTR(moles));
450 }
451 
452 void MultiPhase::setMolesByName(const std::string& x)
453 {
454  // build the composition map from the string, and then set the moles.
456  setMolesByName(xx);
457 }
458 
459 void MultiPhase::getMoles(doublereal* molNum) const
460 {
461  /*
462  * First copy in the mole fractions
463  */
464  copy(m_moleFractions.begin(), m_moleFractions.end(), molNum);
465  size_t ik;
466  doublereal* dtmp = molNum;
467  for (size_t ip = 0; ip < m_np; ip++) {
468  doublereal phasemoles = m_moles[ip];
469  ThermoPhase* p = m_phase[ip];
470  size_t nsp = p->nSpecies();
471  for (ik = 0; ik < nsp; ik++) {
472  *(dtmp++) *= phasemoles;
473  }
474  }
475 }
476 
477 void MultiPhase::setMoles(const doublereal* n)
478 {
479  if (!m_init) {
480  init();
481  }
482  size_t ip, loc = 0;
483  size_t ik, k = 0, nsp;
484  doublereal phasemoles;
485  for (ip = 0; ip < m_np; ip++) {
486  ThermoPhase* p = m_phase[ip];
487  nsp = p->nSpecies();
488  phasemoles = 0.0;
489  for (ik = 0; ik < nsp; ik++) {
490  phasemoles += n[k];
491  k++;
492  }
493  m_moles[ip] = phasemoles;
494  if (nsp > 1) {
495  if (phasemoles > 0.0) {
496  p->setState_TPX(m_temp, m_press, n + loc);
498  } else {
500  }
501  } else {
502  m_moleFractions[loc] = 1.0;
503  }
504  loc += nsp;
505  }
506 }
507 
508 void MultiPhase::addSpeciesMoles(const int indexS, const doublereal addedMoles)
509 {
510  vector_fp tmpMoles(m_nsp, 0.0);
511  getMoles(DATA_PTR(tmpMoles));
512  tmpMoles[indexS] += addedMoles;
513  if (tmpMoles[indexS] < 0.0) {
514  tmpMoles[indexS] = 0.0;
515  }
516  setMoles(DATA_PTR(tmpMoles));
517 }
518 
519 void MultiPhase::setState_TP(const doublereal T, const doublereal Pres)
520 {
521  if (!m_init) {
522  init();
523  }
524  m_temp = T;
525  m_press = Pres;
526  updatePhases();
527 }
528 
529 void MultiPhase::setState_TPMoles(const doublereal T, const doublereal Pres,
530  const doublereal* n)
531 {
532  m_temp = T;
533  m_press = Pres;
534  setMoles(n);
535 }
536 
537 void MultiPhase::getElemAbundances(doublereal* elemAbundances) const
538 {
539  size_t eGlobal;
541  for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
542  elemAbundances[eGlobal] = m_elemAbundances[eGlobal];
543  }
544 }
545 
547 {
548  size_t loc = 0;
549  size_t eGlobal;
550  size_t ik, kGlobal;
551  doublereal spMoles;
552  for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
553  m_elemAbundances[eGlobal] = 0.0;
554  }
555  for (size_t ip = 0; ip < m_np; ip++) {
556  ThermoPhase* p = m_phase[ip];
557  size_t nspPhase = p->nSpecies();
558  doublereal phasemoles = m_moles[ip];
559  for (ik = 0; ik < nspPhase; ik++) {
560  kGlobal = loc + ik;
561  spMoles = m_moleFractions[kGlobal] * phasemoles;
562  for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
563  m_elemAbundances[eGlobal] += m_atoms(eGlobal, kGlobal) * spMoles;
564  }
565  }
566  loc += nspPhase;
567  }
568 }
569 
570 doublereal MultiPhase::volume() const
571 {
572  int i;
573  doublereal sum = 0;
574  for (i = 0; i < int(m_np); i++) {
575  double vol = 1.0/m_phase[i]->molarDensity();
576  sum += m_moles[i] * vol;
577  }
578  return sum;
579 }
580 
581 doublereal MultiPhase::equilibrate(int XY, doublereal err,
582  int maxsteps, int maxiter, int loglevel)
583 {
584  bool strt = false;
585  doublereal dt;
586  doublereal h0;
587  int n;
588  doublereal hnow, herr = 1.0;
589  doublereal snow, serr = 1.0, s0;
590  doublereal Tlow = -1.0, Thigh = -1.0;
591  doublereal Hlow = Undef, Hhigh = Undef, tnew;
592  doublereal dta=0.0, dtmax, cpb;
593  MultiPhaseEquil* e = 0;
594 
595  if (!m_init) {
596  init();
597  }
598  if (loglevel > 0) {
599  beginLogGroup("MultiPhase::equilibrate", loglevel);
600  }
601 
602  if (XY == TP) {
603  if (loglevel > 0) {
604  addLogEntry("problem type","fixed T,P");
605  addLogEntry("Temperature",temperature());
606  addLogEntry("Pressure", pressure());
607  }
608 
609  // create an equilibrium manager
610  e = new MultiPhaseEquil(this);
611  try {
612  e->equilibrate(XY, err, maxsteps, loglevel);
613  } catch (CanteraError& err) {
614  err.save();
615  if (loglevel > 0) {
616  endLogGroup();
617  }
618  delete e;
619  e = 0;
620  throw err;
621  }
622  goto done;
623  }
624 
625  else if (XY == HP) {
626  h0 = enthalpy();
627  Tlow = 0.5*m_Tmin; // lower bound on T
628  Thigh = 2.0*m_Tmax; // upper bound on T
629  if (loglevel > 0) {
630  addLogEntry("problem type","fixed H,P");
631  addLogEntry("H target",fp2str(h0));
632  }
633  for (n = 0; n < maxiter; n++) {
634 
635  // if 'strt' is false, the current composition will be used as
636  // the starting estimate; otherwise it will be estimated
637  // if (e) {
638  // cout << "e should be zero, but it is not!" << endl;
639  // delete e;
640  // }
641  e = new MultiPhaseEquil(this, strt);
642  // start with a loose error tolerance, but tighten it as we get
643  // close to the final temperature
644  if (loglevel > 0) {
645  beginLogGroup("iteration "+int2str(n));
646  }
647 
648  try {
649  e->equilibrate(TP, err, maxsteps, loglevel);
650  hnow = enthalpy();
651  // the equilibrium enthalpy monotonically increases with T;
652  // if the current value is below the target, the we know the
653  // current temperature is too low. Set
654  if (hnow < h0) {
655  if (m_temp > Tlow) {
656  Tlow = m_temp;
657  Hlow = hnow;
658  }
659  }
660  // the current enthalpy is greater than the target; therefore the
661  // current temperature is too high.
662  else {
663  if (m_temp < Thigh) {
664  Thigh = m_temp;
665  Hhigh = hnow;
666  }
667  }
668  if (Hlow != Undef && Hhigh != Undef) {
669  cpb = (Hhigh - Hlow)/(Thigh - Tlow);
670  dt = (h0 - hnow)/cpb;
671  dta = fabs(dt);
672  dtmax = 0.5*fabs(Thigh - Tlow);
673  if (dta > dtmax) {
674  dt *= dtmax/dta;
675  }
676  } else {
677  tnew = sqrt(Tlow*Thigh);
678  dt = tnew - m_temp;
679  //cpb = cp();
680  }
681 
682  herr = fabs((h0 - hnow)/h0);
683  if (loglevel > 0) {
685  addLogEntry("H",fp2str(hnow));
686  addLogEntry("H rel error",fp2str(herr));
687  addLogEntry("lower T bound",fp2str(Tlow));
688  addLogEntry("upper T bound",fp2str(Thigh));
689  endLogGroup(); // iteration
690  }
691 
692 
693  if (herr < err) { // || dta < 1.0e-4) {
694  if (loglevel > 0) {
695  addLogEntry("T iterations",int2str(n));
696  addLogEntry("Final T",fp2str(temperature()));
697  addLogEntry("H rel error",fp2str(herr));
698  }
699  goto done;
700  }
701  tnew = m_temp + dt;
702  if (tnew < 0.0) {
703  tnew = 0.5*m_temp;
704  }
705  //dta = fabs(tnew - m_temp);
706  setTemperature(tnew);
707 
708  // if the size of Delta T is not too large, use
709  // the current composition as the starting estimate
710  if (dta < 100.0) {
711  strt = false;
712  }
713 
714  }
715 
716  catch (CanteraError& err) {
717  err.save();
718  if (!strt) {
719  if (loglevel > 0)
720  addLogEntry("no convergence",
721  "try estimating starting composition");
722  strt = true;
723  } else {
724  tnew = 0.5*(m_temp + Thigh);
725  if (fabs(tnew - m_temp) < 1.0) {
726  tnew = m_temp + 1.0;
727  }
728  setTemperature(tnew);
729  if (loglevel > 0)
730  addLogEntry("no convergence",
731  "trying T = "+fp2str(m_temp));
732  }
733  if (loglevel > 0) {
734  endLogGroup();
735  }
736  }
737  delete e;
738  e = 0;
739  }
740  if (loglevel > 0) {
741  addLogEntry("reached max number of T iterations",int2str(maxiter));
742  endLogGroup();
743  }
744  throw CanteraError("MultiPhase::equilibrate",
745  "No convergence for T");
746  } else if (XY == SP) {
747  s0 = entropy();
748  Tlow = 1.0; // m_Tmin; // lower bound on T
749  Thigh = 1.0e6; // m_Tmax; // upper bound on T
750  if (loglevel > 0) {
751  addLogEntry("problem type","fixed S,P");
752  addLogEntry("S target",fp2str(s0));
753  addLogEntry("min T",fp2str(Tlow));
754  addLogEntry("max T",fp2str(Thigh));
755  }
756  for (n = 0; n < maxiter; n++) {
757  delete e;
758  e = new MultiPhaseEquil(this, strt);
759  if (loglevel > 0) {
760  beginLogGroup("iteration "+int2str(n));
761  }
762 
763  try {
764  e->equilibrate(TP, err, maxsteps, loglevel);
765  snow = entropy();
766  if (snow < s0) {
767  if (m_temp > Tlow) {
768  Tlow = m_temp;
769  }
770  } else {
771  if (m_temp < Thigh) {
772  Thigh = m_temp;
773  }
774  }
775  serr = fabs((s0 - snow)/s0);
776  if (loglevel > 0) {
778  addLogEntry("S",fp2str(snow));
779  addLogEntry("S rel error",fp2str(serr));
780  endLogGroup();
781  }
782  dt = (s0 - snow)*m_temp/cp();
783  dtmax = 0.5*fabs(Thigh - Tlow);
784  dtmax = (dtmax > 500.0 ? 500.0 : dtmax);
785  dta = fabs(dt);
786  if (dta > dtmax) {
787  dt *= dtmax/dta;
788  }
789  if (herr < err || dta < 1.0e-4) {
790  if (loglevel > 0) {
791  addLogEntry("T iterations",int2str(n));
792  addLogEntry("Final T",fp2str(temperature()));
793  addLogEntry("S rel error",fp2str(serr));
794  }
795  goto done;
796  }
797  tnew = m_temp + dt;
798  setTemperature(tnew);
799 
800  // if the size of Delta T is not too large, use
801  // the current composition as the starting estimate
802  if (dta < 100.0) {
803  strt = false;
804  }
805  }
806 
807  catch (CanteraError& err) {
808  err.save();
809  if (!strt) {
810  if (loglevel > 0) {
811  addLogEntry("no convergence",
812  "setting strt to True");
813  }
814  strt = true;
815  } else {
816  tnew = 0.5*(m_temp + Thigh);
817  setTemperature(tnew);
818  if (loglevel > 0) {
819  addLogEntry("no convergence",
820  "trying T = "+fp2str(m_temp));
821  }
822  }
823  if (loglevel > 0) {
824  endLogGroup();
825  }
826  }
827  delete e;
828  e = 0;
829  }
830  if (loglevel > 0) {
831  addLogEntry("reached max number of T iterations",int2str(maxiter));
832  endLogGroup();
833  }
834  throw CanteraError("MultiPhase::equilibrate",
835  "No convergence for T");
836  } else if (XY == TV) {
837  addLogEntry("problem type","fixed T, V");
838  // doublereal dt = 1.0e3;
839  doublereal v0 = volume();
840  doublereal dVdP;
841  int n;
842  bool start = true;
843  doublereal vnow, pnow, verr;
844  for (n = 0; n < maxiter; n++) {
845  pnow = pressure();
846  MultiPhaseEquil e(this, start);
847  start = false;
848  beginLogGroup("iteration "+int2str(n));
849 
850  e.equilibrate(TP, err, maxsteps, loglevel);
851  vnow = volume();
852  verr = fabs((v0 - vnow)/v0);
853  addLogEntry("P",fp2str(pressure()));
854  addLogEntry("V rel error",fp2str(verr));
855  endLogGroup();
856 
857  if (verr < err) {
858  addLogEntry("P iterations",int2str(n));
859  addLogEntry("Final P",fp2str(pressure()));
860  addLogEntry("V rel error",fp2str(verr));
861  goto done;
862  }
863  // find dV/dP
864  setPressure(pnow*1.01);
865  dVdP = (volume() - vnow)/(0.01*pnow);
866  setPressure(pnow + 0.5*(v0 - vnow)/dVdP);
867  }
868  }
869 
870  else {
871  if (loglevel > 0) {
872  endLogGroup();
873  }
874  throw CanteraError("MultiPhase::equilibrate","unknown option");
875  }
876  return -1.0;
877 done:
878  delete e;
879  e = 0;
880  if (loglevel > 0) {
881  endLogGroup();
882  }
883  return err;
884 }
885 
886 #ifdef MULTIPHASE_DEVEL
887 void importFromXML(string infile, string id)
888 {
889  XML_Node* root = get_XML_File(infile);
890  if (id == "-") {
891  id = "";
892  }
893  XML_Node* x = get_XML_Node(string("#")+id, root);
894  if (x.name() != "multiphase")
895  throw CanteraError("MultiPhase::importFromXML",
896  "Current XML_Node is not a multiphase element.");
897  vector<XML_Node*> phases;
898  x.getChildren("phase",phases);
899  int np = phases.size();
900  int n;
901  ThermoPhase* p;
902  for (n = 0; n < np; n++) {
903  XML_Node& ph = *phases[n];
904  srcfile = infile;
905  if (ph.hasAttrib("src")) {
906  srcfile = ph["src"];
907  }
908  idstr = ph["id"];
909  p = newPhase(srcfile, idstr);
910  if (p) {
911  addPhase(p, ph.value());
912  }
913  }
914 }
915 #endif
916 
917 void MultiPhase::setTemperature(const doublereal T)
918 {
919  if (!m_init) {
920  init();
921  }
922  m_temp = T;
923  updatePhases();
924 }
925 
926 void MultiPhase::checkElementIndex(size_t m) const
927 {
928  if (m >= m_nel) {
929  throw IndexError("checkElementIndex", "elements", m, m_nel-1);
930  }
931 }
932 
934 {
935  if (m_nel > mm) {
936  throw ArraySizeError("checkElementArraySize", mm, m_nel);
937  }
938 }
939 
940 std::string MultiPhase::elementName(size_t m) const
941 {
942  return m_enames[m];
943 }
944 
945 size_t MultiPhase::elementIndex(const std::string& name) const
946 {
947  for (size_t e = 0; e < m_nel; e++) {
948  if (m_enames[e] == name) {
949  return e;
950  }
951  }
952  return npos;
953 }
954 
955 void MultiPhase::checkSpeciesIndex(size_t k) const
956 {
957  if (k >= m_nsp) {
958  throw IndexError("checkSpeciesIndex", "species", k, m_nsp-1);
959  }
960 }
961 
963 {
964  if (m_nsp > kk) {
965  throw ArraySizeError("checkSpeciesArraySize", kk, m_nsp);
966  }
967 }
968 
969 std::string MultiPhase::speciesName(const size_t k) const
970 {
971  return m_snames[k];
972 }
973 
974 doublereal MultiPhase::nAtoms(const size_t kGlob, const size_t mGlob) const
975 {
976  return m_atoms(mGlob, kGlob);
977 }
978 
979 void MultiPhase::getMoleFractions(doublereal* const x) const
980 {
981  std::copy(m_moleFractions.begin(), m_moleFractions.end(), x);
982 }
983 
984 std::string MultiPhase::phaseName(const size_t iph) const
985 {
986  const ThermoPhase* tptr = m_phase[iph];
987  return tptr->id();
988 }
989 
990 int MultiPhase::phaseIndex(const std::string& pName) const
991 {
992  std::string tmp;
993  for (int iph = 0; iph < (int) m_np; iph++) {
994  const ThermoPhase* tptr = m_phase[iph];
995  tmp = tptr->id();
996  if (tmp == pName) {
997  return iph;
998  }
999  }
1000  return -1;
1001 }
1002 
1003 doublereal MultiPhase::phaseMoles(const size_t n) const
1004 {
1005  return m_moles[n];
1006 }
1007 
1008 void MultiPhase::setPhaseMoles(const size_t n, const doublereal moles)
1009 {
1010  m_moles[n] = moles;
1011 }
1012 
1013 size_t MultiPhase::speciesPhaseIndex(const size_t kGlob) const
1014 {
1015  return m_spphase[kGlob];
1016 }
1017 
1018 doublereal MultiPhase::moleFraction(const size_t kGlob) const
1019 {
1020  return m_moleFractions[kGlob];
1021 }
1022 
1023 bool MultiPhase::tempOK(const size_t p) const
1024 {
1025  return m_temp_OK[p];
1026 }
1027 
1029 {
1030  size_t ip, loc = 0;
1031  for (ip = 0; ip < m_np; ip++) {
1032  ThermoPhase* p = m_phase[ip];
1034  loc += p->nSpecies();
1035  }
1037 }
1038 
1040 {
1041  size_t p, nsp, loc = 0;
1042  for (p = 0; p < m_np; p++) {
1043  nsp = m_phase[p]->nSpecies();
1044  const doublereal* x = DATA_PTR(m_moleFractions) + loc;
1045  loc += nsp;
1046  m_phase[p]->setState_TPX(m_temp, m_press, x);
1047  m_temp_OK[p] = true;
1048  if (m_temp < m_phase[p]->minTemp()
1049  || m_temp > m_phase[p]->maxTemp()) {
1050  m_temp_OK[p] = false;
1051  }
1052  }
1053 }
1054 }
std::map< std::string, doublereal > compositionMap
Map connecting a string name with a double.
Definition: ct_defs.h:162
XML_Node * get_XML_Node(const std::string &file_ID, XML_Node *root)
This routine will locate an XML node in either the input XML tree or in another input file specified ...
Definition: global.cpp:249
doublereal nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
Definition: Phase.cpp:215
ThermoPhase * newPhase(XML_Node &xmlphase)
Create a new ThermoPhase object and initializes it according to the XML tree.
Array size error.
Definition: ctexceptions.h:123
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
std::vector< size_t > m_spphase
Mapping between the global species number and the phase ID.
Definition: MultiPhase.h:569
doublereal m_temp
Current value of the temperature (kelvin)
Definition: MultiPhase.h:608
XML_Node * get_XML_File(const std::string &file, int debug)
Return a pointer to the XML tree for a Cantera input file.
Definition: global.cpp:144
doublereal equilibrate(int XY, doublereal err=1.0e-9, int maxsteps=1000, int maxiter=200, int loglevel=-99)
Set the mixture to a state of chemical equilibrium.
Definition: MultiPhase.cpp:581
size_t nElements() const
Number of elements.
Definition: Phase.cpp:139
bool tempOK(size_t p) const
Return true if the phase p has valid thermo data for the current temperature.
void beginLogGroup(const std::string &title, int loglevel)
Create a new group for log messages.
Definition: global.cpp:90
MultiPhase()
Constructor.
Definition: MultiPhase.cpp:19
bool m_init
True if the init() routine has been called, and the MultiPhase frozen.
Definition: MultiPhase.h:620
virtual void setState_TPX(doublereal t, doublereal p, const doublereal *x)
Set the temperature (K), pressure (Pa), and mole fractions.
void setPhaseMoles(const size_t n, const doublereal moles)
Set the number of moles of phase with index n.
doublereal speciesMoles(size_t kGlob) const
Returns the moles of global species k. units = kmol.
Definition: MultiPhase.cpp:261
void checkElementIndex(size_t m) const
Check that the specified element index is in range.
Definition: MultiPhase.cpp:926
std::map< std::string, size_t > m_enamemap
Returns the global element index, given the element string name.
Definition: MultiPhase.h:602
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:173
std::string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
Definition: MultiPhase.cpp:969
void checkSpeciesArraySize(size_t kk) const
Check that an array size is at least nSpecies().
Definition: MultiPhase.cpp:962
vector_int m_atomicNumber
Atomic number of each element.
Definition: MultiPhase.h:589
doublereal phaseMoles(const size_t n) const
Return the number of moles in phase n.
doublereal temperature() const
Temperature [K].
Definition: MultiPhase.h:335
void addPhase(ThermoPhase *p, doublereal moles)
Add a phase to the mixture.
Definition: MultiPhase.cpp:98
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:100
doublereal phaseCharge(size_t p) const
Charge (Coulombs) of phase with index p.
Definition: MultiPhase.cpp:309
int phaseIndex(const std::string &pName) const
Returns the index, given the phase name.
Definition: MultiPhase.cpp:990
This file contains definitions for utility functions and text for modules, inputfiles, logs, textlogs, HTML_logs (see Input File Handling, Diagnostic Output, Writing messages to the screen and Writing HTML Logfiles).
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:519
Headers for the MultiPhase object that is used to set up multiphase equilibrium problems (see Equilib...
void getChemPotentials(doublereal *mu) const
Returns a vector of Chemical potentials.
Definition: MultiPhase.cpp:320
virtual ~MultiPhase()
Destructor.
Definition: MultiPhase.cpp:46
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:254
const doublereal Undef
Fairly random number to be used to initialize variables against to see if they are subsequently defin...
Definition: ct_defs.h:147
vector_fp m_moles
Vector of the number of moles in each phase.
Definition: MultiPhase.h:545
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:101
void getMoleFractions(doublereal *const x) const
Returns the global Species mole fractions.
Definition: MultiPhase.cpp:979
size_t speciesPhaseIndex(const size_t kGlob) const
Returns the phase index of the Kth "global" species.
std::vector< std::string > m_enames
String names of the global elements.
Definition: MultiPhase.h:583
doublereal gibbs() const
The Gibbs function of the mixture [J].
Definition: MultiPhase.cpp:360
void addSpeciesMoles(const int indexS, const doublereal addedMoles)
Adds moles of a certain species to the mixture.
Definition: MultiPhase.cpp:508
vector_fp m_moleFractions
Locally stored vector of mole fractions of all species comprising the MultiPhase object.
Definition: MultiPhase.h:562
void uploadMoleFractionsFromPhases()
Update the locally-stored composition within this object to match the current compositions of the pha...
A class for multiphase mixtures.
Definition: MultiPhase.h:54
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
Definition: DenseMatrix.cpp:71
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:247
doublereal charge() const
Total charge summed over all phases (Coulombs).
Definition: MultiPhase.cpp:283
std::vector< bool > m_temp_OK
Vector of bools indicating whether temperatures are ok for phases.
Definition: MultiPhase.h:633
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:529
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range.
Definition: MultiPhase.cpp:955
std::string elementName(size_t m) const
Returns the name of the global element m.
Definition: MultiPhase.cpp:940
size_t m_eloc
Global ID of the element corresponding to the electronic charge.
Definition: MultiPhase.h:626
void getValidChemPotentials(doublereal not_mu, doublereal *mu, bool standard=false) const
Returns a vector of Valid chemical potentials.
Definition: MultiPhase.cpp:330
doublereal entropy() const
The entropy of the mixture [J/K].
Definition: MultiPhase.cpp:399
virtual doublereal maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
Definition: ThermoPhase.h:244
void addPhases(std::vector< ThermoPhase * > &phases, const vector_fp &phaseMoles)
Add a vector of phases to the mixture.
Definition: MultiPhase.cpp:87
std::string fp2str(const double x, const std::string &fmt)
Convert a double into a c++ string.
Definition: stringUtils.cpp:29
DenseMatrix m_atoms
Global Stoichiometric Coefficient array.
Definition: MultiPhase.h:556
int atomicNumber(size_t m) const
Atomic number of element m.
Definition: Phase.cpp:198
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:577
void setPhaseMoleFractions(const size_t n, const doublereal *const x)
Set the Mole fractions of the nth phase.
Definition: MultiPhase.cpp:425
size_t m_nsp
Number of distinct species in all of the phases.
Definition: MultiPhase.h:617
void calcElemAbundances() const
Calculate the element abundance vector.
Definition: MultiPhase.cpp:546
std::string id() const
Return the string id for the phase.
Definition: Phase.cpp:119
size_t m_nel
Number of distinct elements in all of the phases.
Definition: MultiPhase.h:614
void setTemperature(const doublereal T)
Set the temperature [K].
Definition: MultiPhase.cpp:917
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
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.
void checkElementArraySize(size_t mm) const
Check that an array size is at least nElements().
Definition: MultiPhase.cpp:933
size_t elementIndex(const std::string &name) const
Returns the index of the element with name name.
Definition: MultiPhase.cpp:945
doublereal IntEnergy() const
The internal energy of the mixture [J].
Definition: MultiPhase.cpp:386
size_t m_np
Number of phases in the MultiPhase object.
Definition: MultiPhase.h:605
vector_fp m_elemAbundances
Vector of element abundances.
Definition: MultiPhase.h:653
std::vector< ThermoPhase * > m_phase
Vector of the ThermoPhase pointers.
Definition: MultiPhase.h:548
void getMoles(doublereal *molNum) const
Get the mole numbers of all species in the multiphase object.
Definition: MultiPhase.cpp:459
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:974
void setPressure(doublereal P)
Set the pressure [Pa].
Definition: MultiPhase.h:393
doublereal moleFraction(const size_t kGlob) const
Returns the mole fraction of global species k.
void endLogGroup(const std::string &title)
Close the current group of log messages.
Definition: global.cpp:115
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:252
void updatePhases() const
Set the states of the phase objects to the locally-stored state within this MultiPhase object...
void setState_TP(const doublereal T, const doublereal Pres)
Set the state of the underlying ThermoPhase objects in one call.
Definition: MultiPhase.cpp:519
void addLogEntry(const std::string &tag, const std::string &value)
Add an entry to an HTML log file.
Definition: global.cpp:95
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
Definition: ThermoPhase.h:314
void setMolesByName(compositionMap &xMap)
Set the number of moles of species in the mixture.
Definition: MultiPhase.cpp:438
doublereal temperature() const
Temperature (K).
Definition: Phase.h:528
doublereal enthalpy() const
The enthalpy of the mixture [J].
Definition: MultiPhase.cpp:373
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:165
size_t elementIndex(const std::string &name) const
Return the index of element named 'name'.
Definition: Phase.cpp:164
doublereal maxTemp() const
Maximum temperature for which all solution phases have valid thermo data.
Definition: MultiPhase.h:258
doublereal cp() const
Heat capacity at constant pressure [J/K].
Definition: MultiPhase.cpp:412
Headers for the DenseMatrix object, which deals with dense rectangular matrices and description of th...
std::string phaseName(const size_t iph) const
Returns the name of the n'th phase.
Definition: MultiPhase.cpp:984
doublereal elementMoles(size_t m) const
Total moles of global element m, summed over all phases.
Definition: MultiPhase.cpp:267
void init()
Process phases and build atomic composition array.
Definition: MultiPhase.cpp:173
Contains declarations for string manipulation functions within Cantera.
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
std::vector< std::string > m_snames
Vector of species names in the problem.
Definition: MultiPhase.h:596
std::string elementName(size_t m) const
Name of the element with index m.
Definition: Phase.cpp:158
An array index is out of range.
Definition: ctexceptions.h:153
bool solutionSpecies(size_t kGlob) const
Return true is species kGlob is a species in a multicomponent solution phase.
Definition: MultiPhase.cpp:351
thermo_t & phase(size_t n)
Return a reference to phase n.
Definition: MultiPhase.cpp:236
doublereal minTemp() const
Minimum temperature for which all solution phases have valid thermo data.
Definition: MultiPhase.h:251
doublereal m_press
Current value of the pressure (Pa)
Definition: MultiPhase.h:611
doublereal pressure() const
Pressure [Pa].
Definition: MultiPhase.h:378
void getElemAbundances(doublereal *elemAbundances) const
Retrieves a vector of element abundances.
Definition: MultiPhase.cpp:537
size_t nSpecies() const
Number of species, summed over all phases.
Definition: MultiPhase.h:125
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:231
doublereal m_Tmax
Minimum temperature for which thermo parameterizations are valid.
Definition: MultiPhase.h:646
Header file for class ThermoPhase, the base class for phases with thermodynamic properties, and the text for the Module thermoprops (see Thermodynamic Properties and class ThermoPhase).
doublereal m_Tmin
Minimum temperature for which thermo parameterizations are valid.
Definition: MultiPhase.h:640
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:246
void setMoles(const doublereal *n)
Sets all of the global species mole numbers.
Definition: MultiPhase.cpp:477
doublereal volume() const
The total mixture volume [m^3].
Definition: MultiPhase.cpp:570
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:175
size_t nPhases() const
Number of phases.
Definition: MultiPhase.h:415
void save()
Function to put this error onto Cantera's error stack.