Cantera  2.0
ThermoPhase.cpp
Go to the documentation of this file.
1 /**
2  * @file ThermoPhase.cpp
3  * Definition file for class ThermoPhase, the base class for phases with
4  * thermodynamic properties
5  * (see class \link Cantera::ThermoPhase ThermoPhase\endlink).
6  */
7 
8 // Copyright 2002 California Institute of Technology
9 
11 #include "cantera/base/mdp_allo.h"
13 
14 #include <iomanip>
15 #include <fstream>
16 
17 using namespace std;
18 using namespace ctml;
19 
20 namespace Cantera
21 {
22 
23 //! Constructor. Note that ThermoPhase is meant to be used as
24 //! a base class, so this constructor should not be called
25 //! explicitly.
26 ThermoPhase::ThermoPhase() :
27  Phase(),
28  m_spthermo(0), m_speciesData(0),
29  m_phi(0.0),
30  m_hasElementPotentials(false),
31  m_chargeNeutralityNecessary(false),
32  m_ssConvention(cSS_CONVENTION_TEMPERATURE)
33 {
34 }
35 
37 {
38  for (size_t k = 0; k < m_kk; k++) {
39  if (m_speciesData[k]) {
40  delete m_speciesData[k];
41  m_speciesData[k] = 0;
42  }
43  }
44  delete m_spthermo;
45  m_spthermo = 0;
46 }
47 
48 //====================================================================================================================
49 /*
50  * Copy Constructor for the ThermoPhase object.
51  *
52  * Currently, this is implemented, but not tested. If called it will
53  * throw an exception until fully tested.
54  */
56  Phase(),
57  m_spthermo(0),
58  m_speciesData(0),
59  m_phi(0.0),
60  m_hasElementPotentials(false),
61  m_chargeNeutralityNecessary(false),
62  m_ssConvention(cSS_CONVENTION_TEMPERATURE)
63 {
64  /*
65  * Call the assignment operator
66  */
67  *this = operator=(right);
68 }
69 //====================================================================================================================
70 /*
71  * operator=()
72  *
73  * Note this stuff will not work until the underlying phase
74  * has a working assignment operator
75  */
77 operator=(const ThermoPhase& right)
78 {
79  /*
80  * Check for self assignment.
81  */
82  if (this == &right) {
83  return *this;
84  }
85 
86  /*
87  * We need to destruct first
88  */
89  for (size_t k = 0; k < m_kk; k++) {
90  if (m_speciesData[k]) {
91  delete m_speciesData[k];
92  m_speciesData[k] = 0;
93  }
94  }
95  if (m_spthermo) {
96  delete m_spthermo;
97  }
98 
99  /*
100  * Call the base class assignment operator
101  */
102  (void)Phase::operator=(right);
103 
104  /*
105  * Pointer to the species thermodynamic property manager
106  * We own this, so we need to do a deep copy
107  */
108  m_spthermo = (right.m_spthermo)->duplMyselfAsSpeciesThermo();
109 
110  /*
111  * Do a deep copy of species Data, because we own this
112  */
113  m_speciesData.resize(m_kk);
114  for (size_t k = 0; k < m_kk; k++) {
115  m_speciesData[k] = new XML_Node(*(right.m_speciesData[k]));
116  }
117 
118  m_phi = right.m_phi;
119  m_lambdaRRT = right.m_lambdaRRT;
123  return *this;
124 }
125 //====================================================================================================================
126 /*
127  * Duplication routine for objects which inherit from
128  * ThermoPhase.
129  *
130  * This virtual routine can be used to duplicate thermophase objects
131  * inherited from ThermoPhase even if the application only has
132  * a pointer to ThermoPhase to work with.
133  *
134  * Currently, this is not fully implemented. If called, an
135  * exception will be called by the ThermoPhase copy constructor.
136  */
138 {
139  ThermoPhase* tp = new ThermoPhase(*this);
140  return tp;
141 }
142 //====================================================================================================================
144 {
145  return cAC_CONVENTION_MOLAR;
146 }
147 //=================================================================================================================
149 {
150  return m_ssConvention;
151 }
152 //=================================================================================================================
153 doublereal ThermoPhase::logStandardConc(size_t k) const
154 {
155  return log(standardConcentration(k));
156 }
157 //=================================================================================================================
158 void ThermoPhase::getActivities(doublereal* a) const
159 {
161  for (size_t k = 0; k < nSpecies(); k++) {
162  a[k] /= standardConcentration(k);
163  }
164 }
165 //=================================================================================================================
166 void ThermoPhase::getLnActivityCoefficients(doublereal* lnac) const
167 {
169  for (size_t k = 0; k < m_kk; k++) {
170  lnac[k] = std::log(lnac[k]);
171  }
172 }
173 //=================================================================================================================
174 void ThermoPhase::setState_TPX(doublereal t, doublereal p, const doublereal* x)
175 {
176  setMoleFractions(x);
177  setTemperature(t);
178  setPressure(p);
179 }
180 //=================================================================================================================
181 void ThermoPhase::setState_TPX(doublereal t, doublereal p, compositionMap& x)
182 {
184  setTemperature(t);
185  setPressure(p);
186 }
187 //=================================================================================================================
188 void ThermoPhase::setState_TPX(doublereal t, doublereal p, const std::string& x)
189 {
190  compositionMap xx;
191  for (size_t k = 0; k < nSpecies(); k++) {
192  xx[speciesName(k)] = -1.0;
193  }
194  try {
195  parseCompString(x, xx);
196  } catch (CanteraError& err) {
197  err.save();
198  throw CanteraError("setState_TPX",
199  "Unknown species in composition map: "+ x);
200  }
202  setTemperature(t);
203  setPressure(p);
204 }
205 //=================================================================================================================
206 void ThermoPhase::setState_TPY(doublereal t, doublereal p,
207  const doublereal* y)
208 {
209  setMassFractions(y);
210  setTemperature(t);
211  setPressure(p);
212 }
213 //=================================================================================================================
214 void ThermoPhase::setState_TPY(doublereal t, doublereal p,
215  compositionMap& y)
216 {
218  setTemperature(t);
219  setPressure(p);
220 }
221 //=================================================================================================================
222 void ThermoPhase::setState_TPY(doublereal t, doublereal p,
223  const std::string& y)
224 {
225  compositionMap yy;
226  for (size_t k = 0; k < nSpecies(); k++) {
227  yy[speciesName(k)] = -1.0;
228  }
229  try {
230  parseCompString(y, yy);
231  } catch (CanteraError& err) {
232  err.save();
233  throw CanteraError("setState_TPY",
234  "Unknown species in composition map: "+ y);
235  }
237  setTemperature(t);
238  setPressure(p);
239 }
240 //=================================================================================================================
241 
242 void ThermoPhase::setState_TP(doublereal t, doublereal p)
243 {
244  setTemperature(t);
245  setPressure(p);
246 }
247 //=================================================================================================================
248 
249 void ThermoPhase::setState_PX(doublereal p, doublereal* x)
250 {
251  setMoleFractions(x);
252  setPressure(p);
253 }
254 //=================================================================================================================
255 
256 void ThermoPhase::setState_PY(doublereal p, doublereal* y)
257 {
258  setMassFractions(y);
259  setPressure(p);
260 }
261 //=================================================================================================================
262 
263 void ThermoPhase::setState_HP(doublereal Htarget, doublereal p,
264  doublereal dTtol)
265 {
266  setState_HPorUV(Htarget, p, dTtol, false);
267 }
268 //=================================================================================================================
269 
270 void ThermoPhase::setState_UV(doublereal u, doublereal v,
271  doublereal dTtol)
272 {
273  setState_HPorUV(u, v, dTtol, true);
274 }
275 //=================================================================================================================
276 
277 // Do the convergence work
278 /*
279  * We assume here that H at constant P is a monotonically increasing
280  * function of T.
281  * We assume here that U at constant V is a monotonically increasing
282  * function of T.
283  *
284  * Note, the value of dTtol may become important for some applications
285  * where numerical jacobians are being calculated.
286  */
287 void ThermoPhase::setState_HPorUV(doublereal Htarget, doublereal p,
288  doublereal dTtol, bool doUV)
289 {
290  doublereal dt;
291  doublereal Hmax = 0.0, Hmin = 0.0;
292  doublereal v = 0.0;
293 
294  // Assign the specific volume or pressure and make sure it's positive
295  if (doUV) {
296  v = p;
297  if (v < 1.0E-300) {
298  throw CanteraError("setState_HPorUV (UV)",
299  "Input specific volume is too small or negative. v = " + fp2str(v));
300  }
301  setDensity(1.0/v);
302  } else {
303  if (p < 1.0E-300) {
304  throw CanteraError("setState_HPorUV (HP)",
305  "Input pressure is too small or negative. p = " + fp2str(p));
306  }
307  setPressure(p);
308  }
309  double Tmax = maxTemp() + 0.1;
310  double Tmin = minTemp() - 0.1;
311 
312  // Make sure we are within the temperature bounds at the start
313  // of the iteration
314  double Tnew = temperature();
315  double Tinit = Tnew;
316  if (Tnew > Tmax) {
317  Tnew = Tmax - 1.0;
318  if (doUV) {
319  setTemperature(Tnew);
320  } else {
321  setState_TP(Tnew, p);
322  }
323  }
324  if (Tnew < Tmin) {
325  Tnew = Tmin + 1.0;
326  if (doUV) {
327  setTemperature(Tnew);
328  } else {
329  setState_TP(Tnew, p);
330  }
331  }
332 
333  double Hnew = 0.0;
334  double Cpnew = 0.0;
335  if (doUV) {
336  Hnew = intEnergy_mass();
337  Cpnew = cv_mass();
338  } else {
339  Hnew = enthalpy_mass();
340  Cpnew = cp_mass();
341  }
342  double Htop = Hnew;
343  double Ttop = Tnew;
344  double Hbot = Hnew;
345  double Tbot = Tnew;
346  double Told = Tnew;
347  double Hold = Hnew;
348 
349  bool ignoreBounds = false;
350  // Unstable phases are those for which
351  // cp < 0.0. These are possible for cases where
352  // we have passed the spinodal curve.
353  bool unstablePhase = false;
354  // Counter indicating the last temperature point where the
355  // phase was unstable
356  double Tunstable = -1.0;
357  bool unstablePhaseNew = false;
358 
359 
360  // Newton iteration
361  for (int n = 0; n < 500; n++) {
362  Told = Tnew;
363  Hold = Hnew;
364  double cpd = Cpnew;
365  if (cpd < 0.0) {
366  unstablePhase = true;
367  Tunstable = Tnew;
368  }
369  dt = (Htarget - Hold)/cpd;
370 
371  // limit step size to 100 K
372  if (dt > 100.0) {
373  dt = 100.0;
374  } else if (dt < -100.0) {
375  dt = -100.0;
376  }
377 
378  // Calculate the new T
379  Tnew = Told + dt;
380 
381  // Limit the step size so that we are convergent
382  // This is the step that makes it different from a
383  // Newton's algorithm
384  if (dt > 0.0) {
385  if (!unstablePhase) {
386  if (Htop > Htarget) {
387  if (Tnew > (0.75 * Ttop + 0.25 * Told)) {
388  dt = 0.75 * (Ttop - Told);
389  Tnew = Told + dt;
390  }
391  }
392  } else {
393  if (Hbot < Htarget) {
394  if (Tnew < (0.75 * Tbot + 0.25 * Told)) {
395  dt = 0.75 * (Tbot - Told);
396  Tnew = Told + dt;
397  }
398  }
399  }
400  } else {
401  if (!unstablePhase) {
402  if (Hbot < Htarget) {
403  if (Tnew < (0.75 * Tbot + 0.25 * Told)) {
404  dt = 0.75 * (Tbot - Told);
405  Tnew = Told + dt;
406  }
407  }
408  } else {
409  if (Htop > Htarget) {
410  if (Tnew > (0.75 * Ttop + 0.25 * Told)) {
411  dt = 0.75 * (Ttop - Told);
412  Tnew = Told + dt;
413  }
414  }
415  }
416  }
417  // Check Max and Min values
418  if (Tnew > Tmax) {
419  if (!ignoreBounds) {
420  if (doUV) {
421  setTemperature(Tmax);
422  Hmax = intEnergy_mass();
423  } else {
424  setState_TP(Tmax, p);
425  Hmax = enthalpy_mass();
426  }
427  if (Hmax >= Htarget) {
428  if (Htop < Htarget) {
429  Ttop = Tmax;
430  Htop = Hmax;
431  }
432  } else {
433  Tnew = Tmax + 1.0;
434  ignoreBounds = true;
435  }
436  }
437  }
438  if (Tnew < Tmin) {
439  if (!ignoreBounds) {
440  if (doUV) {
441  setTemperature(Tmin);
442  Hmin = intEnergy_mass();
443  } else {
444  setState_TP(Tmin, p);
445  Hmin = enthalpy_mass();
446  }
447  if (Hmin <= Htarget) {
448  if (Hbot > Htarget) {
449  Tbot = Tmin;
450  Hbot = Hmin;
451  }
452  } else {
453  Tnew = Tmin - 1.0;
454  ignoreBounds = true;
455  }
456  }
457  }
458 
459  // Try to keep phase within its region of stability
460  // -> Could do a lot better if I calculate the
461  // spinodal value of H.
462  for (int its = 0; its < 10; its++) {
463  Tnew = Told + dt;
464  if (doUV) {
465  setTemperature(Tnew);
466  Hnew = intEnergy_mass();
467  Cpnew = cv_mass();
468  } else {
469  setState_TP(Tnew, p);
470  Hnew = enthalpy_mass();
471  Cpnew = cp_mass();
472  }
473  if (Cpnew < 0.0) {
474  unstablePhaseNew = true;
475  Tunstable = Tnew;
476  } else {
477  unstablePhaseNew = false;
478  break;
479  }
480  if (unstablePhase == false) {
481  if (unstablePhaseNew == true) {
482  dt *= 0.25;
483  }
484  }
485  }
486 
487  if (Hnew == Htarget) {
488  return;
489  } else if (Hnew > Htarget) {
490  if ((Htop < Htarget) || (Hnew < Htop)) {
491  Htop = Hnew;
492  Ttop = Tnew;
493  }
494  } else if (Hnew < Htarget) {
495  if ((Hbot > Htarget) || (Hnew > Hbot)) {
496  Hbot = Hnew;
497  Tbot = Tnew;
498  }
499  }
500  // Convergence in H
501  double Herr = Htarget - Hnew;
502  double acpd = std::max(fabs(cpd), 1.0E-5);
503  double denom = std::max(fabs(Htarget), acpd * dTtol);
504  double HConvErr = fabs((Herr)/denom);
505  if (HConvErr < 0.00001 *dTtol) {
506  return;
507  }
508  if (fabs(dt) < dTtol) {
509  return;
510  }
511 
512  }
513  // We are here when there hasn't been convergence
514  /*
515  * Formulate a detailed error message, since questions seem to
516  * arise often about the lack of convergence.
517  */
518  string ErrString = "No convergence in 500 iterations\n";
519  if (doUV) {
520  ErrString += "\tTarget Internal Energy = " + fp2str(Htarget) + "\n";
521  ErrString += "\tCurrent Specific Volume = " + fp2str(v) + "\n";
522  ErrString += "\tStarting Temperature = " + fp2str(Tinit) + "\n";
523  ErrString += "\tCurrent Temperature = " + fp2str(Tnew) + "\n";
524  ErrString += "\tCurrent Internal Energy = " + fp2str(Hnew) + "\n";
525  ErrString += "\tCurrent Delta T = " + fp2str(dt) + "\n";
526  } else {
527  ErrString += "\tTarget Enthalpy = " + fp2str(Htarget) + "\n";
528  ErrString += "\tCurrent Pressure = " + fp2str(p) + "\n";
529  ErrString += "\tStarting Temperature = " + fp2str(Tinit) + "\n";
530  ErrString += "\tCurrent Temperature = " + fp2str(Tnew) + "\n";
531  ErrString += "\tCurrent Enthalpy = " + fp2str(Hnew) + "\n";
532  ErrString += "\tCurrent Delta T = " + fp2str(dt) + "\n";
533  }
534  if (unstablePhase) {
535  ErrString += "\t - The phase became unstable (Cp < 0) T_unstable_last = "
536  + fp2str(Tunstable) + "\n";
537  }
538  if (doUV) {
539  throw CanteraError("setState_HPorUV (UV)", ErrString);
540  } else {
541  throw CanteraError("setState_HPorUV (HP)", ErrString);
542  }
543 }
544 //=================================================================================================================
545 
546 void ThermoPhase::setState_SP(doublereal Starget, doublereal p,
547  doublereal dTtol)
548 {
549  setState_SPorSV(Starget, p, dTtol, false);
550 }
551 //=================================================================================================================
552 
553 void ThermoPhase::setState_SV(doublereal Starget, doublereal v,
554  doublereal dTtol)
555 {
556  setState_SPorSV(Starget, v, dTtol, true);
557 }
558 //=================================================================================================================
559 
560 // Do the convergence work for fixed entropy situations
561 /*
562  * We assume here that S at constant P is a monotonically increasing
563  * function of T.
564  * We assume here that S at constant V is a monotonically increasing
565  * function of T.
566  *
567  * Note, the value of dTtol may become important for some applications
568  * where numerical jacobians are being calculated.
569  */
570 void ThermoPhase::setState_SPorSV(doublereal Starget, doublereal p,
571  doublereal dTtol, bool doSV)
572 {
573  doublereal v = 0.0;
574  doublereal dt;
575  if (doSV) {
576  v = p;
577  if (v < 1.0E-300) {
578  throw CanteraError("setState_SPorSV (SV)",
579  "Input specific volume is too small or negative. v = " + fp2str(v));
580  }
581  setDensity(1.0/v);
582  } else {
583  if (p < 1.0E-300) {
584  throw CanteraError("setState_SPorSV (SP)",
585  "Input pressure is too small or negative. p = " + fp2str(p));
586  }
587  setPressure(p);
588  }
589  double Tmax = maxTemp() + 0.1;
590  double Tmin = minTemp() - 0.1;
591 
592  // Make sure we are within the temperature bounds at the start
593  // of the iteration
594  double Tnew = temperature();
595  double Tinit = Tnew;
596  if (Tnew > Tmax) {
597  Tnew = Tmax - 1.0;
598  if (doSV) {
599  setTemperature(Tnew);
600  } else {
601  setState_TP(Tnew, p);
602  }
603  }
604  if (Tnew < Tmin) {
605  Tnew = Tmin + 1.0;
606  if (doSV) {
607  setTemperature(Tnew);
608  } else {
609  setState_TP(Tnew, p);
610  }
611  }
612 
613  double Snew = entropy_mass();
614  double Cpnew = 0.0;
615  if (doSV) {
616  Cpnew = cv_mass();
617  } else {
618  Cpnew = cp_mass();
619  }
620 
621  double Stop = Snew;
622  double Ttop = Tnew;
623  double Sbot = Snew;
624  double Tbot = Tnew;
625  double Told = Tnew;
626  double Sold = Snew;
627 
628  bool ignoreBounds = false;
629  // Unstable phases are those for which
630  // Cp < 0.0. These are possible for cases where
631  // we have passed the spinodal curve.
632  bool unstablePhase = false;
633  double Tunstable = -1.0;
634  bool unstablePhaseNew = false;
635 
636 
637  // Newton iteration
638  for (int n = 0; n < 500; n++) {
639  Told = Tnew;
640  Sold = Snew;
641  double cpd = Cpnew;
642  if (cpd < 0.0) {
643  unstablePhase = true;
644  Tunstable = Tnew;
645  }
646  dt = (Starget - Sold)*Told/cpd;
647 
648  // limit step size to 200 K
649  if (dt > 100.0) {
650  dt = 100.0;
651  } else if (dt < -100.0) {
652  dt = -100.0;
653  }
654  Tnew = Told + dt;
655  // Limit the step size so that we are convergent
656  if (dt > 0.0) {
657  if (!unstablePhase) {
658  if (Stop > Starget) {
659  if (Tnew > Ttop) {
660  dt = 0.75 * (Ttop - Told);
661  Tnew = Told + dt;
662  }
663  }
664  } else {
665  if (Sbot < Starget) {
666  if (Tnew < Tbot) {
667  dt = 0.75 * (Tbot - Told);
668  Tnew = Told + dt;
669  }
670  }
671  }
672  } else {
673  if (!unstablePhase) {
674  if (Sbot < Starget) {
675  if (Tnew < Tbot) {
676  dt = 0.75 * (Tbot - Told);
677  Tnew = Told + dt;
678  }
679  }
680  } else {
681  if (Stop > Starget) {
682  if (Tnew > Ttop) {
683  dt = 0.75 * (Ttop - Told);
684  Tnew = Told + dt;
685  }
686  }
687  }
688  }
689  // Check Max and Min values
690  if (Tnew > Tmax) {
691  if (!ignoreBounds) {
692  if (doSV) {
693  setTemperature(Tmax);
694  } else {
695  setState_TP(Tmax, p);
696  }
697  double Smax = entropy_mass();
698  if (Smax >= Starget) {
699  if (Stop < Starget) {
700  Ttop = Tmax;
701  Stop = Smax;
702  }
703  } else {
704  Tnew = Tmax + 1.0;
705  ignoreBounds = true;
706  }
707  }
708  }
709  if (Tnew < Tmin) {
710  if (!ignoreBounds) {
711  if (doSV) {
712  setTemperature(Tmin);
713  } else {
714  setState_TP(Tmin, p);
715  }
716  double Smin = enthalpy_mass();
717  if (Smin <= Starget) {
718  if (Sbot > Starget) {
719  Sbot = Tmin;
720  Sbot = Smin;
721  }
722  } else {
723  Tnew = Tmin - 1.0;
724  ignoreBounds = true;
725  }
726  }
727  }
728 
729  // Try to keep phase within its region of stability
730  // -> Could do a lot better if I calculate the
731  // spinodal value of H.
732  for (int its = 0; its < 10; its++) {
733  Tnew = Told + dt;
734  if (doSV) {
735  setTemperature(Tnew);
736  Cpnew = cv_mass();
737  } else {
738  setState_TP(Tnew, p);
739  Cpnew = cp_mass();
740  }
741  Snew = entropy_mass();
742  if (Cpnew < 0.0) {
743  unstablePhaseNew = true;
744  Tunstable = Tnew;
745  } else {
746  unstablePhaseNew = false;
747  break;
748  }
749  if (unstablePhase == false) {
750  if (unstablePhaseNew == true) {
751  dt *= 0.25;
752  }
753  }
754  }
755 
756  if (Snew == Starget) {
757  return;
758  } else if (Snew > Starget) {
759  if ((Stop < Starget) || (Snew < Stop)) {
760  Stop = Snew;
761  Ttop = Tnew;
762  }
763  } else if (Snew < Starget) {
764  if ((Sbot > Starget) || (Snew > Sbot)) {
765  Sbot = Snew;
766  Tbot = Tnew;
767  }
768  }
769  // Convergence in S
770  double Serr = Starget - Snew;
771  double acpd = std::max(fabs(cpd), 1.0E-5);
772  double denom = std::max(fabs(Starget), acpd * dTtol);
773  double SConvErr = fabs((Serr * Tnew)/denom);
774  if (SConvErr < 0.00001 *dTtol) {
775  return;
776  }
777  if (fabs(dt) < dTtol) {
778  return;
779  }
780  }
781  // We are here when there hasn't been convergence
782  /*
783  * Formulate a detailed error message, since questions seem to
784  * arise often about the lack of convergence.
785  */
786  string ErrString = "No convergence in 500 iterations\n";
787  if (doSV) {
788  ErrString += "\tTarget Entropy = " + fp2str(Starget) + "\n";
789  ErrString += "\tCurrent Specific Volume = " + fp2str(v) + "\n";
790  ErrString += "\tStarting Temperature = " + fp2str(Tinit) + "\n";
791  ErrString += "\tCurrent Temperature = " + fp2str(Tnew) + "\n";
792  ErrString += "\tCurrent Entropy = " + fp2str(Snew) + "\n";
793  ErrString += "\tCurrent Delta T = " + fp2str(dt) + "\n";
794  } else {
795  ErrString += "\tTarget Entropy = " + fp2str(Starget) + "\n";
796  ErrString += "\tCurrent Pressure = " + fp2str(p) + "\n";
797  ErrString += "\tStarting Temperature = " + fp2str(Tinit) + "\n";
798  ErrString += "\tCurrent Temperature = " + fp2str(Tnew) + "\n";
799  ErrString += "\tCurrent Entropy = " + fp2str(Snew) + "\n";
800  ErrString += "\tCurrent Delta T = " + fp2str(dt) + "\n";
801  }
802  if (unstablePhase) {
803  ErrString += "\t - The phase became unstable (Cp < 0) T_unstable_last = "
804  + fp2str(Tunstable) + "\n";
805  }
806  if (doSV) {
807  throw CanteraError("setState_SPorSV (SV)", ErrString);
808  } else {
809  throw CanteraError("setState_SPorSV (SP)", ErrString);
810  }
811 }
812 //=================================================================================================================
813 
814 doublereal ThermoPhase::err(std::string msg) const
815 {
816  throw CanteraError("ThermoPhase","Base class method "
817  +msg+" called. Equation of state type: "+int2str(eosType()));
818  return 0.0;
819 }
820 
821 /*
822  * Returns the units of the standard and general concentrations
823  * Note they have the same units, as their divisor is
824  * defined to be equal to the activity of the kth species
825  * in the solution, which is unitless.
826  *
827  * This routine is used in print out applications where the
828  * units are needed. Usually, MKS units are assumed throughout
829  * the program and in the XML input files.
830  *
831  * On return uA contains the powers of the units (MKS assumed)
832  * of the standard concentrations and generalized concentrations
833  * for the kth species.
834  *
835  * The base %ThermoPhase class assigns the default quantities
836  * of (kmol/m3).
837  * Inherited classes are responsible for overriding the default
838  * values if necessary.
839  *
840  * uA[0] = kmol units - default = 1
841  * uA[1] = m units - default = -nDim(), the number of spatial
842  * dimensions in the Phase class.
843  * uA[2] = kg units - default = 0;
844  * uA[3] = Pa(pressure) units - default = 0;
845  * uA[4] = Temperature units - default = 0;
846  * uA[5] = time units - default = 0
847  */
848 void ThermoPhase::getUnitsStandardConc(double* uA, int k, int sizeUA) const
849 {
850  for (int i = 0; i < sizeUA; i++) {
851  if (i == 0) {
852  uA[0] = 1.0;
853  }
854  if (i == 1) {
855  uA[1] = -int(nDim());
856  }
857  if (i == 2) {
858  uA[2] = 0.0;
859  }
860  if (i == 3) {
861  uA[3] = 0.0;
862  }
863  if (i == 4) {
864  uA[4] = 0.0;
865  }
866  if (i == 5) {
867  uA[5] = 0.0;
868  }
869  }
870 }
871 //=================================================================================================================
872 // Install a species thermodynamic property manager.
873 /*
874  * The species thermodynamic property manager
875  * computes properties of the pure species for use in
876  * constructing solution properties. It is meant for internal
877  * use, and some classes derived from ThermoPhase may not use
878  * any species thermodynamic property manager. This method is
879  * called by function importPhase() in importCTML.cpp.
880  *
881  * @param spthermo input pointer to the species thermodynamic property
882  * manager.
883  *
884  * @internal
885  */
887 {
888  if (m_spthermo) {
889  if (m_spthermo != spthermo) {
890  delete m_spthermo;
891  }
892  }
893  m_spthermo = spthermo;
894 }
895 //=================================================================================================================
896 // Return a changeable reference to the calculation manager
897 // for species reference-state thermodynamic properties
898 /*
899  *
900  * @param k Speices id. The default is -1, meaning return the default
901  *
902  * @internal
903  */
905 {
906  if (!m_spthermo) {
907  throw CanteraError("ThermoPhase::speciesThermo()",
908  "species reference state thermo manager was not set");
909  }
910  return *m_spthermo;
911 }
912 //=================================================================================================================
913 /*
914  * initThermoFile():
915  *
916  * Initialization of a phase using an xml file.
917  *
918  * This routine is a precursor to initThermoXML(XML_Node*)
919  * routine, which does most of the work.
920  *
921  * @param infile XML file containing the description of the
922  * phase
923  *
924  * @param id Optional parameter identifying the name of the
925  * phase. If none is given, the first XML
926  * phase element will be used.
927  */
928 void ThermoPhase::initThermoFile(std::string inputFile, std::string id)
929 {
930 
931  if (inputFile.size() == 0) {
932  throw CanteraError("ThermoPhase::initThermoFile",
933  "input file is null");
934  }
935  string path = findInputFile(inputFile);
936  ifstream fin(path.c_str());
937  if (!fin) {
938  throw CanteraError("initThermoFile","could not open "
939  +path+" for reading.");
940  }
941  /*
942  * The phase object automatically constructs an XML object.
943  * Use this object to store information.
944  */
945  XML_Node& phaseNode_XML = xml();
946  XML_Node* fxml = new XML_Node();
947  fxml->build(fin);
948  XML_Node* fxml_phase = findXMLPhase(fxml, id);
949  if (!fxml_phase) {
950  throw CanteraError("ThermoPhase::initThermo",
951  "ERROR: Can not find phase named " +
952  id + " in file named " + inputFile);
953  }
954  fxml_phase->copy(&phaseNode_XML);
955  initThermoXML(*fxml_phase, id);
956  delete fxml;
957 }
958 //=================================================================================================================
959 
960 /*
961  * Import and initialize a ThermoPhase object
962  *
963  * This function is called from importPhase()
964  * after the elements and the
965  * species are initialized with default ideal solution
966  * level data.
967  *
968  * @param phaseNode This object must be the phase node of a
969  * complete XML tree
970  * description of the phase, including all of the
971  * species data. In other words while "phase" must
972  * point to an XML phase object, it must have
973  * sibling nodes "speciesData" that describe
974  * the species in the phase.
975  * @param id ID of the phase. If nonnull, a check is done
976  * to see if phaseNode is pointing to the phase
977  * with the correct id.
978  */
979 void ThermoPhase::initThermoXML(XML_Node& phaseNode, std::string id)
980 {
981 
982  /*
983  * and sets the state
984  */
985  if (phaseNode.hasChild("state")) {
986  XML_Node& stateNode = phaseNode.child("state");
987  setStateFromXML(stateNode);
988  }
990 }
991 
992 void ThermoPhase::setReferenceComposition(const doublereal* const x)
993 {
994  xMol_Ref.resize(m_kk);
995  if (x) {
996  for (size_t k = 0; k < m_kk; k++) {
997  xMol_Ref[k] = x[k];
998  }
999  } else {
1001  }
1002  double sum = -1.0;
1003  for (size_t k = 0; k < m_kk; k++) {
1004  sum += xMol_Ref[k];
1005  }
1006  if (fabs(sum) > 1.0E-11) {
1007  throw CanteraError("ThermoPhase::setReferenceComposition",
1008  "input mole fractions don't sum to 1.0");
1009  }
1010 
1011 }
1012 
1013 void ThermoPhase::getReferenceComposition(doublereal* const x) const
1014 {
1015  for (size_t k = 0; k < m_kk; k++) {
1016  x[k] = xMol_Ref[k];
1017  }
1018 }
1019 
1020 /*
1021  * Initialize.
1022  *
1023  * This method is provided to allow
1024  * subclasses to perform any initialization required after all
1025  * species have been added. For example, it might be used to
1026  * resize internal work arrays that must have an entry for
1027  * each species. The base class implementation does nothing,
1028  * and subclasses that do not require initialization do not
1029  * need to overload this method. When importing a CTML phase
1030  * description, this method is called just prior to returning
1031  * from function importPhase.
1032  *
1033  * @see importCTML.cpp
1034  */
1036 {
1037  // Check to see that there is at least one species defined in the phase
1038  if (m_kk == 0) {
1039  throw CanteraError("ThermoPhase::initThermo()",
1040  "Number of species is equal to zero");
1041  }
1042  xMol_Ref.resize(m_kk, 0.0);
1043 }
1044 //====================================================================================================================
1046 {
1047 
1048 }
1049 //====================================================================================================================
1050 void ThermoPhase::saveSpeciesData(const size_t k, const XML_Node* const data)
1051 {
1052  if (m_speciesData.size() < (k + 1)) {
1053  m_speciesData.resize(k+1, 0);
1054  }
1055  m_speciesData[k] = new XML_Node(*data);
1056 }
1057 //====================================================================================================================
1058 // Return a pointer to the XML tree containing the species
1059 // data for this phase.
1060 const std::vector<const XML_Node*> & ThermoPhase::speciesData() const
1061 {
1062  if (m_speciesData.size() != m_kk) {
1063  throw CanteraError("ThermoPhase::speciesData",
1064  "m_speciesData is the wrong size");
1065  }
1066  return m_speciesData;
1067 }
1068 //====================================================================================================================
1069 /*
1070  * Set the thermodynamic state.
1071  */
1073 {
1074  string comp = getChildValue(state,"moleFractions");
1075  if (comp != "") {
1076  setMoleFractionsByName(comp);
1077  } else {
1078  comp = getChildValue(state,"massFractions");
1079  if (comp != "") {
1080  setMassFractionsByName(comp);
1081  }
1082  }
1083  if (state.hasChild("temperature")) {
1084  double t = getFloat(state, "temperature", "temperature");
1085  setTemperature(t);
1086  }
1087  if (state.hasChild("pressure")) {
1088  double p = getFloat(state, "pressure", "pressure");
1089  setPressure(p);
1090  }
1091  if (state.hasChild("density")) {
1092  double rho = getFloat(state, "density", "density");
1093  setDensity(rho);
1094  }
1095 }
1096 //====================================================================================================================
1097 /*
1098  * Called by function 'equilibrate' in ChemEquil.h to transfer
1099  * the element potentials to this object after every successful
1100  * equilibration routine.
1101  * The element potentials are stored in their dimensionless
1102  * forms, calculated by dividing by RT.
1103  * @param lambda vector containing the element potentials.
1104  * Length = nElements. Units are Joules/kmol.
1105  */
1107 {
1108  doublereal rrt = 1.0/(GasConstant* temperature());
1109  size_t mm = nElements();
1110  if (lambda.size() < mm) {
1111  throw CanteraError("setElementPotentials", "lambda too small");
1112  }
1113  if (!m_hasElementPotentials) {
1114  m_lambdaRRT.resize(mm);
1115  }
1116  for (size_t m = 0; m < mm; m++) {
1117  m_lambdaRRT[m] = lambda[m] * rrt;
1118  }
1119  m_hasElementPotentials = true;
1120 }
1121 
1122 /*
1123  * Returns the stored element potentials.
1124  * The element potentials are retrieved from their stored
1125  * dimensionless forms by multiplying by RT.
1126  * @param lambda Vector containing the element potentials.
1127  * Length = nElements. Units are Joules/kmol.
1128  */
1129 bool ThermoPhase::getElementPotentials(doublereal* lambda) const
1130 {
1131  doublereal rt = GasConstant* temperature();
1132  if (m_hasElementPotentials) {
1133  for (size_t m = 0; m < nElements(); m++) {
1134  lambda[m] = m_lambdaRRT[m] * rt;
1135  }
1136  }
1137  return (m_hasElementPotentials);
1138 }
1139 //====================================================================================================================
1140 // Get the array of derivatives of the log activity coefficients with respect to the species mole numbers
1141 /*
1142  * Implementations should take the derivative of the logarithm of the activity coefficient with respect to a
1143  * species mole number (with all other species mole numbers held constant)
1144  *
1145  * units = 1 / kmol
1146  *
1147  * dlnActCoeffdN[ ld * k + m] will contain the derivative of log act_coeff for the <I>m</I><SUP>th</SUP>
1148  * species with respect to the number of moles of the <I>k</I><SUP>th</SUP> species.
1149  *
1150  * \f[
1151  * \frac{d \ln(\gamma_m) }{d n_k }\Bigg|_{n_i}
1152  * \f]
1153  *
1154  * @param ld Number of rows in the matrix
1155  * @param dlnActCoeffdN Output vector of derivatives of the
1156  * log Activity Coefficients. length = m_kk * m_kk
1157  */
1158 void ThermoPhase::getdlnActCoeffdlnN(const size_t ld, doublereal* const dlnActCoeffdlnN)
1159 {
1160  for (size_t m = 0; m < m_kk; m++) {
1161  for (size_t k = 0; k < m_kk; k++) {
1162  dlnActCoeffdlnN[ld * k + m] = 0.0;
1163  }
1164  }
1165  return;
1166 }
1167 //====================================================================================================================
1168 void ThermoPhase::getdlnActCoeffdlnN_numderiv(const size_t ld, doublereal* const dlnActCoeffdlnN)
1169 {
1170  double deltaMoles_j = 0.0;
1171  double pres = pressure();
1172 
1173  /*
1174  * Evaluate the current base activity coefficients if necessary
1175  */
1176  std::vector<double> ActCoeff_Base(m_kk);
1177  getActivityCoefficients(DATA_PTR(ActCoeff_Base));
1178  std::vector<double> Xmol_Base(m_kk);
1179  getMoleFractions(DATA_PTR(Xmol_Base));
1180 
1181  // Make copies of ActCoeff and Xmol_ for use in taking differences
1182  std::vector<double> ActCoeff(m_kk);
1183  std::vector<double> Xmol(m_kk);
1184  double v_totalMoles = 1.0;
1185  double TMoles_base = v_totalMoles;
1186 
1187  /*
1188  * Loop over the columns species to be deltad
1189  */
1190  for (size_t j = 0; j < m_kk; j++) {
1191  /*
1192  * Calculate a value for the delta moles of species j
1193  * -> NOte Xmol_[] and Tmoles are always positive or zero
1194  * quantities.
1195  * -> experience has shown that you always need to make the deltas greater than needed to
1196  * change the other mole fractions in order to capture some effects.
1197  */
1198  double moles_j_base = v_totalMoles * Xmol_Base[j];
1199  deltaMoles_j = 1.0E-7 * moles_j_base + v_totalMoles * 1.0E-13 + 1.0E-150;
1200  /*
1201  * Now, update the total moles in the phase and all of the
1202  * mole fractions based on this.
1203  */
1204  v_totalMoles = TMoles_base + deltaMoles_j;
1205  for (size_t k = 0; k < m_kk; k++) {
1206  Xmol[k] = Xmol_Base[k] * TMoles_base / v_totalMoles;
1207  }
1208  Xmol[j] = (moles_j_base + deltaMoles_j) / v_totalMoles;
1209 
1210  /*
1211  * Go get new values for the activity coefficients.
1212  * -> Note this calls setState_PX();
1213  */
1214  setState_PX(pres, DATA_PTR(Xmol));
1215  getActivityCoefficients(DATA_PTR(ActCoeff));
1216 
1217  /*
1218  * Calculate the column of the matrix
1219  */
1220  double* const lnActCoeffCol = dlnActCoeffdlnN + ld * j;
1221  for (size_t k = 0; k < m_kk; k++) {
1222  lnActCoeffCol[k] = (2*moles_j_base + deltaMoles_j) *(ActCoeff[k] - ActCoeff_Base[k]) /
1223  ((ActCoeff[k] + ActCoeff_Base[k]) * deltaMoles_j);
1224  }
1225  /*
1226  * Revert to the base case Xmol_, v_totalMoles
1227  */
1228  v_totalMoles = TMoles_base;
1229  mdp::mdp_copy_dbl_1(DATA_PTR(Xmol), DATA_PTR(Xmol_Base), (int) m_kk);
1230  }
1231  /*
1232  * Go get base values for the activity coefficients.
1233  * -> Note this calls setState_TPX() again;
1234  * -> Just wanted to make sure that cantera is in sync
1235  * with VolPhase after this call.
1236  */
1237  setState_PX(pres, DATA_PTR(Xmol_Base));
1238 }
1239 //====================================================================================================================
1240 /*
1241  * Format a summary of the mixture state for output.
1242  */
1243 std::string ThermoPhase::report(bool show_thermo) const
1244 {
1245  char p[800];
1246  string s = "";
1247  try {
1248  if (name() != "") {
1249  sprintf(p, " \n %s:\n", name().c_str());
1250  s += p;
1251  }
1252  sprintf(p, " \n temperature %12.6g K\n", temperature());
1253  s += p;
1254  sprintf(p, " pressure %12.6g Pa\n", pressure());
1255  s += p;
1256  sprintf(p, " density %12.6g kg/m^3\n", density());
1257  s += p;
1258  sprintf(p, " mean mol. weight %12.6g amu\n", meanMolecularWeight());
1259  s += p;
1260 
1261  doublereal phi = electricPotential();
1262  if (phi != 0.0) {
1263  sprintf(p, " potential %12.6g V\n", phi);
1264  s += p;
1265  }
1266  if (show_thermo) {
1267  sprintf(p, " \n");
1268  s += p;
1269  sprintf(p, " 1 kg 1 kmol\n");
1270  s += p;
1271  sprintf(p, " ----------- ------------\n");
1272  s += p;
1273  sprintf(p, " enthalpy %12.6g %12.4g J\n",
1275  s += p;
1276  sprintf(p, " internal energy %12.6g %12.4g J\n",
1278  s += p;
1279  sprintf(p, " entropy %12.6g %12.4g J/K\n",
1280  entropy_mass(), entropy_mole());
1281  s += p;
1282  sprintf(p, " Gibbs function %12.6g %12.4g J\n",
1283  gibbs_mass(), gibbs_mole());
1284  s += p;
1285  sprintf(p, " heat capacity c_p %12.6g %12.4g J/K\n",
1286  cp_mass(), cp_mole());
1287  s += p;
1288  try {
1289  sprintf(p, " heat capacity c_v %12.6g %12.4g J/K\n",
1290  cv_mass(), cv_mole());
1291  s += p;
1292  } catch (CanteraError& err) {
1293  err.save();
1294  sprintf(p, " heat capacity c_v <not implemented> \n");
1295  s += p;
1296  }
1297  }
1298 
1299  size_t kk = nSpecies();
1300  vector_fp x(kk);
1301  vector_fp y(kk);
1302  vector_fp mu(kk);
1303  getMoleFractions(&x[0]);
1304  getMassFractions(&y[0]);
1305  getChemPotentials(&mu[0]);
1306  doublereal rt = GasConstant * temperature();
1307  //if (th.nSpecies() > 1) {
1308 
1309  if (show_thermo) {
1310  sprintf(p, " \n X "
1311  " Y Chem. Pot. / RT \n");
1312  s += p;
1313  sprintf(p, " ------------- "
1314  "------------ ------------\n");
1315  s += p;
1316  for (size_t k = 0; k < kk; k++) {
1317  if (x[k] > SmallNumber) {
1318  sprintf(p, "%18s %12.6g %12.6g %12.6g\n",
1319  speciesName(k).c_str(), x[k], y[k], mu[k]/rt);
1320  } else {
1321  sprintf(p, "%18s %12.6g %12.6g \n",
1322  speciesName(k).c_str(), x[k], y[k]);
1323  }
1324  s += p;
1325  }
1326  } else {
1327  sprintf(p, " \n X"
1328  "Y\n");
1329  s += p;
1330  sprintf(p, " -------------"
1331  " ------------\n");
1332  s += p;
1333  for (size_t k = 0; k < kk; k++) {
1334  sprintf(p, "%18s %12.6g %12.6g\n",
1335  speciesName(k).c_str(), x[k], y[k]);
1336  s += p;
1337  }
1338  }
1339  }
1340  //}
1341  catch (CanteraError& err) {
1342  err.save();
1343  }
1344  return s;
1345 }
1346 //====================================================================================================================
1347 /*
1348  * Format a summary of the mixture state for output.
1349  */
1350 void ThermoPhase::reportCSV(std::ofstream& csvFile) const
1351 {
1352 
1353  csvFile.precision(3);
1354  int tabS = 15;
1355  int tabM = 30;
1356  int tabL = 40;
1357  try {
1358  if (name() != "") {
1359  csvFile << "\n"+name()+"\n\n";
1360  }
1361  csvFile << setw(tabL) << "temperature (K) =" << setw(tabS) << temperature() << endl;
1362  csvFile << setw(tabL) << "pressure (Pa) =" << setw(tabS) << pressure() << endl;
1363  csvFile << setw(tabL) << "density (kg/m^3) =" << setw(tabS) << density() << endl;
1364  csvFile << setw(tabL) << "mean mol. weight (amu) =" << setw(tabS) << meanMolecularWeight() << endl;
1365  csvFile << setw(tabL) << "potential (V) =" << setw(tabS) << electricPotential() << endl;
1366  csvFile << endl;
1367 
1368  csvFile << setw(tabL) << "enthalpy (J/kg) = " << setw(tabS) << enthalpy_mass() << setw(tabL)
1369  << "enthalpy (J/kmol) = " << setw(tabS) << enthalpy_mole() << endl;
1370  csvFile << setw(tabL) << "internal E (J/kg) = " << setw(tabS) << intEnergy_mass() << setw(tabL)
1371  << "internal E (J/kmol) = " << setw(tabS) << intEnergy_mole() << endl;
1372  csvFile << setw(tabL) << "entropy (J/kg) = " << setw(tabS) << entropy_mass() << setw(tabL)
1373  << "entropy (J/kmol) = " << setw(tabS) << entropy_mole() << endl;
1374  csvFile << setw(tabL) << "Gibbs (J/kg) = " << setw(tabS) << gibbs_mass() << setw(tabL)
1375  << "Gibbs (J/kmol) = " << setw(tabS) << gibbs_mole() << endl;
1376  csvFile << setw(tabL) << "heat capacity c_p (J/K/kg) = " << setw(tabS) << cp_mass()
1377  << setw(tabL) << "heat capacity c_p (J/K/kmol) = " << setw(tabS) << cp_mole() << endl;
1378  csvFile << setw(tabL) << "heat capacity c_v (J/K/kg) = " << setw(tabS) << cv_mass()
1379  << setw(tabL) << "heat capacity c_v (J/K/kmol) = " << setw(tabS) << cv_mole() << endl;
1380 
1381  csvFile.precision(8);
1382 
1383  size_t kk = nSpecies();
1384  doublereal* x = new doublereal[kk];
1385  doublereal* y = new doublereal[kk];
1386  doublereal* mu = new doublereal[kk];
1387  doublereal* a = new doublereal[kk];
1388  doublereal* ac = new doublereal[kk];
1389  doublereal* hbar = new doublereal[kk];
1390  doublereal* sbar = new doublereal[kk];
1391  doublereal* ubar = new doublereal[kk];
1392  doublereal* cpbar= new doublereal[kk];
1393  doublereal* vbar = new doublereal[kk];
1394  std::vector<std::string> pNames;
1395  std::vector<doublereal*> data;
1396 
1397  getMoleFractions(x);
1398  pNames.push_back("X");
1399  data.push_back(x);
1400  try {
1401  getMassFractions(y);
1402  pNames.push_back("Y");
1403  data.push_back(y);
1404  } catch (CanteraError& err) {
1405  err.save();
1406  }
1407  try {
1408  getChemPotentials(mu);
1409  pNames.push_back("Chem. Pot (J/kmol)");
1410  data.push_back(mu);
1411  } catch (CanteraError& err) {
1412  err.save();
1413  }
1414  try {
1415  getActivities(a);
1416  pNames.push_back("Activity");
1417  data.push_back(a);
1418  } catch (CanteraError& err) {
1419  err.save();
1420  }
1421  try {
1423  pNames.push_back("Act. Coeff.");
1424  data.push_back(ac);
1425  } catch (CanteraError& err) {
1426  err.save();
1427  }
1428  try {
1430  pNames.push_back("Part. Mol Enthalpy (J/kmol)");
1431  data.push_back(hbar);
1432  } catch (CanteraError& err) {
1433  err.save();
1434  }
1435  try {
1437  pNames.push_back("Part. Mol. Entropy (J/K/kmol)");
1438  data.push_back(sbar);
1439  } catch (CanteraError& err) {
1440  err.save();
1441  }
1442  try {
1444  pNames.push_back("Part. Mol. Energy (J/kmol)");
1445  data.push_back(ubar);
1446  } catch (CanteraError& err) {
1447  err.save();
1448  }
1449  try {
1450  getPartialMolarCp(cpbar);
1451  pNames.push_back("Part. Mol. Cp (J/K/kmol");
1452  data.push_back(cpbar);
1453  } catch (CanteraError& err) {
1454  err.save();
1455  }
1456  try {
1457  getPartialMolarVolumes(vbar);
1458  pNames.push_back("Part. Mol. Cv (J/K/kmol)");
1459  data.push_back(vbar);
1460  } catch (CanteraError& err) {
1461  err.save();
1462  }
1463 
1464  csvFile << endl << setw(tabS) << "Species,";
1465  for (size_t i = 0; i < pNames.size(); i++) {
1466  csvFile << setw(tabM) << pNames[i] << ",";
1467  }
1468  csvFile << endl;
1469  /*
1470  csvFile.fill('-');
1471  csvFile << setw(tabS+(tabM+1)*pNames.size()) << "-\n";
1472  csvFile.fill(' ');
1473  */
1474  for (size_t k = 0; k < kk; k++) {
1475  csvFile << setw(tabS) << speciesName(k) + ",";
1476  if (x[k] > SmallNumber) {
1477  for (size_t i = 0; i < pNames.size(); i++) {
1478  csvFile << setw(tabM) << data[i][k] << ",";
1479  }
1480  csvFile << endl;
1481  } else {
1482  for (size_t i = 0; i < pNames.size(); i++) {
1483  csvFile << setw(tabM) << 0 << ",";
1484  }
1485  csvFile << endl;
1486  }
1487  }
1488  delete [] x;
1489  delete [] y;
1490  delete [] mu;
1491  delete [] a;
1492  delete [] ac;
1493  delete [] hbar;
1494  delete [] sbar;
1495  delete [] ubar;
1496  delete [] cpbar;
1497  delete [] vbar;
1498 
1499  } catch (CanteraError& err) {
1500  err.save();
1501  }
1502 }
1503 
1504 std::string report(const ThermoPhase& th, const bool show_thermo)
1505 {
1506  return th.report(show_thermo);
1507 }
1508 
1509 }