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 #include "antioch/vexcl_utils_decl.h"
50
51 #include "antioch/hercourtessen_rate.h"
52
53 #include "antioch/eigen_utils.h"
54 #include "antioch/metaphysicl_utils.h"
55 #include "antioch/valarray_utils.h"
56 #include "antioch/vexcl_utils.h"
57
58 #include <cmath>
59 #include <limits>
60
61 #ifdef ANTIOCH_HAVE_GRVY
62 #include "grvy.h"
63
64 GRVY::GRVY_Timer_Class gt;
65 #endif
66
67
68 template <typename PairScalars>
check_rate_and_derivative(const PairScalars & rate_exact,const PairScalars & derive_exact,const PairScalars & rate,const PairScalars & derive,const PairScalars & T)69 int check_rate_and_derivative(const PairScalars & rate_exact, const PairScalars & derive_exact,
70 const PairScalars & rate, const PairScalars & derive, const PairScalars & T)
71 {
72 typedef typename Antioch::value_type<PairScalars>::type Scalar;
73 const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 2;
74
75 int return_flag(0);
76
77 using std::abs;
78
79 for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
80 {
81 if( abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) > tol )
82 {
83 std::cout << std::scientific << std::setprecision(16)
84 << "Error: Mismatch in rate values." << std::endl
85 << "T = " << T << " K" << std::endl
86 << "rate(T) = " << rate[2*tuple] << std::endl
87 << "rate_exact = " << rate_exact[2*tuple] << std::endl
88 << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple])/rate_exact[2*tuple] ) << std::endl
89 << "tolerance = " << tol << std::endl;
90
91 return_flag = 1;
92 }
93 if( abs( (rate[2*tuple+1] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) > tol )
94 {
95 std::cout << std::scientific << std::setprecision(16)
96 << "Error: Mismatch in rate values." << std::endl
97 << "T = " << T << " K" << std::endl
98 << "rate(T) = " << rate[2*tuple] << std::endl
99 << "rate_exact = " << rate_exact[2*tuple+1] << std::endl
100 << "relative difference = " << abs( (rate[2*tuple] - rate_exact[2*tuple+1])/rate_exact[2*tuple+1] ) << std::endl
101 << "tolerance = " << tol << std::endl;
102
103 return_flag = 1;
104 }
105 if( abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) > tol )
106 {
107 std::cout << std::scientific << std::setprecision(16)
108 << "Error: Mismatch in rate derivative values." << std::endl
109 << "T = " << T << " K" << std::endl
110 << "drate_dT(T) = " << derive[2*tuple] << std::endl
111 << "derive_exact = " << derive_exact[2*tuple] << std::endl
112 << "relative difference = " << abs( (derive[2*tuple] - derive_exact[2*tuple])/derive_exact[2*tuple] ) << std::endl
113 << "tolerance = " << tol << std::endl;
114
115 return_flag = 1;
116 }
117 if( abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) > tol )
118 {
119 std::cout << std::scientific << std::setprecision(16)
120 << "Error: Mismatch in rate derivative values." << std::endl
121 << "T = " << T << " K" << std::endl
122 << "drate_dT(T) = " << derive[2*tuple+1] << std::endl
123 << "derive_exact = " << derive_exact[2*tuple+1] << std::endl
124 << "relative difference = " << abs( (derive[2*tuple+1] - derive_exact[2*tuple+1])/derive_exact[2*tuple+1] ) << std::endl
125 << "tolerance = " << tol << std::endl;
126
127 return_flag = 1;
128 }
129
130 }
131 return return_flag;
132 }
133
134 template <typename PairScalars>
vectester(const PairScalars & example,const std::string & testname)135 int vectester(const PairScalars& example, const std::string & testname)
136 {
137 using std::abs;
138 using std::pow;
139
140 typedef typename Antioch::value_type<PairScalars>::type Scalar;
141
142 const Scalar Cf = 1.4;
143 const Scalar eta = 1.2;
144
145 Antioch::HercourtEssenRate<Scalar> hercourtessen_rate(Cf,eta);
146
147 // Construct from example to avoid resizing issues
148 PairScalars T = example;
149 PairScalars rate_exact = example;
150 PairScalars derive_exact = example;
151 for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
152 {
153 T[2*tuple] = 1500.1L;
154 T[2*tuple+1] = 1600.1L;
155 rate_exact[2*tuple] = Cf*pow(Scalar(1500.1),eta);
156 rate_exact[2*tuple+1] = Cf*pow(Scalar(1600.1),eta);
157 derive_exact[2*tuple] = eta * Cf * pow(Scalar(1500.1),eta)/Scalar(1500.1);
158 derive_exact[2*tuple+1] = eta * Cf * pow(Scalar(1600.1),eta)/Scalar(1600.1);
159 }
160 Antioch::KineticsConditions<PairScalars> cond(T);
161
162 int return_flag = 0;
163
164 #ifdef ANTIOCH_HAVE_GRVY
165 gt.BeginTimer(testname);
166 #endif
167
168 // KineticsConditions
169 PairScalars rate = hercourtessen_rate(cond);
170 PairScalars derive = hercourtessen_rate.derivative(cond);
171
172 #ifdef ANTIOCH_HAVE_GRVY
173 gt.EndTimer(testname);
174 #endif
175
176 return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag;
177
178 #ifdef ANTIOCH_HAVE_GRVY
179 gt.BeginTimer(testname);
180 #endif
181
182 hercourtessen_rate.rate_and_derivative(cond,rate,derive);
183
184 #ifdef ANTIOCH_HAVE_GRVY
185 gt.EndTimer(testname);
186 #endif
187
188 return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag;
189
190 // T
191 #ifdef ANTIOCH_HAVE_GRVY
192 gt.BeginTimer(testname);
193 #endif
194 rate = hercourtessen_rate(T);
195 derive = hercourtessen_rate.derivative(T);
196
197 #ifdef ANTIOCH_HAVE_GRVY
198 gt.EndTimer(testname);
199 #endif
200
201 return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag;
202
203 #ifdef ANTIOCH_HAVE_GRVY
204 gt.BeginTimer(testname);
205 #endif
206 hercourtessen_rate.rate_and_derivative(T,rate,derive);
207
208 #ifdef ANTIOCH_HAVE_GRVY
209 gt.EndTimer(testname);
210 #endif
211
212 return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag;
213
214 std::cout << "Hercourt Essen rate: " << hercourtessen_rate << std::endl;
215
216 return return_flag;
217 }
218
219
main()220 int main()
221 {
222 int returnval = 0;
223
224 returnval = returnval ||
225 vectester (std::valarray<float>(2*ANTIOCH_N_TUPLES), "valarray<float>");
226 returnval = returnval ||
227 vectester (std::valarray<double>(2*ANTIOCH_N_TUPLES), "valarray<double>");
228 returnval = returnval ||
229 vectester (std::valarray<long double>(2*ANTIOCH_N_TUPLES), "valarray<ld>");
230 #ifdef ANTIOCH_HAVE_EIGEN
231 returnval = returnval ||
232 vectester (Eigen::Array<float, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXf");
233 returnval = returnval ||
234 vectester (Eigen::Array<double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXd");
235 returnval = returnval ||
236 vectester (Eigen::Array<long double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXld");
237 #endif
238 #ifdef ANTIOCH_HAVE_METAPHYSICL
239 returnval = returnval ||
240 vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray<float>");
241 returnval = returnval ||
242 vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray<double>");
243 returnval = returnval ||
244 vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray<ld>");
245 #endif
246 #ifdef ANTIOCH_HAVE_VEXCL
247 vex::Context ctx_f (vex::Filter::All);
248 if (!ctx_f.empty())
249 returnval = returnval ||
250 vectester (vex::vector<float> (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector<float>");
251
252 vex::Context ctx_d (vex::Filter::DoublePrecision);
253 if (!ctx_d.empty())
254 returnval = returnval ||
255 vectester (vex::vector<double> (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector<double>");
256 #endif
257
258 #ifdef ANTIOCH_HAVE_GRVY
259 gt.Finalize();
260 gt.Summarize();
261 #endif
262
263 return returnval;
264 }
265