antioch-0.4.0
cea_thermo.h
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // Antioch - A Gas Dynamics Thermochemistry Library
5 //
6 // Copyright (C) 2014-2016 Paul T. Bauman, Benjamin S. Kirk,
7 // Sylvain Plessis, Roy H. Stonger
8 //
9 // Copyright (C) 2013 The PECOS Development Team
10 //
11 // This library is free software; you can redistribute it and/or
12 // modify it under the terms of the Version 2.1 GNU Lesser General
13 // Public License as published by the Free Software Foundation.
14 //
15 // This library is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
23 // Boston, MA 02110-1301 USA
24 //
25 //-----------------------------------------------------------------------el-
26 
27 
28 #ifndef ANTIOCH_CEA_THERMO_H
29 #define ANTIOCH_CEA_THERMO_H
30 
31 // Antioch
32 #include "antioch_config.h"
33 #include "antioch/cea_curve_fit.h"
38 #include "antioch/input_utils.h"
40 
41 // C++
42 #include <iomanip>
43 #include <vector>
44 #include <cmath>
45 
46 namespace Antioch
47 {
48  // Forward declarations
49  template<typename CoeffType>
50  class NASA9CurveFit;
51 
52  template<typename CoeffType=double>
53  class CEAThermodynamics
54  {
55  public:
56 
57  CEAThermodynamics( const ChemicalMixture<CoeffType>& chem_mixture,
58  const std::string & filename = DefaultFilename::thermo_data() );
59 
61 
62  virtual ~CEAThermodynamics();
63 
64  template<typename StateType=CoeffType>
65  class Cache
66  {
67  public:
68  const StateType &T;
69  StateType T2;
70  StateType T3;
71  StateType T4;
72  StateType lnT;
73 
74  explicit Cache(const StateType &T_in)
75  : T(T_in), T2(T*T), T3(T2*T), T4(T2*T2), lnT(T_in) {
76  using std::log;
77  lnT = log(T);
78  }
79 
80  Cache(const StateType &T_in,
81  const StateType &T2_in,
82  const StateType &T3_in,
83  const StateType &T4_in,
84  const StateType &lnT_in)
85  : T(T_in), T2(T2_in), T3(T3_in), T4(T4_in), lnT(lnT_in) {
87  }
88  private:
89  Cache();
90  };
91 
92  void add_curve_fit( const std::string& species_name, const std::vector<CoeffType>& coeffs );
93 
94  // ! compatibility for parsers
95  void add_curve_fit( const std::string& species_name, const std::vector<CoeffType>& coeffs, const std::vector<CoeffType> & /*temps */ )
96  {this->add_curve_fit(species_name,coeffs);}
97 
99  bool check() const;
100 
102  template<typename StateType>
103  StateType
104  cp( const Cache<StateType> &cache, unsigned int species ) const;
105 
106  template<typename StateType, typename VectorStateType>
107  typename enable_if_c<
109  >::type
110  cp( const Cache<StateType> &cache, const VectorStateType& mass_fractions ) const;
111 
112  template<typename StateType>
113  StateType cv( const Cache<StateType> &cache, unsigned int species ) const;
114 
115  template<typename StateType, typename VectorStateType>
116  typename enable_if_c<
117  has_size<VectorStateType>::value, StateType
118  >::type
119  cv( const Cache<StateType> &cache, const VectorStateType& mass_fractions ) const;
120 
121  template<typename StateType>
122  StateType h( const Cache<StateType> &cache, unsigned int species ) const;
123 
124  template<typename StateType, typename VectorStateType>
125  typename enable_if_c<
126  has_size<VectorStateType>::value, void
127  >::type
128  h( const Cache<StateType> &cache, VectorStateType& h ) const;
129 
131  template<typename StateType>
132  StateType
133  h_RT_minus_s_R( const Cache<StateType> &cache, unsigned int species ) const;
134 
135  template<typename StateType, typename VectorStateType>
136  typename enable_if_c<
137  has_size<VectorStateType>::value, void
138  >::type
139  h_RT_minus_s_R( const Cache<StateType> &cache, VectorStateType& h_RT_minus_s_R ) const;
140 
141 
142 
144  template<typename StateType>
145  StateType
146  dh_RT_minus_s_R_dT( const Cache<StateType> &cache, unsigned int species ) const;
147 
148  template<typename StateType, typename VectorStateType>
149  typename enable_if_c<
150  has_size<VectorStateType>::value, void
151  >::type
152  dh_RT_minus_s_R_dT( const Cache<StateType> &cache, VectorStateType& dh_RT_minus_s_R_dT ) const;
153 
154 
155  template<typename StateType>
156  StateType cp_over_R( const Cache<StateType> &cache, unsigned int species ) const;
157 
158  template<typename StateType>
159  StateType h_over_RT( const Cache<StateType> &cache, unsigned int species ) const;
160 
161  template<typename StateType>
162  StateType s_over_R( const Cache<StateType> &cache, unsigned int species ) const;
163 
165 
166  protected:
167 
169 
170  std::vector<CEACurveFit<CoeffType>* > _species_curve_fits;
171 
172  std::vector<CoeffType> _cp_at_200p1;
173 
174  private:
175 
177 
179  };
180 
181 
182  /* ------------------------- Inline Functions -------------------------*/
183 
184  template<typename CoeffType>
185  inline
186  CEAThermodynamics<CoeffType>::CEAThermodynamics( const ChemicalMixture<CoeffType>& chem_mixture, const std::string & filename )
187  : _chem_mixture(chem_mixture),
188  _species_curve_fits(chem_mixture.n_species(), NULL)
189  {
190 
191  // should not use this object
193 
194  // Read in CEA coefficients. Note this assumes chem_mixture is fully constructed.
196  read_cea_mixture_data_ascii(*this, filename);
197 
198  // Cache cp values at small temperatures
199  _cp_at_200p1.reserve( _species_curve_fits.size() );
200  for( unsigned int s = 0; s < _species_curve_fits.size(); s++ )
201  {
202  _cp_at_200p1.push_back(CoeffType(0));
203  _cp_at_200p1[s] = this->cp( Cache<CoeffType>(200.1), s );
204  }
205 
206  return;
207  }
208 
209 
210  template<typename CoeffType>
211  inline
213  {
214  // Clean up all the CEACurveFits we created
215  for( typename std::vector<CEACurveFit<CoeffType>* >::iterator it = _species_curve_fits.begin();
216  it < _species_curve_fits.end(); ++it )
217  {
218  delete (*it);
219  }
220 
221  return;
222  }
223 
224 
225  template<typename CoeffType>
226  inline
227  void CEAThermodynamics<CoeffType>::add_curve_fit( const std::string& species_name,
228  const std::vector<CoeffType>& coeffs )
229  {
230  antioch_assert( _chem_mixture.species_name_map().find(species_name) !=
231  _chem_mixture.species_name_map().end() );
232 
233  unsigned int s = _chem_mixture.species_name_map().find(species_name)->second;
234 
235  antioch_assert_less_equal( s, _species_curve_fits.size() );
236  antioch_assert( !_species_curve_fits[s] );
237 
238  _species_curve_fits[s] = new CEACurveFit<CoeffType>(coeffs);
239  return;
240  }
241 
242 
243  template<typename CoeffType>
244  inline
246  {
247  bool valid = true;
248 
249  for( typename std::vector<CEACurveFit<CoeffType>* >::const_iterator it = _species_curve_fits.begin();
250  it != _species_curve_fits.end(); ++ it )
251  {
252  if( !(*it) )
253  valid = false;
254  }
255 
256  return valid;
257  }
258 
259 
260  template<typename CoeffType>
261  template<typename StateType>
262  inline
263  StateType
264  CEAThermodynamics<CoeffType>::cp( const Cache<StateType> &cache, unsigned int species ) const
265  {
266  typedef typename Antioch::value_type<StateType>::type ScalarType;
267 
268  // T < 200.1 ? cp_at_200p1 : R * cp_over_R
269  return
271  (cache.T < ScalarType(200.1),
273  (cache.T,_cp_at_200p1[species]),
274  StateType
275  (this->_chem_mixture.R(species) *
276  this->cp_over_R(cache, species)));
277  }
278 
279 
280  template<typename CoeffType>
281  template<typename StateType, typename VectorStateType>
282  inline
283  typename enable_if_c<
285  >::type
287  const VectorStateType& mass_fractions ) const
288  {
289  antioch_assert_equal_to( mass_fractions.size(), _species_curve_fits.size() );
290  antioch_assert_greater( mass_fractions.size(), 0 );
291 
292  StateType cp = mass_fractions[0]*this->cp(cache,0);
293 
294  for( unsigned int s = 1; s < _species_curve_fits.size(); s++ )
295  {
296  cp += mass_fractions[s]*this->cp(cache,s);
297  }
298 
299  return cp;
300  }
301 
302 
303  template<typename CoeffType>
304  template<typename StateType>
305  inline
306  StateType CEAThermodynamics<CoeffType>::h( const Cache<StateType> &cache, unsigned int species ) const
307  {
308  return this->_chem_mixture.R(species)*cache.T*this->h_over_RT(cache,species);
309  }
310 
311 
312  template<typename CoeffType>
313  template<typename StateType, typename VectorStateType>
314  inline
315  typename enable_if_c<
316  has_size<VectorStateType>::value, void
317  >::type
318  CEAThermodynamics<CoeffType>::h( const Cache<StateType> &cache, VectorStateType& h ) const
319  {
320  antioch_assert_equal_to( h.size(), _species_curve_fits.size() );
321 
322  for( unsigned int s = 0; s < _species_curve_fits.size(); s++ )
323  {
324  h[s] = this->_chem_mixture.R(s)*cache.T*this->h_over_RT(cache,s);
325  }
326 
327  return;
328  }
329 
330 
331  template<typename CoeffType>
332  template<typename StateType>
333  inline
334  StateType
335  CEAThermodynamics<CoeffType>::cp_over_R( const Cache<StateType> &cache, unsigned int species ) const
336  {
337  antioch_assert_less( species, _species_curve_fits.size() );
338  // FIXME - we need assert_less to be vectorizable
339  // antioch_assert_less( _species_curve_fits[species]->interval(cache.T),
340  // _species_curve_fits[species]->n_intervals() );
341 
342  typedef typename
344  const UIntType interval = this->_species_curve_fits[species]->interval(cache.T);
345 
346  const unsigned int begin_interval = Antioch::min(interval);
347  const unsigned int end_interval = Antioch::max(interval)+1;
348 
349  // FIXME - this needs expression templates to be faster...
350 
351  StateType returnval = Antioch::zero_clone(cache.T);
352 
353  /* cp/R = a0*T^-2 + a1*T^-1 + a2 + a3*T + a4*T^2 + a5*T^3 + a6*T^4 */
354  for (unsigned int i=begin_interval; i != end_interval; ++i)
355  {
356  const CoeffType * const a =
357  this->_species_curve_fits[species]->coefficients(i);
358  returnval = Antioch::if_else
359  (interval == i,
360  StateType(a[0]/cache.T2 + a[1]/cache.T + a[2] + a[3]*cache.T +
361  a[4]*cache.T2 + a[5]*cache.T3 + a[6]*cache.T4),
362  returnval);
363  }
364 
365  return returnval;
366  }
367 
368 
369  template<typename CoeffType>
370  template<typename StateType>
371  inline
372  StateType CEAThermodynamics<CoeffType>::h_over_RT( const Cache<StateType> &cache, unsigned int species ) const
373  {
374  antioch_assert_less( species, _species_curve_fits.size() );
375  antioch_assert_less( _species_curve_fits[species]->interval(cache.T),
376  _species_curve_fits[species]->n_intervals() );
377 
378  const unsigned int interval = this->_species_curve_fits[species]->interval(cache.T);
379 
380  const CoeffType *a = this->_species_curve_fits[species]->coefficients(interval);
381 
382  typedef typename
384 
385  /* h/RT = -a0*T^-2 + a1*T^-1*lnT + a2 + a3*T/2 + a4*T^2/3 + a5*T^3/4 + a6*T^4/5 + a7/T */
386  return -a[0]/cache.T2 + a[1]*cache.lnT/cache.T + a[2] +
387  a[3]*cache.T/raw_type(2) + a[4]*cache.T2/raw_type(3) +
388  a[5]*cache.T3/raw_type(4) + a[6]*cache.T4/raw_type(5) +
389  a[7]/cache.T;
390  }
391 
392 
393  template<typename CoeffType>
394  template<typename StateType>
395  inline
396  StateType CEAThermodynamics<CoeffType>::s_over_R( const Cache<StateType> &cache, unsigned int species ) const
397  {
398  antioch_assert_less( species, _species_curve_fits.size() );
399  antioch_assert_less( _species_curve_fits[species]->interval(cache.T),
400  _species_curve_fits[species]->n_intervals() );
401 
402  const unsigned int interval = this->_species_curve_fits[species]->interval(cache.T);
403 
404  const CoeffType *a = this->_species_curve_fits[species]->coefficients(interval);
405 
406  typedef typename
408 
409  /* s/R = -a0*T^-2/2 - a1*T^-1 + a2*lnT + a3*T + a4*T^2/2 + a5*T^3/3 + a6*T^4/4 + a8 */
410  return -a[0]/cache.T2/raw_type(2) - a[1]/cache.T +
411  a[2]*cache.lnT + a[3]*cache.T + a[4]*cache.T2/raw_type(2) +
412  a[5]*cache.T3/raw_type(3) + a[6]*cache.T4/raw_type(4) + a[8];
413  }
414 
415 
416  template<typename CoeffType>
417  template<typename StateType>
418  inline
419  StateType
420  CEAThermodynamics<CoeffType>::h_RT_minus_s_R( const Cache<StateType> &cache, unsigned int species ) const
421  {
422  antioch_assert_less( species, _species_curve_fits.size() );
423  // FIXME - we need assert_less to be vectorizable
424 // antioch_assert_less( _species_curve_fits[species]->interval(cache.T),
425 // _species_curve_fits[species]->n_intervals() );
426 
427  typedef typename
429  const UIntType interval = this->_species_curve_fits[species]->interval(cache.T);
430  const unsigned int begin_interval = Antioch::min(interval);
431  const unsigned int end_interval = Antioch::max(interval)+1;
432 
433  // FIXME - this needs expression templates to be faster...
434 
435  StateType returnval = Antioch::zero_clone(cache.T);
436 
437  typedef typename
439 
440  /* h/RT = -a[0]/T2 + a[1]*lnT/T + a[2] + a[3]*T/2. + a[4]*T2/3. + a[5]*T3/4. + a[6]*T4/5. + a[7]/T,
441  s/R = -a[0]/T2/2. - a[1]/T + a[2]*lnT + a[3]*T + a[4]*T2/2. + a[5]*T3/3. + a[6]*T4/4. + a[8] */
442  for (unsigned int i=begin_interval; i != end_interval; ++i)
443  {
444  const CoeffType * const a =
445  this->_species_curve_fits[species]->coefficients(i);
446  returnval = Antioch::if_else
447  (interval == i,
448  StateType(-a[0]/cache.T2/raw_type(2) +
449  (a[1] + a[7])/cache.T +
450  a[1]*cache.lnT/cache.T - a[2]*cache.lnT +
451  (a[2] - a[8]) - a[3]*cache.T/raw_type(2) -
452  a[4]*cache.T2/raw_type(6) -
453  a[5]*cache.T3/raw_type(12) -
454  a[6]*cache.T4/raw_type(20)),
455  returnval);
456  }
457 
458  return returnval;
459  }
460 
461 
462  template<typename CoeffType>
463  template<typename StateType, typename VectorStateType>
464  inline
465  typename enable_if_c<
466  has_size<VectorStateType>::value, void
467  >::type
469  VectorStateType& h_RT_minus_s_R ) const
470  {
471  antioch_assert_equal_to( h_RT_minus_s_R.size(), _species_curve_fits.size() );
472 
473  for( unsigned int s = 0; s < _species_curve_fits.size(); s++ )
474  {
475  h_RT_minus_s_R[s] = this->h_RT_minus_s_R(cache,s);
476  }
477 
478  return;
479  }
480 
481 
482 
483  template<typename CoeffType>
484  template<typename StateType>
485  inline
486  StateType
487  CEAThermodynamics<CoeffType>::dh_RT_minus_s_R_dT( const Cache<StateType> &cache, unsigned int species ) const
488  {
489  antioch_assert_less( species, _species_curve_fits.size() );
490  // FIXME - we need assert_less to be vectorizable
491  // antioch_assert_less( _species_curve_fits[species]->interval(cache.T),
492  // _species_curve_fits[species]->n_intervals() );
493 
494  typedef typename
496  const UIntType interval = this->_species_curve_fits[species]->interval(cache.T);
497 
498  const unsigned int begin_interval = Antioch::min(interval);
499  const unsigned int end_interval = Antioch::max(interval)+1;
500 
501  // FIXME - this needs expression templates to be faster...
502 
503  StateType returnval = Antioch::zero_clone(cache.T);
504 
505  typedef typename
507 
508  /* h/RT = -a[0]/T2 + a[1]*lnT/T + a[2] + a[3]*T/2. + a[4]*T2/3. + a[5]*T3/4. + a[6]*T4/5. + a[7]/T,
509  s/R = -a[0]/T2/2. - a[1]/T + a[2]*lnT + a[3]*T + a[4]*T2/2. + a[5]*T3/3. + a[6]*T4/4. + a[8] */
510  for (unsigned int i=begin_interval; i != end_interval; ++i)
511  {
512  const CoeffType * const a =
513  this->_species_curve_fits[species]->coefficients(i);
514  returnval = Antioch::if_else
515  (interval == i,
516  StateType(a[0]/cache.T3 - a[7]/cache.T2 -
517  a[1]*cache.lnT/cache.T2 - a[2]/cache.T -
518  a[3]/raw_type(2) - a[4]*cache.T/raw_type(3) -
519  a[5]*cache.T2/raw_type(4) -
520  a[6]*cache.T3/raw_type(5)),
521  returnval);
522  }
523 
524  return returnval;
525  }
526 
527 
528  template<typename CoeffType>
529  template<typename StateType, typename VectorStateType>
530  inline
531  typename enable_if_c<
532  has_size<VectorStateType>::value, void
533  >::type
535  VectorStateType& dh_RT_minus_s_R_dT ) const
536  {
537  antioch_assert_equal_to( dh_RT_minus_s_R_dT.size(), _species_curve_fits.size() );
538 
539  for( unsigned int s = 0; s < _species_curve_fits.size(); s++ )
540  {
541  dh_RT_minus_s_R_dT[s] = this->dh_RT_minus_s_R_dT(cache,s);
542  }
543 
544  return;
545  }
546 
547 
548 
549 
550  template<typename CoeffType>
551  template<typename StateType>
552  inline
553  StateType CEAThermodynamics<CoeffType>::cv( const Cache<StateType> &cache, unsigned int species ) const
554  {
555  return this->cp(cache,species) - _chem_mixture.R(species);
556  }
557 
558  template<typename CoeffType>
559  template<typename StateType, typename VectorStateType>
560  inline
561  typename enable_if_c<
562  has_size<VectorStateType>::value, StateType
563  >::type
565  const VectorStateType& mass_fractions ) const
566  {
567  return this->cp(cache,mass_fractions) - _chem_mixture.R(mass_fractions);
568  }
569 
570  template<typename CoeffType>
571  inline
573  {
574  return _chem_mixture;
575  }
576 
577 } // end namespace Antioch
578 
579 #endif //ANTIOCH_CEA_THERMO_H
Cache(const StateType &T_in, const StateType &T2_in, const StateType &T3_in, const StateType &T4_in, const StateType &lnT_in)
Definition: cea_thermo.h:80
Antioch::enable_if_c< is_eigen< T >::value, typename value_type< T >::type >::type max(const T &in)
Definition: eigen_utils.h:88
#define antioch_assert(asserted)
std::vector< CEACurveFit< CoeffType > * > _species_curve_fits
Definition: cea_thermo.h:170
#define antioch_assert_equal_to(expr1, expr2)
#define antioch_assert_greater(expr1, expr2)
const ChemicalMixture< CoeffType > & _chem_mixture
Definition: cea_thermo.h:168
StateType cp_over_R(const Cache< StateType > &cache, unsigned int species) const
Definition: cea_thermo.h:335
virtual ~CEAThermodynamics()
Destructor.
Definition: cea_thermo.h:212
#define antioch_assert_less(expr1, expr2)
#define antioch_assert_less_equal(expr1, expr2)
Scalar cp(Scalar T, Scalar a0, Scalar a1, Scalar a2, Scalar a3, Scalar a4, Scalar a5, Scalar a6)
Antioch::enable_if_c< is_eigen< T >::value, typename value_type< T >::type >::type min(const T &in)
Definition: eigen_utils.h:98
StateType cv(const Cache< StateType > &cache, unsigned int species) const
Definition: cea_thermo.h:553
_Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > constant_clone(const _Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > &ex, const Scalar &value)
Definition: eigen_utils.h:181
Cache(const StateType &T_in)
Definition: cea_thermo.h:74
enable_if_c< is_eigen< T1 >::value &&is_eigen< T2 >::value, typename state_type< T1 >::type >::type if_else(const Condition &condition, const T1 &if_true, const T2 &if_false)
Definition: eigen_utils.h:250
StateType dh_RT_minus_s_R_dT(const Cache< StateType > &cache, unsigned int species) const
We currently need different specializations for scalar vs vector inputs here.
Definition: cea_thermo.h:487
This class only differs from NASA9CurveFit in the construction.
Definition: ascii_parser.h:72
void add_curve_fit(const std::string &species_name, const std::vector< CoeffType > &coeffs, const std::vector< CoeffType > &)
Definition: cea_thermo.h:95
static const std::string & thermo_data()
const ChemicalMixture< CoeffType > & chemical_mixture() const
Definition: cea_thermo.h:572
#define antioch_deprecated()
StateType h(const Cache< StateType > &cache, unsigned int species) const
Definition: cea_thermo.h:306
std::vector< CoeffType > _cp_at_200p1
Definition: cea_thermo.h:172
CEAThermodynamics()
Default constructor.
void add_curve_fit(const std::string &species_name, const std::vector< CoeffType > &coeffs)
Definition: cea_thermo.h:227
Class storing chemical mixture properties.
The parameters are reduced parameters.
StateType h_RT_minus_s_R(const Cache< StateType > &cache, unsigned int species) const
We currently need different specializations for scalar vs vector inputs here.
Definition: cea_thermo.h:420
void read_cea_mixture_data_ascii(CEAThermoMixture< NumericType > &thermo, const std::string &filename)
bool check() const
Checks that curve fits have been specified for all species in the mixture.
Definition: cea_thermo.h:245
_Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > zero_clone(const _Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > &ex)
Definition: eigen_utils.h:145
StateType h_over_RT(const Cache< StateType > &cache, unsigned int species) const
Definition: cea_thermo.h:372
StateType s_over_R(const Cache< StateType > &cache, unsigned int species) const
Definition: cea_thermo.h:396
StateType cp(const Cache< StateType > &cache, unsigned int species) const
We currently need different specializations for scalar vs vector inputs here.
Definition: cea_thermo.h:264

Generated on Thu Jul 7 2016 11:09:46 for antioch-0.4.0 by  doxygen 1.8.8