antioch-0.4.0
cea_evaluator_vec_unit.C
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 // $Id$
28 //
29 //--------------------------------------------------------------------------
30 //--------------------------------------------------------------------------
31 
32 #include "antioch_config.h"
33 
34 #include <valarray>
35 
36 #ifdef ANTIOCH_HAVE_EIGEN
37 #include "Eigen/Dense"
38 #endif
39 
40 #ifdef ANTIOCH_HAVE_METAPHYSICL
41 #include "metaphysicl/numberarray.h"
42 #endif
43 
44 #ifdef ANTIOCH_HAVE_VEXCL
45 #include "vexcl/vexcl.hpp"
46 #endif
47 
48 // Antioch
49 // Declare metaprogramming overloads before they're used
50 #include "antioch/eigen_utils_decl.h"
51 #include "antioch/metaphysicl_utils_decl.h"
52 #include "antioch/valarray_utils_decl.h"
53 #include "antioch/vexcl_utils_decl.h"
54 
55 #include "antioch/chemical_mixture.h"
56 #include "antioch/cea_evaluator.h"
57 #include "antioch/default_filename.h"
58 #include "antioch/physical_constants.h"
59 #include "antioch/cea_mixture_ascii_parsing.h"
60 
61 #include "antioch/eigen_utils.h"
62 #include "antioch/metaphysicl_utils.h"
63 #include "antioch/valarray_utils.h"
64 #include "antioch/vexcl_utils.h"
65 
66 #ifdef ANTIOCH_HAVE_GRVY
67 #include "grvy.h"
68 
69 GRVY::GRVY_Timer_Class gt;
70 #endif
71 
72 // C++
73 #include <cmath>
74 #include <limits>
75 
76 template <typename Scalar, typename TrioScalars>
77 int test_cp( const std::string& species_name, unsigned int species,
78  TrioScalars cp_exact, TrioScalars T,
79  const Antioch::CEAEvaluator<Scalar>& thermo,
80  const std::string& testname )
81 {
82  using std::abs;
83 
84  int return_flag = 0;
85 
86  const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 25;
87 
88  typedef typename Antioch::
89  template TempCache<TrioScalars> Cache;
90 
91 #ifdef ANTIOCH_HAVE_GRVY
92  gt.BeginTimer(testname);
93 #endif
94 
95  const TrioScalars cp = thermo.cp(Cache(T), species);
96 
97 #ifdef ANTIOCH_HAVE_GRVY
98  gt.EndTimer(testname);
99 #endif
100 
101  // Workaround for a non-standard gcc definition of valarray binary operator-
102  const TrioScalars diff = cp_exact - cp;
103  const TrioScalars rel_cp_error = abs(diff/cp_exact);
104 
105  if( Antioch::max(rel_cp_error) > tol )
106  {
107  std::cerr << "Error: Mismatch in species specific heat."
108  << std::setprecision
109  (std::numeric_limits<Scalar>::digits10 + 1)
110  << "\nspecies = " << species_name
111  << "\ncp = " << cp
112  << "\ncp_exact = " << cp_exact
113  << "\ndifference = " << diff
114  << "\nrelative = " << rel_cp_error
115  << "\ntolerance = " << tol
116  << "\nT = " << T << std::endl;
117  return_flag = 1;
118  }
119 
120  return return_flag;
121 }
122 
123 
124 template <typename Scalar>
125 Scalar cp( Scalar T, Scalar a0, Scalar a1, Scalar a2,
126  Scalar a3, Scalar a4, Scalar a5, Scalar a6 )
127 {
128  if( T < 200.1)
129  T = 200.1;
130 
131  return a0/(T*T) + a1/T + a2 + a3*T + a4*(T*T) + a5*(T*T*T) + a6*(T*T*T*T);
132 }
133 
134 
135 template <typename TrioScalars>
136 int vectester(const TrioScalars& example, const std::string& testname)
137 {
138  typedef typename Antioch::value_type<TrioScalars>::type Scalar;
139 
140  std::vector<std::string> species_str_list;
141  const unsigned int n_species = 5;
142  species_str_list.reserve(n_species);
143  species_str_list.push_back( "N2" );
144  species_str_list.push_back( "O2" );
145  species_str_list.push_back( "N" );
146  species_str_list.push_back( "O" );
147  species_str_list.push_back( "NO" );
148 
149  const Scalar Mm_N = 14.008e-3L;
150  const Scalar Mm_O = 16.000e-3L;
151  const Scalar Mm_N2 = 2.L * Mm_N;
152  const Scalar Mm_O2 = 2.L * Mm_O;
153  const Scalar Mm_NO = Mm_N + Mm_O;
154 
155  Antioch::ChemicalMixture<Scalar> chem_mixture( species_str_list );
156 
157  Antioch::CEAThermoMixture<Scalar> cea_mixture( chem_mixture );
158  Antioch::read_cea_mixture_data_ascii( cea_mixture, Antioch::DefaultFilename::thermo_data() );
159  Antioch::CEAEvaluator<Scalar> thermo( cea_mixture );
160 
161  //const Scalar P = 100000.0;
162  TrioScalars T = example;
163  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
164  {
165  T[3*tuple] = 190.0;
166  T[3*tuple+1] = 1500.0;
167  T[3*tuple+2] = 10000.0;
168  }
169 
170  const Scalar R_N2 = Antioch::Constants::R_universal<Scalar>()/Mm_N2;
171  const Scalar R_O2 = Antioch::Constants::R_universal<Scalar>()/Mm_O2;
172  const Scalar R_N = Antioch::Constants::R_universal<Scalar>()/Mm_N;
173  const Scalar R_O = Antioch::Constants::R_universal<Scalar>()/Mm_O;
174  const Scalar R_NO = Antioch::Constants::R_universal<Scalar>()/Mm_NO;
175 
176  int return_flag = 0;
177 
178  // Test N2 cp
179  {
180  unsigned int index = 0;
181  TrioScalars cp_N2 = example;
182  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
183  {
184  cp_N2[3*tuple ] = R_N2*cp( Scalar(T[0]), Scalar(2.21037122e+04), Scalar(-3.81846145e+02), Scalar(6.08273815e+00),
185  Scalar(-8.53091381e-03), Scalar(1.38464610e-05), Scalar(-9.62579293e-09), Scalar(2.51970560e-12));
186 
187  cp_N2[3*tuple+1] = R_N2*cp( Scalar(T[1]), Scalar(5.87709908e+05), Scalar(-2.23924255e+03), Scalar(6.06694267e+00),
188  Scalar(-6.13965296e-04), Scalar(1.49179819e-07), Scalar(-1.92309442e-11), Scalar(1.06194871e-15) );
189 
190  cp_N2[3*tuple+2] = R_N2*cp( Scalar(T[2]), Scalar(8.30971200e+08), Scalar(-6.42048187e+05), Scalar(2.02020507e+02), Scalar(-3.06501961e-02),
191  Scalar(2.48685558e-06), Scalar(-9.70579208e-11), Scalar(1.43751673e-15));
192  }
193 
194  const Antioch::Species species = chem_mixture.species_list()[index];
195  const std::string species_name = chem_mixture.species_inverse_name_map().find(species)->second;
196 
197  int return_flag_temp = 0;
198  return_flag_temp = test_cp( species_name, index, cp_N2, T, thermo, testname );
199  if( return_flag_temp != 0 ) return_flag = 1;
200  }
201 
202  // Test O2 cp
203  {
204  unsigned int index = 1;
205  TrioScalars cp_O2 = example;
206  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
207  {
208  cp_O2[3*tuple ] = R_O2*cp( Scalar(T[0]), Scalar(-3.42556269e+04), Scalar(4.84699986e+02), Scalar(1.11901159e+00),
209  Scalar(4.29388743e-03), Scalar(-6.83627313e-07), Scalar(-2.02337478e-09) , Scalar(1.03904064e-12) );
210 
211  cp_O2[3*tuple+1] = R_O2*cp( Scalar(T[1]), Scalar(-1.03793994e+06), Scalar(2.34483275e+03), Scalar(1.81972949e+00), Scalar(1.26784887e-03),
212  Scalar(-2.18807142e-07), Scalar(2.05372411e-11), Scalar(-8.19349062e-16) );
213 
214  cp_O2[3*tuple+2] = R_O2*cp( Scalar(T[2]), Scalar(4.97515261e+08), Scalar(-2.86602339e+05), Scalar(6.69015464e+01), Scalar(-6.16971869e-03),
215  Scalar(3.01623757e-07), Scalar(-7.42087888e-12), Scalar(7.27744063e-17));
216  }
217 
218  const Antioch::Species species = chem_mixture.species_list()[index];
219  const std::string species_name = chem_mixture.species_inverse_name_map().find(species)->second;
220 
221  int return_flag_temp = 0;
222  return_flag_temp = test_cp( species_name, index, cp_O2, T, thermo, testname );
223  if( return_flag_temp != 0 ) return_flag = 1;
224  }
225 
226  // Test N cp
227  {
228  unsigned int index = 2;
229  TrioScalars cp_N = example;
230  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
231  {
232  cp_N[3*tuple ] = R_N*cp( Scalar(T[0]), Scalar(0.00000000e+00), Scalar(0.00000000e+00), Scalar(2.50000000e+00), Scalar(0.00000000e+00),
233  Scalar(0.00000000e+00), Scalar(0.00000000e+00), Scalar(0.00000000e+00));
234 
235  cp_N[3*tuple+1] = R_N*cp( Scalar(T[1]), Scalar(8.87650138e+04), Scalar(-1.07123150e+02), Scalar(2.36218829e+00), Scalar(2.91672008e-04),
236  Scalar(-1.72951510e-07), Scalar(4.01265788e-11), Scalar(-2.67722757e-15) );
237 
238  cp_N[3*tuple+2] = R_N*cp( Scalar(T[2]), Scalar(5.47518105e+08), Scalar(-3.10757498e+05), Scalar(6.91678274e+01), Scalar(-6.84798813e-03),
239  Scalar(3.82757240e-07), Scalar(-1.09836771e-11), Scalar(1.27798602e-16));
240  }
241 
242  const Antioch::Species species = chem_mixture.species_list()[index];
243  const std::string species_name = chem_mixture.species_inverse_name_map().find(species)->second;
244 
245  int return_flag_temp = 0;
246  return_flag_temp = test_cp( species_name, index, cp_N, T, thermo, testname );
247  if( return_flag_temp != 0 ) return_flag = 1;
248  }
249 
250 
251  // Test O cp
252  {
253  unsigned int index = 3;
254  TrioScalars cp_O = example;
255  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
256  {
257  cp_O[3*tuple ] = R_O*cp( Scalar(T[0]), Scalar(-7.95361130e+03) , Scalar(1.60717779e+02), Scalar(1.96622644e+00), Scalar(1.01367031e-03),
258  Scalar(-1.11041542e-06), Scalar(6.51750750e-10), Scalar(-1.58477925e-13) );
259 
260  cp_O[3*tuple+1] = R_O*cp( Scalar(T[1]), Scalar(2.61902026e+05), Scalar(-7.29872203e+02), Scalar(3.31717727e+00), Scalar(-4.28133436e-04),
261  Scalar(1.03610459e-07), Scalar(-9.43830433e-12), Scalar(2.72503830e-16) );
262 
263  cp_O[3*tuple+2] = R_O*cp( Scalar(T[2]), Scalar(1.77900426e+08), Scalar(-1.08232826e+05), Scalar(2.81077837e+01), Scalar(-2.97523226e-03),
264  Scalar(1.85499753e-07), Scalar(-5.79623154e-12), Scalar(7.19172016e-17) );
265  }
266 
267  const Antioch::Species species = chem_mixture.species_list()[index];
268  const std::string species_name = chem_mixture.species_inverse_name_map().find(species)->second;
269 
270  int return_flag_temp = 0;
271  return_flag_temp = test_cp( species_name, index, cp_O, T, thermo, testname );
272  if( return_flag_temp != 0 ) return_flag = 1;
273  }
274 
275 
276  // Test NO cp
277  {
278  unsigned int index = 4;
279  TrioScalars cp_NO = example;
280  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
281  {
282  cp_NO[3*tuple ] = R_NO*cp( Scalar(T[0]), Scalar(-1.14391658e+04), Scalar(1.53646774e+02), Scalar(3.43146865e+00), Scalar(-2.66859213e-03),
283  Scalar(8.48139877e-06), Scalar(-7.68511079e-09), Scalar(2.38679758e-12) );
284 
285  cp_NO[3*tuple+1] = R_NO*cp( Scalar(T[1]), Scalar(2.23903708e+05), Scalar(-1.28965624e+03), Scalar(5.43394039e+00), Scalar(-3.65605546e-04),
286  Scalar(9.88101763e-08), Scalar(-1.41608327e-11), Scalar(9.38021642e-16) );
287 
288  cp_NO[3*tuple+2] = R_NO*cp( Scalar(T[2]), Scalar(-9.57530764e+08), Scalar(5.91243671e+05), Scalar(-1.38456733e+02), Scalar(1.69433998e-02),
289  Scalar(-1.00735146e-06), Scalar(2.91258526e-11), Scalar(-3.29511091e-16) );
290  }
291 
292  const Antioch::Species species = chem_mixture.species_list()[index];
293  const std::string species_name = chem_mixture.species_inverse_name_map().find(species)->second;
294 
295  int return_flag_temp = 0;
296  return_flag_temp = test_cp( species_name, index, cp_NO, T, thermo, testname );
297  if( return_flag_temp != 0 ) return_flag = 1;
298  }
299 
300  return return_flag;
301 }
302 
303 
304 int main()
305 {
306  int returnval = 0;
307 
308  returnval = returnval ||
309  vectester (std::valarray<float>(3*ANTIOCH_N_TUPLES), "valarray<float>");
310  returnval = returnval ||
311  vectester (std::valarray<double>(3*ANTIOCH_N_TUPLES), "valarray<double>");
312 // We're not getting the full long double precision yet?
313 // returnval = returnval ||
314 // vectester<long double, std::valarray<long double> >
315 // (std::valarray<long double>(3*ANTIOCH_N_TUPLES), "valarray<ld>");
316 #ifdef ANTIOCH_HAVE_EIGEN
317  returnval = returnval ||
318  vectester (Eigen::Array<float, 3*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXf");
319  returnval = returnval ||
320  vectester (Eigen::Array<double, 3*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXd");
321 // returnval = returnval ||
322 // vectester (Eigen::Array<long double, 3*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXld");
323 #endif
324 #ifdef ANTIOCH_HAVE_METAPHYSICL
325  returnval = returnval ||
326  vectester (MetaPhysicL::NumberArray<3*ANTIOCH_N_TUPLES, float> (0), "NumberArray<float>");
327  returnval = returnval ||
328  vectester (MetaPhysicL::NumberArray<3*ANTIOCH_N_TUPLES, double> (0), "NumberArray<double>");
329 // returnval = returnval ||
330 // vectester (MetaPhysicL::NumberArray<3*ANTIOCH_N_TUPLES, long double> (0)), "NumberArray<ld>");
331 #endif
332 #ifdef ANTIOCH_HAVE_VEXCL
333  vex::Context ctx_f (vex::Filter::All);
334  if (!ctx_f.empty())
335  returnval = returnval ||
336  vectester (vex::vector<float> (ctx_f, 3*ANTIOCH_N_TUPLES), "vex::vector<float>");
337 
338  vex::Context ctx_d (vex::Filter::DoublePrecision);
339  if (!ctx_d.empty())
340  returnval = returnval ||
341  vectester (vex::vector<double> (ctx_d, 3*ANTIOCH_N_TUPLES), "vex::vector<double>");
342 #endif
343 
344 #ifdef ANTIOCH_HAVE_GRVY
345  gt.Finalize();
346  gt.Summarize();
347 #endif
348 
349  return returnval;
350 }
Antioch::max
Antioch::enable_if_c< is_eigen< T >::value, typename value_type< T >::type >::type max(const T &in)
Definition: eigen_utils.h:88
Antioch::Species
unsigned int Species
Definition: chemical_mixture.h:49
Antioch::ChemicalMixture::species_list
const std::vector< Species > & species_list() const
Definition: chemical_mixture.h:245
vectester
int vectester(const TrioScalars &example, const std::string &testname)
Definition: cea_evaluator_vec_unit.C:136
main
int main()
Definition: cea_evaluator_vec_unit.C:304
cp
Scalar cp(Scalar T, Scalar a0, Scalar a1, Scalar a2, Scalar a3, Scalar a4, Scalar a5, Scalar a6)
Definition: cea_evaluator_vec_unit.C:125
Antioch::ChemicalMixture::species_inverse_name_map
const std::map< Species, std::string > & species_inverse_name_map() const
Definition: chemical_mixture.h:266
Antioch::DefaultSourceFilename::thermo_data
static const std::string & thermo_data()
Definition: default_filename.h:76
Antioch::ChemicalMixture
Class storing chemical mixture properties.
Definition: chemical_mixture.h:75
Antioch::read_cea_mixture_data_ascii
void read_cea_mixture_data_ascii(CEAThermoMixture< NumericType > &thermo, const std::string &filename)
Definition: cea_mixture_ascii_parsing.C:45
test_cp
int test_cp(const std::string &species_name, unsigned int species, TrioScalars cp_exact, TrioScalars T, const Antioch::CEAEvaluator< Scalar > &thermo, const std::string &testname)
Definition: cea_evaluator_vec_unit.C:77
Antioch::NASAEvaluator< CoeffType, CEACurveFit< CoeffType > >::cp
StateType cp(const TempCache< StateType > &cache, unsigned int species) const
We currently need different specializations for scalar vs vector inputs here.
Generated on Thu Jul 7 2016 11:09:46 for antioch-0.4.0 by  doxygen 1.8.8