Cantera  2.0
SpeciesThermoFactory.cpp
Go to the documentation of this file.
1 /**
2  * @file SpeciesThermoFactory.cpp
3  * Definitions for factory to build instances of classes that manage the
4  * standard-state thermodynamic properties of a set of species
5  * (see \ref spthermo and class \link Cantera::SpeciesThermoFactory SpeciesThermoFactory\endlink);
6  */
7 // Copyright 2001 California Institute of Technology
8 
10 using namespace std;
11 
13 #include "NasaThermo.h"
14 #include "ShomateThermo.h"
17 #include "cantera/thermo/Mu0Poly.h"
20 
24 #include "cantera/thermo/VPSSMgr.h"
26 
27 #include "cantera/base/xml.h"
28 #include "cantera/base/ctml.h"
29 
30 #include <cmath>
31 
32 using namespace ctml;
33 
34 namespace Cantera
35 {
36 
37 SpeciesThermoFactory* SpeciesThermoFactory::s_factory = 0;
38 mutex_t SpeciesThermoFactory::species_thermo_mutex;
39 
40 //! Examine the types of species thermo parameterizations,
41 //! and return a flag indicating the type of reference state thermo manager
42 //! that will be needed in order to evaluate them all.
43 /*!
44  *
45  * @param spDataNodeList This vector contains a list
46  * of species XML nodes that will be in the phase
47  * @param has_nasa Return int that indicates whether the phase has a NASA polynomial form for one of its species
48  * @param has_shomate Return int that indicates whether the phase has a SHOMATE polynomial form for one of its species
49  * @param has_simple Return int that indicates whether the phase has a SIMPLE polynomial form for one of its species
50  * @param has_other Return int that indicates whether the phase has a form for one of its species that is not one of the ones listed above.
51  *
52  * @todo Make sure that spDadta_node is species Data XML node by checking its name is speciesData
53  */
54 static void getSpeciesThermoTypes(std::vector<XML_Node*> & spDataNodeList,
55  int& has_nasa, int& has_shomate, int& has_simple,
56  int& has_other)
57 {
58  size_t ns = spDataNodeList.size();
59  for (size_t n = 0; n < ns; n++) {
60  XML_Node* spNode = spDataNodeList[n];
61  if (spNode->hasChild("standardState")) {
62  const XML_Node& ss = spNode->child("standardState");
63  string mname = ss["model"];
64  if (mname == "water" || mname == "waterIAPWS") {
65  has_other = 1;
66  continue;
67  }
68  }
69  if (spNode->hasChild("thermo")) {
70  const XML_Node& th = spNode->child("thermo");
71  if (th.hasChild("NASA")) {
72  has_nasa = 1;
73  } else if (th.hasChild("Shomate")) {
74  has_shomate = 1;
75  } else if (th.hasChild("MinEQ3")) {
76  has_shomate = 1;
77  } else if (th.hasChild("const_cp")) {
78  has_simple = 1;
79  } else if (th.hasChild("poly")) {
80  if (th.child("poly")["order"] == "1") {
81  has_simple = 1;
82  } else throw CanteraError("newSpeciesThermo",
83  "poly with order > 1 not yet supported");
84  } else if (th.hasChild("Mu0")) {
85  has_other = 1;
86  } else if (th.hasChild("NASA9")) {
87  has_other = 1;
88  } else if (th.hasChild("NASA9MULTITEMP")) {
89  has_other = 1;
90  } else if (th.hasChild("adsorbate")) {
91  has_other = 1;
92  } else {
93  has_other = 1;
94  //throw UnknownSpeciesThermoModel("getSpeciesThermoTypes:",
95  // spNode->attrib("name"), "missing");
96  }
97  } else {
98  throw CanteraError("getSpeciesThermoTypes:",
99  spNode->attrib("name") + " is missing the thermo XML node");
100  }
101  }
102 }
103 
104 //! Static method to return an instance of this class
105 /*!
106  * This class is implemented as a singleton -- one in which
107  * only one instance is needed. The recommended way to access
108  * the factory is to call this static method, which
109  * instantiates the class if it is the first call, but
110  * otherwise simply returns the pointer to the existing
111  * instance.
112  */
113 SpeciesThermoFactory* SpeciesThermoFactory::factory()
114 {
115  ScopedLock lock(species_thermo_mutex);
116  if (!s_factory) {
117  s_factory = new SpeciesThermoFactory;
118  }
119  return s_factory;
120 }
121 
122 // Delete static instance of this class
123 /*
124  * If it is necessary to explicitly delete the factory before
125  * the process terminates (for example, when checking for
126  * memory leaks) then this method can be called to delete it.
127  */
128 void SpeciesThermoFactory::deleteFactory()
129 {
130  ScopedLock lock(species_thermo_mutex);
131  if (s_factory) {
132  delete s_factory;
133  s_factory = 0;
134  }
135 }
136 
137 // Destructor
138 /*
139  * Doesn't do anything. We do not delete statically
140  * created single instance of this class here, because it would
141  * create an infinite loop if destructor is called for that
142  * single instance.
143  */
144 SpeciesThermoFactory::~SpeciesThermoFactory()
145 {
146 }
147 
148 /*
149  * Return a species thermo manager to handle the parameterizations
150  * specified in a CTML phase specification.
151  */
152 SpeciesThermo* SpeciesThermoFactory::newSpeciesThermo(std::vector<XML_Node*> & spDataNodeList) const
153 {
154  int inasa = 0, ishomate = 0, isimple = 0, iother = 0;
155  try {
156  getSpeciesThermoTypes(spDataNodeList, inasa, ishomate, isimple, iother);
157  } catch (UnknownSpeciesThermoModel) {
158  iother = 1;
159  popError();
160  }
161  if (iother) {
162  //writelog("returning new GeneralSpeciesThermo");
163  return new GeneralSpeciesThermo();
164  }
165  return newSpeciesThermo(NASA*inasa
166  + SHOMATE*ishomate + SIMPLE*isimple);
167 }
168 
169 
170 /*
171  * @todo is this used?
172  */
173 SpeciesThermo* SpeciesThermoFactory::
174 newSpeciesThermoOpt(std::vector<XML_Node*> & spDataNodeList) const
175 {
176  int inasa = 0, ishomate = 0, isimple = 0, iother = 0;
177  try {
178  getSpeciesThermoTypes(spDataNodeList, inasa, ishomate, isimple, iother);
179  } catch (UnknownSpeciesThermoModel) {
180  iother = 1;
181  popError();
182  }
183 
184  if (iother) {
185  return new GeneralSpeciesThermo();
186  }
187  return newSpeciesThermo(NASA*inasa
188  + SHOMATE*ishomate + SIMPLE*isimple);
189 }
190 
191 SpeciesThermo* SpeciesThermoFactory::newSpeciesThermo(int type) const
192 {
193  switch (type) {
194  case NASA:
195  return new NasaThermo;
196  case SHOMATE:
197  return new ShomateThermo;
198  case SIMPLE:
199  return new SimpleThermo;
200  case NASA + SHOMATE:
202  case NASA + SIMPLE:
204  case SHOMATE + SIMPLE:
206  default:
207  throw UnknownSpeciesThermo("SpeciesThermoFactory::newSpeciesThermo",
208  type);
209  return 0;
210  }
211 }
212 
213 SpeciesThermo* SpeciesThermoFactory::newSpeciesThermoManager(std::string& stype) const
214 {
215  std::string ltype = lowercase(stype);
216  if (ltype == "nasa") {
217  return new NasaThermo;
218  } else if (ltype == "shomate") {
219  return new ShomateThermo;
220  } else if (ltype == "simple" || ltype == "constant_cp") {
221  return new SimpleThermo;
222  } else if (ltype == "nasa_shomate_duo") {
224  } else if (ltype == "nasa_simple_duo") {
226  } else if (ltype == "shomate_simple_duo") {
228  } else if (ltype == "general") {
229  return new GeneralSpeciesThermo();
230  } else if (ltype == "") {
231  return (SpeciesThermo*) 0;
232  } else {
233  throw UnknownSpeciesThermo("SpeciesThermoFactory::newSpeciesThermoManager",
234  stype);
235  }
236  return (SpeciesThermo*) 0;
237 }
238 
239 /*
240  * Check the continuity of properties at the midpoint
241  * temperature.
242  */
243 void NasaThermo::checkContinuity(std::string name, double tmid, const doublereal* clow,
244  doublereal* chigh)
245 {
246  // heat capacity
247  doublereal cplow = poly4(tmid, clow);
248  doublereal cphigh = poly4(tmid, chigh);
249  doublereal delta = cplow - cphigh;
250  if (fabs(delta/(fabs(cplow)+1.0E-4)) > 0.001) {
251  writelog("\n\n**** WARNING ****\nFor species "+name+
252  ", discontinuity in cp/R detected at Tmid = "
253  +fp2str(tmid)+"\n");
254  writelog("\tValue computed using low-temperature polynomial: "
255  +fp2str(cplow)+".\n");
256  writelog("\tValue computed using high-temperature polynomial: "
257  +fp2str(cphigh)+".\n");
258  }
259 
260  // enthalpy
261  doublereal hrtlow = enthalpy_RT(tmid, clow);
262  doublereal hrthigh = enthalpy_RT(tmid, chigh);
263  delta = hrtlow - hrthigh;
264  if (fabs(delta/(fabs(hrtlow)+cplow*tmid)) > 0.001) {
265  writelog("\n\n**** WARNING ****\nFor species "+name+
266  ", discontinuity in h/RT detected at Tmid = "
267  +fp2str(tmid)+"\n");
268  writelog("\tValue computed using low-temperature polynomial: "
269  +fp2str(hrtlow)+".\n");
270  writelog("\tValue computed using high-temperature polynomial: "
271  +fp2str(hrthigh)+".\n");
272  }
273 
274  // entropy
275  doublereal srlow = entropy_R(tmid, clow);
276  doublereal srhigh = entropy_R(tmid, chigh);
277  delta = srlow - srhigh;
278  if (fabs(delta/(fabs(srlow)+cplow)) > 0.001) {
279  writelog("\n\n**** WARNING ****\nFor species "+name+
280  ", discontinuity in s/R detected at Tmid = "
281  +fp2str(tmid)+"\n");
282  writelog("\tValue computed using low-temperature polynomial: "
283  +fp2str(srlow)+".\n");
284  writelog("\tValue computed using high-temperature polynomial: "
285  +fp2str(srhigh)+".\n");
286  }
287 }
288 
289 
290 //! Install a NASA polynomial thermodynamic property parameterization for species k into a SpeciesThermo instance.
291 /*!
292  * This is called by method installThermoForSpecies if a NASA block is found in the XML input.
293  *
294  * @param speciesName String name of the species
295  * @param sp SpeciesThermo object that will receive the nasa polynomial object
296  * @param k Species index within the phase
297  * @param f0ptr Ptr to the first XML_Node for the first NASA polynomial
298  * @param f1ptr Ptr to the first XML_Node for the first NASA polynomial
299  */
300 static void installNasaThermoFromXML(std::string speciesName,
301  SpeciesThermo& sp, size_t k,
302  const XML_Node* f0ptr, const XML_Node* f1ptr)
303 {
304  doublereal tmin0, tmax0, tmin1, tmax1, tmin, tmid, tmax;
305 
306  const XML_Node& f0 = *f0ptr;
307 
308  // default to a single temperature range
309  bool dualRange = false;
310 
311  // but if f1ptr is suppled, then it is a two-range
312  // parameterization
313  if (f1ptr) {
314  dualRange = true;
315  }
316 
317  tmin0 = fpValue(f0["Tmin"]);
318  tmax0 = fpValue(f0["Tmax"]);
319 
320  doublereal p0 = OneAtm;
321  if (f0.hasAttrib("P0")) {
322  p0 = fpValue(f0["P0"]);
323  }
324  if (f0.hasAttrib("Pref")) {
325  p0 = fpValue(f0["Pref"]);
326  }
327  p0 = OneAtm;
328 
329  tmin1 = tmax0;
330  tmax1 = tmin1 + 0.0001;
331  if (dualRange) {
332  tmin1 = fpValue((*f1ptr)["Tmin"]);
333  tmax1 = fpValue((*f1ptr)["Tmax"]);
334  }
335 
336  vector_fp c0, c1;
337  if (fabs(tmax0 - tmin1) < 0.01) {
338  // f0 has the lower T data, and f1 the higher T data
339  tmin = tmin0;
340  tmid = tmax0;
341  tmax = tmax1;
342  getFloatArray(f0.child("floatArray"), c0, false);
343  if (dualRange) {
344  getFloatArray(f1ptr->child("floatArray"), c1, false);
345  } else {
346  // if there is no higher range data, then copy c0 to c1.
347  c1.resize(7,0.0);
348  copy(c0.begin(), c0.end(), c1.begin());
349  }
350  } else if (fabs(tmax1 - tmin0) < 0.01) {
351  // f1 has the lower T data, and f0 the higher T data
352  tmin = tmin1;
353  tmid = tmax1;
354  tmax = tmax0;
355  getFloatArray(f1ptr->child("floatArray"), c0, false);
356  getFloatArray(f0.child("floatArray"), c1, false);
357  } else {
358  throw CanteraError("installNasaThermo",
359  "non-continuous temperature ranges.");
360  }
361 
362  // The NasaThermo species property manager expects the
363  // coefficients in a different order, so rearrange them.
364  vector_fp c(15);
365  c[0] = tmid;
366 
367  c[1] = c0[5];
368  c[2] = c0[6];
369  copy(c0.begin(), c0.begin()+5, c.begin() + 3);
370  c[8] = c1[5];
371  c[9] = c1[6];
372  copy(c1.begin(), c1.begin()+5, c.begin() + 10);
373  sp.install(speciesName, k, NASA, &c[0], tmin, tmax, p0);
374 }
375 
376 //! Look up the elemental reference state entropies
377 /*!
378  * @param elemName String name of the element
379  * @param th_ptr Pointer to the thermophase.
380  */
381 static doublereal LookupGe(const std::string& elemName, ThermoPhase* th_ptr)
382 {
383  size_t iE = th_ptr->elementIndex(elemName);
384  if (iE == npos) {
385  throw CanteraError("PDSS_HKFT::LookupGe", "element " + elemName + " not found");
386  }
387  doublereal geValue = th_ptr->entropyElement298(iE);
388  if (geValue == ENTROPY298_UNKNOWN) {
389  throw CanteraError("PDSS_HKFT::LookupGe",
390  "element " + elemName + " doesn not have a supplied entropy298");
391  }
392  geValue *= (-298.15);
393  return geValue;
394 }
395 
396 //! Convert delta G formulation
397 /*!
398  * Calculates the sum of the elemental reference state entropies
399  *
400  * @param k species index
401  * @param th_ptr Pointer to the ThermoPhase
402  */
403 static doublereal convertDGFormation(size_t k, ThermoPhase* th_ptr)
404 {
405  /*
406  * Ok let's get the element compositions and conversion factors.
407  */
408  size_t ne = th_ptr->nElements();
409  doublereal na;
410  doublereal ge;
411  string ename;
412 
413  doublereal totalSum = 0.0;
414  for (size_t m = 0; m < ne; m++) {
415  na = th_ptr->nAtoms(k, m);
416  if (na > 0.0) {
417  ename = th_ptr->elementName(m);
418  ge = LookupGe(ename, th_ptr);
419  totalSum += na * ge;
420  }
421  }
422  return totalSum;
423 }
424 
425 
426 //! Install a NASA96 polynomial thermodynamic property parameterization for species k into a SpeciesThermo instance.
427 /*!
428  * This is called by method installThermoForSpecies if a MinEQ3node block is found in the XML input.
429  *
430  * @param speciesName String name of the species
431  * @param th_ptr Pointer to the %ThermoPhase object
432  * @param sp SpeciesThermo object that will receive the nasa polynomial object
433  * @param k Species index within the phase
434  * @param MinEQ3node Ptr to the first XML_Node for the first MinEQ3 parameterization
435  */
436 static void installMinEQ3asShomateThermoFromXML(std::string speciesName,
437  ThermoPhase* th_ptr,
438  SpeciesThermo& sp, size_t k,
439  const XML_Node* MinEQ3node)
440 {
441 
442  vector_fp coef(15), c0(7, 0.0);
443  std::string astring = (*MinEQ3node)["Tmin"];
444  doublereal tmin0 = strSItoDbl(astring);
445  astring = (*MinEQ3node)["Tmax"];
446  doublereal tmax0 = strSItoDbl(astring);
447  astring = (*MinEQ3node)["Pref"];
448  doublereal p0 = strSItoDbl(astring);
449 
450  doublereal deltaG_formation_pr_tr =
451  getFloatDefaultUnits(*MinEQ3node, "DG0_f_Pr_Tr", "cal/gmol", "actEnergy");
452  doublereal deltaH_formation_pr_tr =
453  getFloatDefaultUnits(*MinEQ3node, "DH0_f_Pr_Tr", "cal/gmol", "actEnergy");
454  doublereal Entrop_pr_tr = getFloatDefaultUnits(*MinEQ3node, "S0_Pr_Tr", "cal/gmol/K");
455  doublereal a = getFloatDefaultUnits(*MinEQ3node, "a", "cal/gmol/K");
456  doublereal b = getFloatDefaultUnits(*MinEQ3node, "b", "cal/gmol/K2");
457  doublereal c = getFloatDefaultUnits(*MinEQ3node, "c", "cal-K/gmol");
458  doublereal dg = deltaG_formation_pr_tr * 4.184 * 1.0E3;
459  doublereal fac = convertDGFormation(k, th_ptr);
460  doublereal Mu0_tr_pr = fac + dg;
461  doublereal e = Entrop_pr_tr * 1.0E3 * 4.184;
462  doublereal Hcalc = Mu0_tr_pr + 298.15 * e;
463  doublereal DHjmol = deltaH_formation_pr_tr * 1.0E3 * 4.184;
464 
465  // If the discrepency is greater than 100 cal gmol-1, print
466  // an error and exit.
467  if (fabs(Hcalc -DHjmol) > 10.* 1.0E6 * 4.184) {
468  throw CanteraError("installMinEQ3asShomateThermoFromXML()",
469  "DHjmol is not consistent with G and S" +
470  fp2str(Hcalc) + " vs " + fp2str(DHjmol));
471  }
472 
473  /*
474  * Now calculate the shomate polynomials
475  *
476  * Cp first
477  *
478  * Shomate: (Joules / gmol / K)
479  * Cp = As + Bs * t + Cs * t*t + Ds * t*t*t + Es / (t*t)
480  * where
481  * t = temperature(Kelvin) / 1000
482  */
483  double As = a * 4.184;
484  double Bs = b * 4.184 * 1000.;
485  double Cs = 0.0;
486  double Ds = 0.0;
487  double Es = c * 4.184 / (1.0E6);
488 
489  double t = 298.15 / 1000.;
490  double H298smFs = As * t + Bs * t * t / 2.0 - Es / t;
491 
492  double HcalcS = Hcalc / 1.0E6;
493  double Fs = HcalcS - H298smFs;
494 
495  double S298smGs = As * log(t) + Bs * t - Es/(2.0*t*t);
496  double ScalcS = e / 1.0E3;
497  double Gs = ScalcS - S298smGs;
498 
499  c0[0] = As;
500  c0[1] = Bs;
501  c0[2] = Cs;
502  c0[3] = Ds;
503  c0[4] = Es;
504  c0[5] = Fs;
505  c0[6] = Gs;
506 
507  coef[0] = tmax0 - 0.001;
508  copy(c0.begin(), c0.begin()+7, coef.begin() + 1);
509  copy(c0.begin(), c0.begin()+7, coef.begin() + 8);
510  sp.install(speciesName, k, SHOMATE, &coef[0], tmin0, tmax0, p0);
511 }
512 
513 //! Install a Shomate polynomial thermodynamic property parameterization for species k into a SpeciesThermo instance.
514 /*!
515  * This is called by method installThermoForSpecies if a Shomate block is found in the XML input.
516  *
517  * @param speciesName String name of the species
518  * @param sp SpeciesThermo object that will receive the nasa polynomial object
519  * @param k Species index within the phase
520  * @param f0ptr Ptr to the first XML_Node for the first NASA polynomial
521  * @param f1ptr Ptr to the first XML_Node for the first NASA polynomial
522  */
523 static void installShomateThermoFromXML(std::string speciesName, SpeciesThermo& sp, size_t k,
524  const XML_Node* f0ptr, const XML_Node* f1ptr)
525 {
526  doublereal tmin0, tmax0, tmin1, tmax1, tmin, tmid, tmax;
527  const XML_Node& f0 = *f0ptr;
528  bool dualRange = false;
529  if (f1ptr) {
530  dualRange = true;
531  }
532  tmin0 = fpValue(f0["Tmin"]);
533  tmax0 = fpValue(f0["Tmax"]);
534 
535  doublereal p0 = OneAtm;
536  if (f0.hasAttrib("P0")) {
537  p0 = fpValue(f0["P0"]);
538  }
539  if (f0.hasAttrib("Pref")) {
540  p0 = fpValue(f0["Pref"]);
541  }
542  p0 = OneAtm;
543 
544  tmin1 = tmax0;
545  tmax1 = tmin1 + 0.0001;
546  if (dualRange) {
547  tmin1 = fpValue((*f1ptr)["Tmin"]);
548  tmax1 = fpValue((*f1ptr)["Tmax"]);
549  }
550 
551  vector_fp c0, c1;
552  if (fabs(tmax0 - tmin1) < 0.01) {
553  tmin = tmin0;
554  tmid = tmax0;
555  tmax = tmax1;
556  getFloatArray(f0.child("floatArray"), c0, false);
557  if (dualRange) {
558  getFloatArray(f1ptr->child("floatArray"), c1, false);
559  } else {
560  c1.resize(7,0.0);
561  copy(c0.begin(), c0.begin()+7, c1.begin());
562  }
563  } else if (fabs(tmax1 - tmin0) < 0.01) {
564  tmin = tmin1;
565  tmid = tmax1;
566  tmax = tmax0;
567  getFloatArray(f1ptr->child("floatArray"), c0, false);
568  getFloatArray(f0.child("floatArray"), c1, false);
569  } else {
570  throw CanteraError("installShomateThermoFromXML",
571  "non-continuous temperature ranges.");
572  }
573  vector_fp c(15);
574  c[0] = tmid;
575  copy(c0.begin(), c0.begin()+7, c.begin() + 1);
576  copy(c1.begin(), c1.begin()+7, c.begin() + 8);
577  sp.install(speciesName, k, SHOMATE, &c[0], tmin, tmax, p0);
578 }
579 
580 //! Install a Simple thermodynamic property parameterization for species k into a SpeciesThermo instance.
581 /*!
582  * This is called by method installThermoForSpecies if a SimpleThermo block is found
583  *
584  * @param speciesName String name of the species
585  * @param sp SpeciesThermo object that will receive the nasa polynomial object
586  * @param k Species index within the phase
587  * @param f XML_Node for the SimpleThermo block
588  */
589 static void installSimpleThermoFromXML(std::string speciesName,
590  SpeciesThermo& sp, size_t k,
591  const XML_Node& f)
592 {
593  doublereal tmin, tmax;
594  tmin = fpValue(f["Tmin"]);
595  tmax = fpValue(f["Tmax"]);
596  if (tmax == 0.0) {
597  tmax = 1.0e30;
598  }
599 
600  vector_fp c(4);
601  c[0] = getFloat(f, "t0", "toSI");
602  c[1] = getFloat(f, "h0", "toSI");
603  c[2] = getFloat(f, "s0", "toSI");
604  c[3] = getFloat(f, "cp0", "toSI");
605  doublereal p0 = OneAtm;
606  sp.install(speciesName, k, SIMPLE, &c[0], tmin, tmax, p0);
607 }
608 
609 
610 //! Install a NASA9 polynomial thermodynamic property parameterization for species k into a SpeciesThermo instance.
611 /*!
612  * This is called by method installThermoForSpecies if a NASA9 block is found in the XML input.
613  *
614  * @param speciesName String name of the species
615  * @param sp SpeciesThermo object that will receive the nasa polynomial object
616  * @param k Species index within the phase
617  * @param tp Vector of XML Nodes that make up the parameterization
618  */
619 static void installNasa9ThermoFromXML(std::string speciesName,
620  SpeciesThermo& sp, size_t k,
621  const std::vector<XML_Node*>& tp)
622 {
623  const XML_Node* fptr = tp[0];
624  int nRegions = 0;
625  vector_fp cPoly;
626  Nasa9Poly1* np_ptr = 0;
627  std::vector<Nasa9Poly1*> regionPtrs;
628  doublereal tmin, tmax, pref = OneAtm;
629  // Loop over all of the possible temperature regions
630  for (size_t i = 0; i < tp.size(); i++) {
631  fptr = tp[i];
632  if (fptr) {
633  if (fptr->name() == "NASA9") {
634  if (fptr->hasChild("floatArray")) {
635 
636  tmin = fpValue((*fptr)["Tmin"]);
637  tmax = fpValue((*fptr)["Tmax"]);
638  if ((*fptr).hasAttrib("P0")) {
639  pref = fpValue((*fptr)["P0"]);
640  }
641  if ((*fptr).hasAttrib("Pref")) {
642  pref = fpValue((*fptr)["Pref"]);
643  }
644 
645  getFloatArray(fptr->child("floatArray"), cPoly, false);
646  if (cPoly.size() != 9) {
647  throw CanteraError("installNasa9ThermoFromXML",
648  "Expected 9 coeff polynomial");
649  }
650  np_ptr = new Nasa9Poly1(k, tmin, tmax, pref,
651  DATA_PTR(cPoly));
652  regionPtrs.push_back(np_ptr);
653  nRegions++;
654  }
655  }
656  }
657  }
658  if (nRegions == 0) {
659  throw UnknownSpeciesThermoModel("installThermoForSpecies",
660  speciesName, " ");
661  } else if (nRegions == 1) {
662  sp.install_STIT(np_ptr);
663  } else {
664  Nasa9PolyMultiTempRegion* npMulti_ptr = new Nasa9PolyMultiTempRegion(regionPtrs);
665  sp.install_STIT(npMulti_ptr);
666  }
667 }
668 
669 
670 //! Install a Adsorbate polynomial thermodynamic property parameterization for species k into a SpeciesThermo instance.
671 /*!
672  * This is called by method installThermoForSpecies if a Adsorbate block is found in the XML input.
673  *
674  * @param speciesName String name of the species
675  * @param sp SpeciesThermo object that will receive the nasa polynomial object
676  * @param k Species index within the phase
677  * @param f XML Node that contains the parameterization
678  */
679 static void installAdsorbateThermoFromXML(std::string speciesName,
680  SpeciesThermo& sp, size_t k,
681  const XML_Node& f)
682 {
683  vector_fp freqs;
684  doublereal tmin, tmax, pref = OneAtm;
685  size_t nfreq = 0;
686  tmin = fpValue(f["Tmin"]);
687  tmax = fpValue(f["Tmax"]);
688  if (f.hasAttrib("P0")) {
689  pref = fpValue(f["P0"]);
690  }
691  if (f.hasAttrib("Pref")) {
692  pref = fpValue(f["Pref"]);
693  }
694  if (tmax == 0.0) {
695  tmax = 1.0e30;
696  }
697 
698  if (f.hasChild("floatArray")) {
699  getFloatArray(f.child("floatArray"), freqs, false);
700  nfreq = freqs.size();
701  }
702  for (size_t n = 0; n < nfreq; n++) {
703  freqs[n] *= 3.0e10;
704  }
705  vector_fp coeffs(nfreq + 2);
706  coeffs[0] = static_cast<double>(nfreq);
707  coeffs[1] = getFloat(f, "binding_energy", "toSI");
708  copy(freqs.begin(), freqs.end(), coeffs.begin() + 2);
709  (&sp)->install(speciesName, k, ADSORBATE, &coeffs[0], tmin, tmax, pref);
710 }
711 
712 //================================================================================================
713 // Install a species thermodynamic property parameterization
714 // for the reference state for one species into a species thermo manager.
715 /*
716  * @param k Species number
717  * @param speciesNode Reference to the XML node specifying the species standard
718  * state information
719  * @param th_ptr Pointer to the %ThermoPhase object for the species
720  * @param spthermo Species reference state thermo manager
721  * @param phaseNode_ptr Optional pointer to the XML phase
722  * information for the phase in which the species
723  * resides
724  */
725 void SpeciesThermoFactory::installThermoForSpecies
726 (size_t k, const XML_Node& speciesNode, ThermoPhase* th_ptr,
727  SpeciesThermo& spthermo, const XML_Node* phaseNode_ptr) const
728 {
729  /*
730  * Check to see that the species block has a thermo block
731  * before processing. Throw an error if not there.
732  */
733  if (!(speciesNode.hasChild("thermo"))) {
734  throw UnknownSpeciesThermoModel("installThermoForSpecies",
735  speciesNode["name"], "<nonexistent>");
736  }
737  const XML_Node& thermo = speciesNode.child("thermo");
738 
739  // Get the children of the thermo XML node. In the next bit of code we take out the comments that
740  // may have been children of the thermo XML node by doing a selective copy.
741  // These shouldn't interfere with the algorithm at any point.
742  const std::vector<XML_Node*>& tpWC = thermo.children();
743  std::vector<XML_Node*> tp;
744  for (int i = 0; i < static_cast<int>(tpWC.size()); i++) {
745  if (!(tpWC[i])->isComment()) {
746  tp.push_back(tpWC[i]);
747  }
748  }
749  int nc = static_cast<int>(tp.size());
750  string mname = thermo["model"];
751 
752  if (mname == "MineralEQ3") {
753  const XML_Node* f = tp[0];
754  if (f->name() != "MinEQ3") {
755  throw CanteraError("SpeciesThermoFactory::installThermoForSpecies",
756  "confused: expedted MinEQ3");
757  }
758  installMinEQ3asShomateThermoFromXML(speciesNode["name"], th_ptr, spthermo, k, f);
759  } else {
760  if (nc == 1) {
761  const XML_Node* f = tp[0];
762  if (f->name() == "Shomate") {
763  installShomateThermoFromXML(speciesNode["name"], spthermo, k, f, 0);
764  } else if (f->name() == "const_cp") {
765  installSimpleThermoFromXML(speciesNode["name"], spthermo, k, *f);
766  } else if (f->name() == "NASA") {
767  installNasaThermoFromXML(speciesNode["name"], spthermo, k, f, 0);
768  } else if (f->name() == "Mu0") {
769  installMu0ThermoFromXML(speciesNode["name"], spthermo, k, f);
770  } else if (f->name() == "NASA9") {
771  installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp);
772  }
773  else if (f->name() == "adsorbate") {
774  installAdsorbateThermoFromXML(speciesNode["name"], spthermo, k, *f);
775  }
776  else {
777  throw UnknownSpeciesThermoModel("installThermoForSpecies",
778  speciesNode["name"], f->name());
779  }
780  } else if (nc == 2) {
781  const XML_Node* f0 = tp[0];
782  const XML_Node* f1 = tp[1];
783  if (f0->name() == "NASA" && f1->name() == "NASA") {
784  installNasaThermoFromXML(speciesNode["name"], spthermo, k, f0, f1);
785  } else if (f0->name() == "Shomate" && f1->name() == "Shomate") {
786  installShomateThermoFromXML(speciesNode["name"], spthermo, k, f0, f1);
787  } else if (f0->name() == "NASA9" && f1->name() == "NASA9") {
788  installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp);
789  } else {
790  throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"],
791  f0->name() + " and " + f1->name());
792  }
793  } else if (nc > 2) {
794  const XML_Node* f0 = tp[0];
795  if (f0->name() == "NASA9") {
796  installNasa9ThermoFromXML(speciesNode["name"], spthermo, k, tp);
797  } else {
798  throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"],
799  "multiple");
800  }
801  } else {
802  throw UnknownSpeciesThermoModel("installThermoForSpecies", speciesNode["name"],
803  "multiple");
804  }
805  }
806 }
807 //================================================================================================
808 // Install a species thermodynamic property parameterization
809 // for the standard state for one species into a species thermo manager, VPSSMgr
810 /*
811  * This is a wrapper around the createInstallVPSS() function in the
812  * VPStandardStateTP object.
813  *
814  * This serves to install the species into vpss_ptr, create a PDSS file. We also
815  * read the xml database to extract the constants for these steps.
816  *
817  * @param k species number
818  * @param speciesNode Reference to the XML node specifying the species standard
819  * state information
820  * @param vp_ptr variable pressure ThermoPhase object
821  * @param vpss_ptr Pointer to the Manager for calculating variable pressure
822  * substances.
823  * @param spthermo_ptr Species reference state thermo manager
824  * @param phaseNode_ptr Optional Pointer to the XML phase
825  * information for the phase in which the species
826  * resides
827  */
828 void SpeciesThermoFactory::
829 installVPThermoForSpecies(size_t k, const XML_Node& speciesNode,
830  VPStandardStateTP* vp_ptr,
831  VPSSMgr* vpssmgr_ptr,
832  SpeciesThermo* spthermo_ptr,
833  const XML_Node* phaseNode_ptr) const
834 {
835 
836  // Call the VPStandardStateTP object to install the pressure dependent species
837  // standard state into the object.
838  //
839  // We don't need to pass spthermo_ptr down, because it's already installed
840  // into vp_ptr.
841  //
842  // We don't need to pass vpssmgr_ptr down, because it's already installed
843  // into vp_ptr.
844  vp_ptr->createInstallPDSS(k, speciesNode, phaseNode_ptr);
845 }
846 
847 // Create a new species thermo manager instance, by specifying
848 // the type and (optionally) a pointer to the factory to use to create it.
849 /*
850  * This utility program will look through species nodes. It will discover what
851  * each species needs for its species property managers. Then,
852  * it will malloc and return the proper species property manager to use.
853  *
854  * These functions allow using a different factory class that
855  * derives from SpeciesThermoFactory.
856  *
857  * @param type Species thermo type.
858  * @param f Pointer to a SpeciesThermoFactory. optional parameter.
859  * Defaults to NULL.
860  */
862 {
863  if (f == 0) {
864  f = SpeciesThermoFactory::factory();
865  }
866  SpeciesThermo* sptherm = f->newSpeciesThermo(type);
867  return sptherm;
868 }
869 
870 // Create a new species thermo manager instance, by specifying
871 //the type and (optionally) a pointer to the factory to use to create it.
872 /*
873  * This utility program is a basic factory operation for spawning a
874  * new species reference-state thermo manager
875  *
876  * These functions allows for using a different factory class that
877  * derives from SpeciesThermoFactory. However, no applications of this
878  * have been done yet.
879  *
880  * @param stype String specifying the species thermo type
881  * @param f Pointer to a SpeciesThermoFactory. optional parameter.
882  * Defaults to NULL.
883  */
884 SpeciesThermo* newSpeciesThermoMgr(std::string& stype,
886 {
887  if (f == 0) {
888  f = SpeciesThermoFactory::factory();
889  }
890  SpeciesThermo* sptherm = f->newSpeciesThermoManager(stype);
891  return sptherm;
892 }
893 
894 // Function to return SpeciesThermo manager
895 /*
896  * This utility program will look through species nodes. It will discover what
897  * each species needs for its species property managers. Then,
898  * it will malloc and return the proper species property manager to use.
899  *
900  * These functions allow using a different factory class that
901  * derives from SpeciesThermoFactory.
902  *
903  * @param spData_nodes Vector of XML_Nodes, each of which is a speciesData XML Node.
904  * Each %speciesData node contains a list of XML species elements
905  * e.g., <speciesData id="Species_Data">
906  * @param f Pointer to a SpeciesThermoFactory. optional parameter.
907  * Defaults to NULL.
908  * @param opt Boolean defaults to false.
909  */
910 SpeciesThermo* newSpeciesThermoMgr(std::vector<XML_Node*> spData_nodes,
911  SpeciesThermoFactory* f, bool opt)
912 {
913  if (f == 0) {
914  f = SpeciesThermoFactory::factory();
915  }
916  SpeciesThermo* sptherm;
917  if (opt) {
918  sptherm = f->newSpeciesThermoOpt(spData_nodes);
919  } else {
920  sptherm = f->newSpeciesThermo(spData_nodes);
921  }
922  return sptherm;
923 }
924 
925 }