antioch-0.4.0
hercourtessen_rate_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 // valarray has to be declared before Antioch or gcc can't find the
28 // right versions of exp() and pow() to use??
29 
30 #include "antioch_config.h"
31 
32 #include <valarray>
33 
34 #ifdef ANTIOCH_HAVE_EIGEN
35 #include "Eigen/Dense"
36 #endif
37 
38 #ifdef ANTIOCH_HAVE_METAPHYSICL
39 #include "metaphysicl/numberarray.h"
40 #endif
41 
42 #ifdef ANTIOCH_HAVE_VEXCL
43 #include "vexcl/vexcl.hpp"
44 #endif
45 
46 #include "antioch/eigen_utils_decl.h"
47 #include "antioch/metaphysicl_utils_decl.h"
48 #include "antioch/valarray_utils_decl.h"
49 
50 #include "antioch/hercourtessen_rate.h"
51 
52 #include "antioch/eigen_utils.h"
53 #include "antioch/metaphysicl_utils.h"
54 #include "antioch/valarray_utils.h"
55 
56 #include <cmath>
57 #include <limits>
58 
59 #ifdef ANTIOCH_HAVE_GRVY
60 #include "grvy.h"
61 
62 GRVY::GRVY_Timer_Class gt;
63 #endif
64 
65 
66 template <typename PairScalars>
67 int check_rate_and_derivative(const PairScalars & rate_exact, const PairScalars & derive_exact,
68  const PairScalars & rate, const PairScalars & derive, const PairScalars & T)
69 {
70  typedef typename Antioch::value_type<PairScalars>::type Scalar;
71  const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 2;
72 
73  int return_flag(0);
74  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
75  {
76  if( abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) > tol )
77  {
78  std::cout << std::scientific << std::setprecision(16)
79  << "Error: Mismatch in rate values." << std::endl
80  << "T = " << T << " K" << std::endl
81  << "rate(T) = " << rate[2*tuple] << std::endl
82  << "rate_exact = " << rate_exact[2*tuple] << std::endl
83  << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) << std::endl
84  << "tolerance = " << tol << std::endl;
85 
86  return_flag = 1;
87  }
88  if( abs( (rate[2*tuple+1] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) > tol )
89  {
90  std::cout << std::scientific << std::setprecision(16)
91  << "Error: Mismatch in rate values." << std::endl
92  << "T = " << T << " K" << std::endl
93  << "rate(T) = " << rate[2*tuple] << std::endl
94  << "rate_exact = " << rate_exact[2*tuple+1] << std::endl
95  << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) << std::endl
96  << "tolerance = " << tol << std::endl;
97 
98  return_flag = 1;
99  }
100  if( abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) > tol )
101  {
102  std::cout << std::scientific << std::setprecision(16)
103  << "Error: Mismatch in rate derivative values." << std::endl
104  << "T = " << T << " K" << std::endl
105  << "drate_dT(T) = " << derive[2*tuple] << std::endl
106  << "derive_exact = " << derive_exact[2*tuple] << std::endl
107  << "relative difference = " << abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) << std::endl
108  << "tolerance = " << tol << std::endl;
109 
110  return_flag = 1;
111  }
112  if( abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) > tol )
113  {
114  std::cout << std::scientific << std::setprecision(16)
115  << "Error: Mismatch in rate derivative values." << std::endl
116  << "T = " << T << " K" << std::endl
117  << "drate_dT(T) = " << derive[2*tuple+1] << std::endl
118  << "derive_exact = " << derive_exact[2*tuple+1] << std::endl
119  << "relative difference = " << abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) << std::endl
120  << "tolerance = " << tol << std::endl;
121 
122  return_flag = 1;
123  }
124 
125  }
126  return return_flag;
127 }
128 
129 template <typename PairScalars>
130 int vectester(const PairScalars& example, const std::string & testname)
131 {
132  using std::abs;
133  using std::pow;
134 
135  typedef typename Antioch::value_type<PairScalars>::type Scalar;
136 
137  const Scalar Cf = 1.4;
138  const Scalar eta = 1.2;
139 
140  Antioch::HercourtEssenRate<Scalar> hercourtessen_rate(Cf,eta);
141 
142  // Construct from example to avoid resizing issues
143  PairScalars T = example;
144  PairScalars rate_exact = example;
145  PairScalars derive_exact = example;
146  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
147  {
148  T[2*tuple] = 1500.1L;
149  T[2*tuple+1] = 1600.1L;
150  rate_exact[2*tuple] = Cf*pow(Scalar(1500.1),eta);
151  rate_exact[2*tuple+1] = Cf*pow(Scalar(1600.1),eta);
152  derive_exact[2*tuple] = eta * Cf * pow(Scalar(1500.1),eta)/Scalar(1500.1);
153  derive_exact[2*tuple+1] = eta * Cf * pow(Scalar(1600.1),eta)/Scalar(1600.1);
154  }
155  Antioch::KineticsConditions<PairScalars> cond(T);
156 
157  int return_flag = 0;
158 
159 #ifdef ANTIOCH_HAVE_GRVY
160  gt.BeginTimer(testname);
161 #endif
162 
163 // KineticsConditions
164  PairScalars rate = hercourtessen_rate(cond);
165  PairScalars derive = hercourtessen_rate.derivative(cond);
166 
167 #ifdef ANTIOCH_HAVE_GRVY
168  gt.EndTimer(testname);
169 #endif
170 
171  return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag;
172 
173 #ifdef ANTIOCH_HAVE_GRVY
174  gt.BeginTimer(testname);
175 #endif
176 
177  hercourtessen_rate.rate_and_derivative(cond,rate,derive);
178 
179 #ifdef ANTIOCH_HAVE_GRVY
180  gt.EndTimer(testname);
181 #endif
182 
183  return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag;
184 
185 // T
186 #ifdef ANTIOCH_HAVE_GRVY
187  gt.BeginTimer(testname);
188 #endif
189  rate = hercourtessen_rate(T);
190  derive = hercourtessen_rate.derivative(T);
191 
192 #ifdef ANTIOCH_HAVE_GRVY
193  gt.EndTimer(testname);
194 #endif
195 
196  return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag;
197 
198 #ifdef ANTIOCH_HAVE_GRVY
199  gt.BeginTimer(testname);
200 #endif
201  hercourtessen_rate.rate_and_derivative(T,rate,derive);
202 
203 #ifdef ANTIOCH_HAVE_GRVY
204  gt.EndTimer(testname);
205 #endif
206 
207  return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag;
208 
209  std::cout << "Hercourt Essen rate: " << hercourtessen_rate << std::endl;
210 
211  return return_flag;
212 }
213 
214 
215 int main()
216 {
217  int returnval = 0;
218 
219  returnval = returnval ||
220  vectester (std::valarray<float>(2*ANTIOCH_N_TUPLES), "valarray<float>");
221  returnval = returnval ||
222  vectester (std::valarray<double>(2*ANTIOCH_N_TUPLES), "valarray<double>");
223  returnval = returnval ||
224  vectester (std::valarray<long double>(2*ANTIOCH_N_TUPLES), "valarray<ld>");
225 #ifdef ANTIOCH_HAVE_EIGEN
226  returnval = returnval ||
227  vectester (Eigen::Array<float, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXf");
228  returnval = returnval ||
229  vectester (Eigen::Array<double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXd");
230  returnval = returnval ||
231  vectester (Eigen::Array<long double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXld");
232 #endif
233 #ifdef ANTIOCH_HAVE_METAPHYSICL
234  returnval = returnval ||
235  vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray<float>");
236  returnval = returnval ||
237  vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray<double>");
238  returnval = returnval ||
239  vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray<ld>");
240 #endif
241 #ifdef ANTIOCH_HAVE_VEXCL
242  vex::Context ctx_f (vex::Filter::All);
243  if (!ctx_f.empty())
244  returnval = returnval ||
245  vectester (vex::vector<float> (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector<float>");
246 
247  vex::Context ctx_d (vex::Filter::DoublePrecision);
248  if (!ctx_d.empty())
249  returnval = returnval ||
250  vectester (vex::vector<double> (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector<double>");
251 #endif
252 
253 #ifdef ANTIOCH_HAVE_GRVY
254  gt.Finalize();
255  gt.Summarize();
256 #endif
257 
258  return returnval;
259 }
Antioch::HercourtEssenRate::derivative
_Cf VectorStateType VectorStateType derivative(const KineticsConditions< StateType, VectorStateType > &T) const ANTIOCH_AUTOFUNC(StateType
std::pow
Antioch::enable_if_c< Antioch::is_valarray< T >::value, typename Antioch::state_type< T >::type >::type pow(const T &in, const T2 &n)
Definition: valarray_utils.h:71
vectester
int vectester(const PairScalars &example, const std::string &testname)
Definition: hercourtessen_rate_vec_unit.C:130
Antioch::HercourtEssenRate
Hercourt-Essen rate equation.
Definition: hercourtessen_rate.h:57
check_rate_and_derivative
int check_rate_and_derivative(const PairScalars &rate_exact, const PairScalars &derive_exact, const PairScalars &rate, const PairScalars &derive, const PairScalars &T)
Definition: hercourtessen_rate_vec_unit.C:67
Antioch::KineticsConditions
This class contains the conditions of the chemistry.
Definition: kinetics_conditions.h:64
Antioch::HercourtEssenRate::rate_and_derivative
_Cf VectorStateType VectorStateType(*) VectorStateType voi rate_and_derivative)(const KineticsConditions< StateType, VectorStateType > &T, StateType &rate, StateType &drate_dT) const
Definition: hercourtessen_rate.h:149
Generated on Thu Jul 7 2016 11:09:46 for antioch-0.4.0 by  doxygen 1.8.8