Cantera  2.0
NasaThermo.h
Go to the documentation of this file.
1 /**
2  * @file NasaThermo.h
3  * Header for the 2 regime 7 coefficient Nasa thermodynamic
4  * polynomials for multiple species in a phase, derived from the
5  * \link Cantera::SpeciesThermo SpeciesThermo\endlink base class (see \ref mgrsrefcalc and
6  * \link Cantera::NasaThermo NasaThermo\endlink).
7  */
8 // Copyright 2003 California Institute of Technology
9 
10 #ifndef CT_NASATHERMO_H
11 #define CT_NASATHERMO_H
12 #include <string>
13 
17 //#include "cantera/numerics/polyfit.h"
18 #include "cantera/base/global.h"
19 
20 namespace Cantera
21 {
22 
23 /**
24  * A species thermodynamic property manager for the NASA
25  * polynomial parameterization with two temperature ranges.
26  *
27  * This class is designed to efficiently evaluate the properties
28  * of a large number of species with the NASA parameterization.
29  *
30  * The original NASA polynomial parameterization expressed the
31  * heat capacity as a fourth-order polynomial in temperature, with
32  * separate coefficients for each of two temperature ranges. (The
33  * newer NASA format adds coefficients for 1/T and 1/T^2, and
34  * allows multiple temperature ranges.) This class is designed for
35  * use with the original parameterization, which is used, for
36  * example, by the Chemkin software package.
37  *
38  * In many cases, the midpoint temperature is the same for many
39  * species. To take advantage of this, class NasaThermo groups
40  * species with a common midpoint temperature, so that checking
41  * which range the desired temperature is in need be done only
42  * once for each group.
43  *
44  * @note There is a special CTML element for entering the
45  * coefficients of this parameterization.
46  * @see importCTML
47  *
48  * @ingroup mgrsrefcalc
49  */
50 class NasaThermo : public SpeciesThermo
51 {
52 
53 public:
54 
55  //! Initialized to the type of parameterization
56  /*!
57  * Note, this value is used in some template functions
58  */
59  const int ID;
60 
61  //! constructor
63  ID(NASA),
64  m_tlow_max(0.0),
65  m_thigh_min(1.e30),
66  m_p0(-1.0),
67  m_ngroups(0) {
68  m_t.resize(6);
69  }
70 
71  //! Copy constructor
72  /*!
73  * @param right NasaThermo object to be copied.
74  */
75  NasaThermo(const NasaThermo& right) :
76  ID(NASA),
77  m_tlow_max(0.0),
78  m_thigh_min(1.e30),
79  m_p0(-1.0),
80  m_ngroups(0) {
81  *this = operator=(right);
82  }
83 
84  //! Assignment operator
85  /*!
86  * @param right NasaThermo object to be copied.
87  */
88  NasaThermo& operator=(const NasaThermo& right) {
89  /*
90  * Check for self assignment.
91  */
92  if (this == &right) {
93  return *this;
94  }
95 
96  m_high = right.m_high;
97  m_low = right.m_low;
98  m_index = right.m_index;
99  m_tmid = right.m_tmid;
100  m_tlow_max = right.m_tlow_max;
101  m_thigh_min = right.m_thigh_min;
102  m_tlow = right.m_tlow;
103  m_thigh = right.m_thigh;
104  m_p0 = right.m_p0;
105  m_ngroups = right.m_ngroups;
106  m_t = right.m_t;
107  m_group_map = right.m_group_map;
109  m_name = right.m_name;
110 
111  return *this;
112  }
113 
114 
115  //! destructor
116  virtual ~NasaThermo() {}
117 
118  //! Duplication routine for objects which inherit from
119  //! %SpeciesThermo
120  /*!
121  * This virtual routine can be used to duplicate %SpeciesThermo objects
122  * inherited from %SpeciesThermo even if the application only has
123  * a pointer to %SpeciesThermo to work with.
124  * ->commented out because we first need to add copy constructors
125  * and assignment operators to all of the derived classes.
126  */
128  NasaThermo* nt = new NasaThermo(*this);
129  return (SpeciesThermo*) nt;
130  }
131 
132  //! install a new species thermodynamic property
133  //! parameterization for one species.
134  /*!
135  *
136  * @param name Name of the species
137  * @param index The 'update' method will update the property
138  * values for this species
139  * at position i index in the property arrays.
140  * @param type int flag specifying the type of parameterization to be
141  * installed.
142  * @param c vector of coefficients for the parameterization.
143  * - c[0] midpoint temperature
144  * - c[1] - c[7] coefficients for low T range
145  * - c[8] - c[14] coefficients for high T range
146  * @param minTemp minimum temperature for which this parameterization
147  * is valid.
148  * @param maxTemp maximum temperature for which this parameterization
149  * is valid.
150  * @param refPressure standard-state pressure for this
151  * parameterization.
152  * @see speciesThermoTypes.h
153  */
154  virtual void install(std::string name, size_t index, int type,
155  const doublereal* c,
156  doublereal minTemp, doublereal maxTemp,
157  doublereal refPressure) {
158 
159  m_name[index] = name;
160  int imid = int(c[0]); // midpoint temp converted to integer
161  int igrp = m_index[imid]; // has this value been seen before?
162  if (igrp == 0) { // if not, prepare new group
163  std::vector<NasaPoly1> v;
164  m_high.push_back(v);
165  m_low.push_back(v);
166  m_tmid.push_back(c[0]);
167  m_index[imid] = igrp = static_cast<int>(m_high.size());
168  m_ngroups++;
169  }
170 
171  m_group_map[index] = igrp;
172  m_posInGroup_map[index] = (int) m_low[igrp-1].size();
173 
174  doublereal tlow = minTemp;
175  doublereal tmid = c[0];
176  doublereal thigh = maxTemp;
177  const doublereal* clow = c + 1;
178 
179  vector_fp chigh(7);
180  copy(c + 8, c + 15, chigh.begin());
181 
182  m_high[igrp-1].push_back(NasaPoly1(index, tmid, thigh,
183  refPressure, &chigh[0]));
184  m_low[igrp-1].push_back(NasaPoly1(index, tlow, tmid,
185  refPressure, clow));
186 
187  vector_fp clu(7), chu(7);
188  clu[5] = clow[0];
189  clu[6] = clow[1];
190  copy(clow+2, clow+7, clu.begin());
191  chu[5] = chigh[0];
192  chu[6] = chigh[1];
193  copy(chigh.begin()+2, chigh.begin()+7, chu.begin());
194 
195  checkContinuity(name, tmid, &clu[0], &chu[0]);
196 
197  if (tlow > m_tlow_max) {
198  m_tlow_max = tlow;
199  }
200  if (thigh < m_thigh_min) {
201  m_thigh_min = thigh;
202  }
203  if (m_tlow.size() < index + 1) {
204  m_tlow.resize(index + 1, tlow);
205  m_thigh.resize(index + 1, thigh);
206  }
207  m_tlow[index] = tlow;
208  m_thigh[index] = thigh;
209  if (m_p0 < 0.0) {
210  m_p0 = refPressure;
211  } else if (fabs(m_p0 - refPressure) > 0.1) {
212  std::string logmsg = " ERROR NasaThermo: New Species, " + name + ", has a different reference pressure, "
213  + fp2str(refPressure) + ", than existing reference pressure, " + fp2str(m_p0) + "\n";
214  writelog(logmsg);
215  logmsg = " This is now a fatal error\n";
216  writelog(logmsg);
217  throw CanteraError("install()", "species have different reference pressures");
218  }
219  m_p0 = refPressure;
220  }
221 
222  //! Install a new species thermodynamic property
223  //! parameterization for one species.
224  /*!
225  * @param stit_ptr Pointer to the SpeciesThermoInterpType object
226  * This will set up the thermo for one species
227  */
228  virtual void install_STIT(SpeciesThermoInterpType* stit_ptr) {
229  throw CanteraError("install_STIT", "not implemented");
230  }
231 
232 
233  //! Like update(), but only updates the single species k.
234  /*!
235  * @param k species index
236  * @param t Temperature (Kelvin)
237  * @param cp_R Vector of Dimensionless heat capacities.
238  * (length m_kk).
239  * @param h_RT Vector of Dimensionless enthalpies.
240  * (length m_kk).
241  * @param s_R Vector of Dimensionless entropies.
242  * (length m_kk).
243  *
244  */
245  virtual void update_one(size_t k, doublereal t, doublereal* cp_R,
246  doublereal* h_RT, doublereal* s_R) const {
247 
248  m_t[0] = t;
249  m_t[1] = t*t;
250  m_t[2] = m_t[1]*t;
251  m_t[3] = m_t[2]*t;
252  m_t[4] = 1.0/t;
253  m_t[5] = log(t);
254 
255  size_t grp = m_group_map[k];
256  size_t pos = m_posInGroup_map[k];
257  const std::vector<NasaPoly1> &mlg = m_low[grp-1];
258  const NasaPoly1* nlow = &(mlg[pos]);
259 
260  doublereal tmid = nlow->maxTemp();
261  if (t < tmid) {
262  nlow->updateProperties(&m_t[0], cp_R, h_RT, s_R);
263  } else {
264  const std::vector<NasaPoly1> &mhg = m_high[grp-1];
265  const NasaPoly1* nhigh = &(mhg[pos]);
266  nhigh->updateProperties(&m_t[0], cp_R, h_RT, s_R);
267  }
268  }
269 
270  //! Compute the reference-state properties for all species.
271  /*!
272  * Given temperature T in K, this method updates the values of
273  * the non-dimensional heat capacity at constant pressure,
274  * enthalpy, and entropy, at the reference pressure, Pref
275  * of each of the standard states.
276  *
277  * @param t Temperature (Kelvin)
278  * @param cp_R Vector of Dimensionless heat capacities.
279  * (length m_kk).
280  * @param h_RT Vector of Dimensionless enthalpies.
281  * (length m_kk).
282  * @param s_R Vector of Dimensionless entropies.
283  * (length m_kk).
284  */
285  virtual void update(doublereal t, doublereal* cp_R,
286  doublereal* h_RT, doublereal* s_R) const {
287  int i;
288 
289  // load functions of temperature into m_t vector
290  m_t[0] = t;
291  m_t[1] = t*t;
292  m_t[2] = m_t[1]*t;
293  m_t[3] = m_t[2]*t;
294  m_t[4] = 1.0/t;
295  m_t[5] = log(t);
296 
297  // iterate over the groups
298  std::vector<NasaPoly1>::const_iterator _begin, _end;
299  for (i = 0; i != m_ngroups; i++) {
300  if (t > m_tmid[i]) {
301  _begin = m_high[i].begin();
302  _end = m_high[i].end();
303  } else {
304  _begin = m_low[i].begin();
305  _end = m_low[i].end();
306  }
307  for (; _begin != _end; ++_begin) {
308  _begin->updateProperties(&m_t[0], cp_R, h_RT, s_R);
309  }
310  }
311  }
312 
313  //! Minimum temperature.
314  /*!
315  * If no argument is supplied, this
316  * method returns the minimum temperature for which \e all
317  * parameterizations are valid. If an integer index k is
318  * supplied, then the value returned is the minimum
319  * temperature for species k in the phase.
320  *
321  * @param k Species index
322  */
323  virtual doublereal minTemp(size_t k=npos) const {
324  if (k == npos) {
325  return m_tlow_max;
326  } else {
327  return m_tlow[k];
328  }
329  }
330 
331  //! Maximum temperature.
332  /*!
333  * If no argument is supplied, this
334  * method returns the maximum temperature for which \e all
335  * parameterizations are valid. If an integer index k is
336  * supplied, then the value returned is the maximum
337  * temperature for parameterization k.
338  *
339  * @param k Species index
340  */
341  virtual doublereal maxTemp(size_t k=npos) const {
342  if (k == npos) {
343  return m_thigh_min;
344  } else {
345  return m_thigh[k];
346  }
347  }
348 
349  //! The reference-state pressure for species k.
350  /*!
351  *
352  * returns the reference state pressure in Pascals for
353  * species k. If k is left out of the argument list,
354  * it returns the reference state pressure for the first
355  * species.
356  * Note that some SpeciesThermo implementations, such
357  * as those for ideal gases, require that all species
358  * in the same phase have the same reference state pressures.
359  *
360  * @param k Species index
361  */
362  virtual doublereal refPressure(size_t k=npos) const {
363  return m_p0;
364  }
365 
366  //! This utility function reports the type of parameterization
367  //! used for the species with index number index.
368  /*!
369  *
370  * @param index Species index
371  */
372  virtual int reportType(size_t index) const {
373  return NASA;
374  }
375 
376  /*!
377  * This utility function reports back the type of
378  * parameterization and all of the parameters for the
379  * species, index.
380  *
381  * @param index Species index
382  * @param type Integer type of the standard type
383  * @param c Vector of coefficients used to set the
384  * parameters for the standard state.
385  * For the NASA object, there are 15 coefficients.
386  * @param minTemp output - Minimum temperature
387  * @param maxTemp output - Maximum temperature
388  * @param refPressure output - reference pressure (Pa).
389  */
390  virtual void reportParams(size_t index, int& type,
391  doublereal* const c,
392  doublereal& minTemp,
393  doublereal& maxTemp,
394  doublereal& refPressure) const {
395  type = reportType(index);
396  if (type == NASA) {
397  size_t grp = m_group_map[index];
398  size_t pos = m_posInGroup_map[index];
399  const std::vector<NasaPoly1> &mlg = m_low[grp-1];
400  const std::vector<NasaPoly1> &mhg = m_high[grp-1];
401  const NasaPoly1* lowPoly = &(mlg[pos]);
402  const NasaPoly1* highPoly = &(mhg[pos]);
403  int itype = NASA;
404  doublereal tmid = lowPoly->maxTemp();
405  c[0] = tmid;
406  size_t n;
407  double ttemp;
408  lowPoly->reportParameters(n, itype, minTemp, ttemp, refPressure,
409  c + 1);
410  if (n != index) {
411  throw CanteraError(" ", "confused");
412  }
413  if (itype != NASA1) {
414  throw CanteraError(" ", "confused");
415  }
416  highPoly->reportParameters(n, itype, ttemp, maxTemp, refPressure,
417  c + 8);
418  if (n != index) {
419  throw CanteraError(" ", "confused");
420  }
421  if (itype != NASA1) {
422  throw CanteraError(" ", "confused");
423  }
424  } else {
425  throw CanteraError(" ", "confused");
426  }
427  }
428 
429  //! Modify parameters for the standard state
430  /*!
431  * This utility function modifies the array of coefficients.
432  * The array is the same as that returned by reportParams, so
433  * a call can first be made to reportParams to populate the
434  * array, and then modifyParams can be called to alter
435  * selected values. For the NASA object, there are 15
436  * coefficients.
437 
438  * @param index Species index
439  * @param c Vector of coefficients used to set the
440  * parameters for the standard state.
441  */
442  virtual void modifyParams(size_t index, doublereal* c) {
443  int type = reportType(index);
444  if (type == NASA) {
445  size_t grp = m_group_map[index];
446  size_t pos = m_posInGroup_map[index];
447  std::vector<NasaPoly1> &mlg = m_low[grp-1];
448  std::vector<NasaPoly1> &mhg = m_high[grp-1];
449  NasaPoly1* lowPoly = &(mlg[pos]);
450  NasaPoly1* highPoly = &(mhg[pos]);
451  doublereal tmid = lowPoly->maxTemp();
452  if (c[0] != tmid) {
453  throw CanteraError(" ", "Tmid cannot be changed");
454  }
455  lowPoly->modifyParameters(c + 1);
456  highPoly->modifyParameters(c + 8);
457  checkContinuity(m_name[index], c[0], c + 1, c + 8);
458  } else {
459  throw CanteraError(" ", "confused");
460  }
461  }
462 
463 #ifdef H298MODIFY_CAPABILITY
464  virtual doublereal reportOneHf298(const int k) const {
465 
466  int grp = m_group_map[k];
467  int pos = m_posInGroup_map[k];
468  const std::vector<NasaPoly1> &mlg = m_low[grp-1];
469  const NasaPoly1* nlow = &(mlg[pos]);
470  doublereal tmid = nlow->maxTemp();
471  double h;
472  if (298.15 <= tmid) {
473  h = nlow->reportHf298(0);
474  } else {
475  const std::vector<NasaPoly1> &mhg = m_high[grp-1];
476  const NasaPoly1* nhigh = &(mhg[pos]);
477  h = nhigh->reportHf298(0);
478  }
479  return h;
480  }
481 
482  virtual void modifyOneHf298(const int k, const doublereal Hf298New) {
483  int grp = m_group_map[k];
484  int pos = m_posInGroup_map[k];
485  std::vector<NasaPoly1> &mlg = m_low[grp-1];
486  NasaPoly1* nlow = &(mlg[pos]);
487  std::vector<NasaPoly1> &mhg = m_high[grp-1];
488  NasaPoly1* nhigh = &(mhg[pos]);
489  doublereal tmid = nlow->maxTemp();
490 
491  double hnow = reportOneHf298(k);
492  double delH = Hf298New - hnow;
493  if (298.15 <= tmid) {
494  nlow->modifyOneHf298(k, Hf298New);
495  double h = nhigh->reportHf298(0);
496  double hnew = h + delH;
497  nhigh->modifyOneHf298(k, hnew);
498  } else {
499  nhigh->modifyOneHf298(k, Hf298New);
500  double h = nlow->reportHf298(0);
501  double hnew = h + delH;
502  nlow->modifyOneHf298(k, hnew);
503  }
504 
505  }
506 #endif
507 
508 protected:
509  //! Vector of vector of NasaPoly1's for the high temp region.
510  /*!
511  * This is the high temp region representation.
512  * The first Length is equal to the number of groups.
513  * The second vector is equal to the number of species
514  * in that particular group.
515  */
516  std::vector<std::vector<NasaPoly1> > m_high;
517 
518  //! Vector of vector of NasaPoly1's for the low temp region.
519  /*!
520  * This is the low temp region representation.
521  * The first Length is equal to the number of groups.
522  * The second vector is equal to the number of species
523  * in that particular group.
524  */
525  std::vector<std::vector<NasaPoly1> > m_low;
526 
527  //! Map between the midpoint temperature, as an int, to the group number
528  /*!
529  * Length is equal to the number of groups. Only used in the setup.
530  */
531  std::map<int, int> m_index;
532 
533  //! Vector of log temperature limits
534  /*!
535  * Length is equal to the number of groups.
536  */
538 
539  //! Maximum value of the low temperature limit
540  doublereal m_tlow_max;
541 
542  //! Minimum value of the high temperature limit
543  doublereal m_thigh_min;
544 
545  //! Vector of low temperature limits (species index)
546  /*!
547  * Length is equal to number of species
548  */
550 
551  //! Vector of low temperature limits (species index)
552  /*!
553  * Length is equal to number of species
554  */
556 
557  //! Reference pressure (Pa)
558  /*!
559  * all species must have the same reference pressure.
560  */
561  doublereal m_p0;
562 
563  //! number of groups
565 
566  //! Vector of temperature polynomials
567  mutable vector_fp m_t;
568 
569  /*!
570  * This map takes as its index, the species index in the phase.
571  * It returns the group index, where the temperature polynomials
572  * for that species are stored. group indices start at 1,
573  * so a decrement is always performed to access vectors.
574  */
575  mutable std::map<size_t, size_t> m_group_map;
576 
577  /*!
578  * This map takes as its index, the species index in the phase.
579  * It returns the position index within the group, where the
580  * temperature polynomials for that species are stored.
581  */
582  mutable std::map<size_t, size_t> m_posInGroup_map;
583 
584  //! Species name as a function of the species index
585  mutable std::map<size_t, std::string> m_name;
586 
587 private:
588 
589  //! see SpeciesThermoFactory.cpp for the definition
590  /*!
591  * @param name string name of species
592  * @param tmid Mid temperature, between the two temperature regions
593  * @param clow coefficients for lower temperature region
594  * @param chigh coefficients for higher temperature region
595  */
596  void checkContinuity(std::string name, double tmid, const doublereal* clow,
597  doublereal* chigh);
598 
599  //! for internal use by checkContinuity
600  /*!
601  * @param t temperature
602  * @param c coefficient array
603  */
604  doublereal enthalpy_RT(double t, const doublereal* c) {
605  return c[0] + 0.5*c[1]*t + OneThird*c[2]*t*t
606  + 0.25*c[3]*t*t*t + 0.2*c[4]*t*t*t*t
607  + c[5]/t;
608  }
609 
610  //! for internal use by checkContinuity
611  /*!
612  * @param t temperature
613  * @param c coefficient array
614  */
615  doublereal entropy_R(double t, const doublereal* c) {
616  return c[0]*log(t) + c[1]*t + 0.5*c[2]*t*t
617  + OneThird*c[3]*t*t*t + 0.25*c[4]*t*t*t*t
618  + c[6];
619  }
620 
621 
622 
623 };
624 
625 }
626 
627 #endif
628