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