Cantera  2.0
MolarityIonicVPSSTP.cpp
Go to the documentation of this file.
1 /**
2  * @file MolarityIonicVPSSTP.cpp
3  * Definitions for intermediate ThermoPhase object for phases which
4  * employ excess gibbs free energy formulations
5  * (see \ref thermoprops
6  * and class \link Cantera::MolarityIonicVPSSTP MolarityIonicVPSSTP\endlink).
7  *
8  * Header file for a derived class of ThermoPhase that handles
9  * variable pressure standard state methods for calculating
10  * thermodynamic properties that are further based upon expressions
11  * for the excess gibbs free energy expressed as a function of
12  * the mole fractions.
13  */
14 /*
15  * Copyright (2009) Sandia Corporation. Under the terms of
16  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
17  * U.S. Government retains certain rights in this software.
18  */
22 
23 #include <cmath>
24 #include <fstream>
25 
26 using namespace std;
27 
28 namespace Cantera
29 {
30 static const double xxSmall = 1.0E-150;
31 //====================================================================================================================
32 /*
33  * Default constructor.
34  *
35  */
36 MolarityIonicVPSSTP::MolarityIonicVPSSTP() :
38  PBType_(PBTYPE_PASSTHROUGH),
39  numPBSpecies_(m_kk),
40  indexSpecialSpecies_(npos),
41  numCationSpecies_(0),
42  numAnionSpecies_(0),
43  numPassThroughSpecies_(0),
44  neutralPBindexStart(0)
45 {
46 }
47 //====================================================================================================================
48 /*
49  * Working constructors
50  *
51  * The two constructors below are the normal way
52  * the phase initializes itself. They are shells that call
53  * the routine initThermo(), with a reference to the
54  * XML database to get the info for the phase.
55  */
56 MolarityIonicVPSSTP::MolarityIonicVPSSTP(std::string inputFile, std::string id) :
58  PBType_(PBTYPE_PASSTHROUGH),
59  numPBSpecies_(m_kk),
60  indexSpecialSpecies_(npos),
61  numCationSpecies_(0),
62  numAnionSpecies_(0),
63  numPassThroughSpecies_(0),
64  neutralPBindexStart(0)
65 {
66  constructPhaseFile(inputFile, id);
67 }
68 //====================================================================================================================
71  PBType_(PBTYPE_PASSTHROUGH),
72  numPBSpecies_(m_kk),
73  indexSpecialSpecies_(npos),
74  numCationSpecies_(0),
75  numAnionSpecies_(0),
76  numPassThroughSpecies_(0),
77  neutralPBindexStart(0)
78 {
79  constructPhaseXML(phaseRoot, id);
80 }
81 //====================================================================================================================
82 /*
83  * Copy Constructor:
84  *
85  * Note this stuff will not work until the underlying phase
86  * has a working copy constructor
87  */
90  PBType_(PBTYPE_PASSTHROUGH),
91  numPBSpecies_(m_kk),
92  indexSpecialSpecies_(npos),
93  numCationSpecies_(0),
94  numAnionSpecies_(0),
95  numPassThroughSpecies_(0),
96  neutralPBindexStart(0)
97 {
98  *this = operator=(b);
99 }
100 //====================================================================================================================
101 /*
102  * operator=()
103  *
104  * Note this stuff will not work until the underlying phase
105  * has a working assignment operator
106  */
109 {
110  if (&b != this) {
112  }
113 
114  PBType_ = b.PBType_;
117  PBMoleFractions_ = b.PBMoleFractions_;
120  anionList_ = b.anionList_;
121  numAnionSpecies_ = b.numAnionSpecies_;
122  passThroughList_ = b.passThroughList_;
123  numPassThroughSpecies_ = b.numPassThroughSpecies_;
124  neutralPBindexStart = b.neutralPBindexStart;
125  moleFractionsTmp_ = b.moleFractionsTmp_;
126 
127  return *this;
128 }
129 //====================================================================================================================
130 /**
131  *
132  * ~MolarityIonicVPSSTP(): (virtual)
133  *
134  * Destructor: does nothing:
135  *
136  */
138 {
139 }
140 
141 /*
142  * This routine duplicates the current object and returns
143  * a pointer to ThermoPhase.
144  */
147 {
148  MolarityIonicVPSSTP* mtp = new MolarityIonicVPSSTP(*this);
149  return (ThermoPhase*) mtp;
150 }
151 
152 /*
153  * -------------- Utilities -------------------------------
154  */
155 //====================================================================================================================
156 
157 // Equation of state type flag.
158 /*
159  * The ThermoPhase base class returns
160  * zero. Subclasses should define this to return a unique
161  * non-zero value. Known constants defined for this purpose are
162  * listed in mix_defs.h. The MolarityIonicVPSSTP class also returns
163  * zero, as it is a non-complete class.
164  */
166 {
167  return 0;
168 }
169 
170 //====================================================================================================================
171 /*
172  * Import, construct, and initialize a phase
173  * specification from an XML tree into the current object.
174  *
175  * This routine is a precursor to constructPhaseXML(XML_Node*)
176  * routine, which does most of the work.
177  *
178  * @param infile XML file containing the description of the
179  * phase
180  *
181  * @param id Optional parameter identifying the name of the
182  * phase. If none is given, the first XML
183  * phase element will be used.
184  */
185 void MolarityIonicVPSSTP::constructPhaseFile(std::string inputFile, std::string id)
186 {
187 
188  if ((int) inputFile.size() == 0) {
189  throw CanteraError("MolarityIonicVPSSTP:constructPhaseFile",
190  "input file is null");
191  }
192  string path = findInputFile(inputFile);
193  std::ifstream fin(path.c_str());
194  if (!fin) {
195  throw CanteraError("MolarityIonicVPSSTP:constructPhaseFile","could not open "
196  +path+" for reading.");
197  }
198  /*
199  * The phase object automatically constructs an XML object.
200  * Use this object to store information.
201  */
202  XML_Node& phaseNode_XML = xml();
203  XML_Node* fxml = new XML_Node();
204  fxml->build(fin);
205  XML_Node* fxml_phase = findXMLPhase(fxml, id);
206  if (!fxml_phase) {
207  throw CanteraError("MolarityIonicVPSSTP:constructPhaseFile",
208  "ERROR: Can not find phase named " +
209  id + " in file named " + inputFile);
210  }
211  fxml_phase->copy(&phaseNode_XML);
212  constructPhaseXML(*fxml_phase, id);
213  delete fxml;
214 }
215 //====================================================================================================================
216 /*
217  * Import, construct, and initialize a HMWSoln phase
218  * specification from an XML tree into the current object.
219  *
220  * Most of the work is carried out by the cantera base
221  * routine, importPhase(). That routine imports all of the
222  * species and element data, including the standard states
223  * of the species.
224  *
225  * Then, In this routine, we read the information
226  * particular to the specification of the activity
227  * coefficient model for the Pitzer parameterization.
228  *
229  * We also read information about the molar volumes of the
230  * standard states if present in the XML file.
231  *
232  * @param phaseNode This object must be the phase node of a
233  * complete XML tree
234  * description of the phase, including all of the
235  * species data. In other words while "phase" must
236  * point to an XML phase object, it must have
237  * sibling nodes "speciesData" that describe
238  * the species in the phase.
239  * @param id ID of the phase. If nonnull, a check is done
240  * to see if phaseNode is pointing to the phase
241  * with the correct id.
242  */
243 void MolarityIonicVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string id)
244 {
245  string stemp;
246  if ((int) id.size() > 0) {
247  string idp = phaseNode.id();
248  if (idp != id) {
249  throw CanteraError("MolarityIonicVPSSTP::constructPhaseXML",
250  "phasenode and Id are incompatible");
251  }
252  }
253 
254  /*
255  * Find the Thermo XML node
256  */
257  if (!phaseNode.hasChild("thermo")) {
258  throw CanteraError("MolarityIonicVPSSTP::constructPhaseXML",
259  "no thermo XML node");
260  }
261  XML_Node& thermoNode = phaseNode.child("thermo");
262 
263  /*
264  * Make sure that the thermo model is MolarityIonic
265  */
266  stemp = thermoNode.attrib("model");
267  string formString = lowercase(stemp);
268  if (formString != "molarityionicvpss" && formString != "molarityionicvpsstp") {
269  throw CanteraError("MolarityIonicVPSSTP::constructPhaseXML",
270  "model name isn't MolarityIonicVPSSTP: " + formString);
271  }
272 
273  /*
274  * Call the Cantera importPhase() function. This will import
275  * all of the species into the phase. This will also handle
276  * all of the solvent and solute standard states
277  */
278  bool m_ok = importPhase(phaseNode, this);
279  if (!m_ok) {
280  throw CanteraError("MolarityIonicVPSSTP::constructPhaseXML","importPhase failed ");
281  }
282 
283 }
284 //====================================================================================================================
285 /*
286  * ------------ Molar Thermodynamic Properties ----------------------
287  */
288 //====================================================================================================================
289 /*
290  * - Activities, Standard States, Activity Concentrations -----------
291  */
292 //====================================================================================================================
293 
295 {
296  /*
297  * Update the activity coefficients
298  */
300 
301  /*
302  * take the exp of the internally stored coefficients.
303  */
304  for (size_t k = 0; k < m_kk; k++) {
305  lnac[k] = lnActCoeff_Scaled_[k];
306  }
307 }
308 //====================================================================================================================
309 void MolarityIonicVPSSTP::getChemPotentials(doublereal* mu) const
310 {
311  doublereal xx;
312  /*
313  * First get the standard chemical potentials in
314  * molar form.
315  * -> this requires updates of standard state as a function
316  * of T and P
317  */
319  /*
320  * Update the activity coefficients
321  */
323  /*
324  *
325  */
326  doublereal RT = GasConstant * temperature();
327  for (size_t k = 0; k < m_kk; k++) {
328  xx = std::max(moleFractions_[k], xxSmall);
329  mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]);
330  }
331 }
332 //====================================================================================================================
333 
335 {
336  getChemPotentials(mu);
337  double ve = Faraday * electricPotential();
338  for (size_t k = 0; k < m_kk; k++) {
339  mu[k] += ve*charge(k);
340  }
341 }
342 
343 //====================================================================================================================
344 // Returns an array of partial molar enthalpies for the species
345 // in the mixture.
346 /*
347  * Units (J/kmol)
348  *
349  * For this phase, the partial molar enthalpies are equal to the
350  * standard state enthalpies modified by the derivative of the
351  * molality-based activity coefficient wrt temperature
352  *
353  * \f[
354  * \bar h_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
355  * \f]
356  *
357  */
359 {
360  /*
361  * Get the nondimensional standard state enthalpies
362  */
363  getEnthalpy_RT(hbar);
364  /*
365  * dimensionalize it.
366  */
367  double T = temperature();
368  double RT = GasConstant * T;
369  for (size_t k = 0; k < m_kk; k++) {
370  hbar[k] *= RT;
371  }
372  /*
373  * Update the activity coefficients, This also update the
374  * internally stored molalities.
375  */
378  double RTT = RT * T;
379  for (size_t k = 0; k < m_kk; k++) {
380  hbar[k] -= RTT * dlnActCoeffdT_Scaled_[k];
381  }
382 }
383 //====================================================================================================================
384 // Returns an array of partial molar heat capacities for the species
385 // in the mixture.
386 /*
387  * Units (J/kmol)
388  *
389  * For this phase, the partial molar enthalpies are equal to the
390  * standard state enthalpies modified by the derivative of the
391  * activity coefficient wrt temperature
392  *
393  * \f[
394  * ??????????? \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
395  * \f]
396  *
397  */
398 void MolarityIonicVPSSTP::getPartialMolarCp(doublereal* cpbar) const
399 {
400  /*
401  * Get the nondimensional standard state entropies
402  */
403  getCp_R(cpbar);
404  double T = temperature();
405  /*
406  * Update the activity coefficients, This also update the
407  * internally stored molalities.
408  */
411 
412  for (size_t k = 0; k < m_kk; k++) {
413  cpbar[k] -= 2 * T * dlnActCoeffdT_Scaled_[k] + T * T * d2lnActCoeffdT2_Scaled_[k];
414  }
415  /*
416  * dimensionalize it.
417  */
418  for (size_t k = 0; k < m_kk; k++) {
419  cpbar[k] *= GasConstant;
420  }
421 }
422 //====================================================================================================================
423 // Returns an array of partial molar entropies for the species
424 // in the mixture.
425 /*
426  * Units (J/kmol)
427  *
428  * For this phase, the partial molar enthalpies are equal to the
429  * standard state enthalpies modified by the derivative of the
430  * activity coefficient wrt temperature
431  *
432  * \f[
433  * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
434  * \f]
435  *
436  */
438 {
439  double xx;
440  /*
441  * Get the nondimensional standard state entropies
442  */
443  getEntropy_R(sbar);
444  double T = temperature();
445  /*
446  * Update the activity coefficients, This also update the
447  * internally stored molalities.
448  */
451 
452  for (size_t k = 0; k < m_kk; k++) {
453  xx = std::max(moleFractions_[k], xxSmall);
454  sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
455  }
456  /*
457  * dimensionalize it.
458  */
459  for (size_t k = 0; k < m_kk; k++) {
460  sbar[k] *= GasConstant;
461  }
462 }
463 // Return an array of partial molar volumes for the
464 // species in the mixture. Units: m^3/kmol.
465 /*
466  * Frequently, for this class of thermodynamics representations,
467  * the excess Volume due to mixing is zero. Here, we set it as
468  * a default. It may be overridden in derived classes.
469  *
470  * @param vbar Output vector of species partial molar volumes.
471  * Length = m_kk. units are m^3/kmol.
472  */
474 {
475  /*
476  * Get the standard state values in m^3 kmol-1
477  */
478  getStandardVolumes(vbar);
479  for (size_t iK = 0; iK < m_kk; iK++) {
480  vbar[iK] += 0.0;
481  }
482 }
483 //====================================================================================================================
485 {
486  size_t k;
487  size_t kCat;
488  size_t kMax;
489  doublereal sumCat;
490  doublereal sumAnion;
491  doublereal chP, chM;
492  doublereal sum = 0.0;
493  doublereal sumMax;
494  switch (PBType_) {
495  case PBTYPE_PASSTHROUGH:
496  for (k = 0; k < m_kk; k++) {
497  PBMoleFractions_[k] = moleFractions_[k];
498  }
499  break;
500  case PBTYPE_SINGLEANION:
501  sumCat = 0.0;
502  sumAnion = 0.0;
503  for (k = 0; k < m_kk; k++) {
504  moleFractionsTmp_[k] = moleFractions_[k];
505  }
506  kMax = npos;
507  sumMax = 0.0;
508  for (k = 0; k < cationList_.size(); k++) {
509  kCat = cationList_[k];
510  chP = m_speciesCharge[kCat];
511  if (moleFractions_[kCat] > sumMax) {
512  kMax = k;
513  sumMax = moleFractions_[kCat];
514  }
515  sumCat += chP * moleFractions_[kCat];
516  }
517  k = anionList_[0];
518  chM = m_speciesCharge[k];
519  sumAnion = moleFractions_[k] * chM;
520  sum = sumCat - sumAnion;
521  if (fabs(sum) > 1.0E-16) {
522  moleFractionsTmp_[cationList_[kMax]] -= sum / m_speciesCharge[kMax];
523  sum = 0.0;
524  for (k = 0; k < numCationSpecies_; k++) {
525  sum += moleFractionsTmp_[k];
526  }
527  for (k = 0; k < numCationSpecies_; k++) {
528  moleFractionsTmp_[k]/= sum;
529  }
530  }
531 
532  for (k = 0; k < numCationSpecies_; k++) {
533  PBMoleFractions_[k] = moleFractionsTmp_[cationList_[k]];
534  }
535  for (k = 0; k < numPassThroughSpecies_; k++) {
536  PBMoleFractions_[neutralPBindexStart + k] = moleFractions_[passThroughList_[k]];
537  }
538 
539  sum = std::max(0.0, PBMoleFractions_[0]);
540  for (k = 1; k < numPBSpecies_; k++) {
541  sum += PBMoleFractions_[k];
542 
543  }
544  for (k = 0; k < numPBSpecies_; k++) {
545  PBMoleFractions_[k] /= sum;
546  }
547 
548  break;
549  case PBTYPE_SINGLECATION:
550  throw CanteraError("eosType", "Unknown type");
551 
552  break;
553 
554  case PBTYPE_MULTICATIONANION:
555  throw CanteraError("eosType", "Unknown type");
556 
557  break;
558  default:
559  throw CanteraError("eosType", "Unknown type");
560  break;
561 
562  }
563 }
564 //====================================================================================================================
565 
566 // Update the activity coefficients
567 /*
568  * This function will be called to update the internally stored
569  * natural logarithm of the activity coefficients
570  *
571  */
573 {
574  for (size_t k = 0; k < m_kk; k++) {
575  lnActCoeff_Scaled_[k] = 0.0;
576  }
577 }
578 //====================================================================================================================
580 {
581 
582 
583 }
584 //====================================================================================================================
585 // Internal routine that calculates the derivative of the activity coefficients wrt
586 // the mole fractions.
587 /*
588  * This routine calculates the the derivative of the activity coefficients wrt to mole fraction
589  * with all other mole fractions held constant. This is strictly not permitted. However, if the
590  * resulting matrix is multiplied by a permissible deltaX vector then everything is ok.
591  *
592  * This is the natural way to handle concentration derivatives in this routine.
593  */
595 {
596 
597 }
598 //====================================================================================================================
599 /*
600  * ------------ Partial Molar Properties of the Solution ------------
601  */
602 //====================================================================================================================
603 doublereal MolarityIonicVPSSTP::err(std::string msg) const
604 {
605  throw CanteraError("MolarityIonicVPSSTP","Base class method "
606  +msg+" called. Equation of state type: "+int2str(eosType()));
607  return 0;
608 }
609 //====================================================================================================================
610 /*
611  * @internal Initialize. This method is provided to allow
612  * subclasses to perform any initialization required after all
613  * species have been added. For example, it might be used to
614  * resize internal work arrays that must have an entry for
615  * each species. The base class implementation does nothing,
616  * and subclasses that do not require initialization do not
617  * need to overload this method. When importing a CTML phase
618  * description, this method is called just prior to returning
619  * from function importPhase.
620  *
621  * @see importCTML.cpp
622  */
624 {
626  initLengths();
627  /*
628  * Go find the list of cations and anions
629  */
630  double ch;
631  numCationSpecies_ = 0;
632  cationList_.clear();
633  anionList_.clear();
634  passThroughList_.clear();
635  for (size_t k = 0; k < m_kk; k++) {
636  ch = m_speciesCharge[k];
637  if (ch > 0.0) {
638  cationList_.push_back(k);
640  } else if (ch < 0.0) {
641  anionList_.push_back(k);
642  numAnionSpecies_++;
643  } else {
644  passThroughList_.push_back(k);
645  numPassThroughSpecies_++;
646  }
647  }
648  numPBSpecies_ = numCationSpecies_ + numAnionSpecies_ - 1;
649  neutralPBindexStart = numPBSpecies_;
650  PBType_ = PBTYPE_MULTICATIONANION;
651  if (numAnionSpecies_ == 1) {
652  PBType_ = PBTYPE_SINGLEANION;
653  } else if (numCationSpecies_ == 1) {
654  PBType_ = PBTYPE_SINGLECATION;
655  }
656  if (numAnionSpecies_ == 0 && numCationSpecies_ == 0) {
657  PBType_ = PBTYPE_PASSTHROUGH;
658  }
659 }
660 //====================================================================================================================
661 // Initialize lengths of local variables after all species have been identified.
663 {
664  m_kk = nSpecies();
665  moleFractionsTmp_.resize(m_kk);
666 }
667 //====================================================================================================================
668 /*
669  * initThermoXML() (virtual from ThermoPhase)
670  * Import and initialize a ThermoPhase object
671  *
672  * @param phaseNode This object must be the phase node of a
673  * complete XML tree
674  * description of the phase, including all of the
675  * species data. In other words while "phase" must
676  * point to an XML phase object, it must have
677  * sibling nodes "speciesData" that describe
678  * the species in the phase.
679  * @param id ID of the phase. If nonnull, a check is done
680  * to see if phaseNode is pointing to the phase
681  * with the correct id.
682  */
683 void MolarityIonicVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id)
684 {
685  std::string subname = "MolarityIonicVPSSTP::initThermoXML";
686  std::string stemp;
687  /*
688  * Check on the thermo field. Must have:
689  * <thermo model="MolarityIonic" />
690  */
691 
692  XML_Node& thermoNode = phaseNode.child("thermo");
693  std::string mStringa = thermoNode.attrib("model");
694  std::string mString = lowercase(mStringa);
695  if (mString != "molarityionicvpss" && mString != "molarityionicvpsstp") {
696  throw CanteraError(subname.c_str(),
697  "Unknown thermo model: " + mStringa + " - This object only knows \"MolarityIonicVPSSTP\" ");
698  }
699  /*
700  * Go get all of the coefficients and factors in the
701  * activityCoefficients XML block
702  */
703  /*
704  * Go get all of the coefficients and factors in the
705  * activityCoefficients XML block
706  */
707  XML_Node* acNodePtr = 0;
708  if (thermoNode.hasChild("activityCoefficients")) {
709  XML_Node& acNode = thermoNode.child("activityCoefficients");
710  acNodePtr = &acNode;
711  std::string mStringa = acNode.attrib("model");
712  std::string mString = lowercase(mStringa);
713  // if (mString != "redlich-kister") {
714  // throw CanteraError(subname.c_str(),
715  // "Unknown activity coefficient model: " + mStringa);
716  //}
717  size_t n = acNodePtr->nChildren();
718  for (size_t i = 0; i < n; i++) {
719  XML_Node& xmlACChild = acNodePtr->child(i);
720  stemp = xmlACChild.name();
721  std::string nodeName = lowercase(stemp);
722  /*
723  * Process a binary interaction
724  */
725  if (nodeName == "binaryneutralspeciesparameters") {
726  readXMLBinarySpecies(xmlACChild);
727  }
728  }
729  }
730 
731 
732  /*
733  * Go down the chain
734  */
735  GibbsExcessVPSSTP::initThermoXML(phaseNode, id);
736 }
737 //====================================================================================================================
738 // Process an XML node called "binaryNeutralSpeciesParameters"
739 /*
740  * This node contains all of the parameters necessary to describe
741  * a single binary interaction. This function reads the XML file and writes the coefficients
742  * it finds to an internal data structures.
743  */
745 {
746  std::string xname = xmLBinarySpecies.name();
747 
748 }
749 //====================================================================================================================
750 /*
751  * Format a summary of the mixture state for output.
752  */
753 std::string MolarityIonicVPSSTP::report(bool show_thermo) const
754 {
755  char p[800];
756  string s = "";
757  try {
758  if (name() != "") {
759  sprintf(p, " \n %s:\n", name().c_str());
760  s += p;
761  }
762  sprintf(p, " \n temperature %12.6g K\n", temperature());
763  s += p;
764  sprintf(p, " pressure %12.6g Pa\n", pressure());
765  s += p;
766  sprintf(p, " density %12.6g kg/m^3\n", density());
767  s += p;
768  sprintf(p, " mean mol. weight %12.6g amu\n", meanMolecularWeight());
769  s += p;
770 
771  doublereal phi = electricPotential();
772  sprintf(p, " potential %12.6g V\n", phi);
773  s += p;
774 
775  size_t kk = nSpecies();
776  vector_fp x(kk);
777  vector_fp molal(kk);
778  vector_fp mu(kk);
779  vector_fp muss(kk);
780  vector_fp acMolal(kk);
781  vector_fp actMolal(kk);
782  getMoleFractions(&x[0]);
783 
784  getChemPotentials(&mu[0]);
785  getStandardChemPotentials(&muss[0]);
786  getActivities(&actMolal[0]);
787 
788 
789  if (show_thermo) {
790  sprintf(p, " \n");
791  s += p;
792  sprintf(p, " 1 kg 1 kmol\n");
793  s += p;
794  sprintf(p, " ----------- ------------\n");
795  s += p;
796  sprintf(p, " enthalpy %12.6g %12.4g J\n",
798  s += p;
799  sprintf(p, " internal energy %12.6g %12.4g J\n",
801  s += p;
802  sprintf(p, " entropy %12.6g %12.4g J/K\n",
804  s += p;
805  sprintf(p, " Gibbs function %12.6g %12.4g J\n",
806  gibbs_mass(), gibbs_mole());
807  s += p;
808  sprintf(p, " heat capacity c_p %12.6g %12.4g J/K\n",
809  cp_mass(), cp_mole());
810  s += p;
811  try {
812  sprintf(p, " heat capacity c_v %12.6g %12.4g J/K\n",
813  cv_mass(), cv_mole());
814  s += p;
815  } catch (CanteraError& err) {
816  err.save();
817  sprintf(p, " heat capacity c_v <not implemented> \n");
818  s += p;
819  }
820  }
821 
822  } catch (CanteraError& err) {
823  err.save();
824  }
825  return s;
826 }
827 //====================================================================================================================
828 }
829