antioch-0.4.0
constant_rate_vec_unit.C
Go to the documentation of this file.
1 
2 //-----------------------------------------------------------------------bl-
3 //--------------------------------------------------------------------------
4 //
5 // Antioch - A Gas Dynamics Thermochemistry Library
6 //
7 // Copyright (C) 2014-2016 Paul T. Bauman, Benjamin S. Kirk,
8 // Sylvain Plessis, Roy H. Stonger
9 //
10 // Copyright (C) 2013 The PECOS Development Team
11 //
12 // This library is free software; you can redistribute it and/or
13 // modify it under the terms of the Version 2.1 GNU Lesser General
14 // Public License as published by the Free Software Foundation.
15 //
16 // This library is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 // Lesser General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License along with this library; if not, write to the Free Software
23 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
24 // Boston, MA 02110-1301 USA
25 //
26 //-----------------------------------------------------------------------el-
27 
28 // valarray has to be declared before Antioch or gcc can't find the
29 // right versions of exp() and pow() to use??
30 
31 #include "antioch_config.h"
32 
33 #include <valarray>
34 
35 #ifdef ANTIOCH_HAVE_EIGEN
36 #include "Eigen/Dense"
37 #endif
38 
39 #ifdef ANTIOCH_HAVE_METAPHYSICL
40 #include "metaphysicl/numberarray.h"
41 #endif
42 
43 #ifdef ANTIOCH_HAVE_VEXCL
44 #include "vexcl/vexcl.hpp"
45 #endif
46 
47 #include "antioch/eigen_utils_decl.h"
48 #include "antioch/metaphysicl_utils_decl.h"
49 #include "antioch/valarray_utils_decl.h"
50 
51 #include "antioch/constant_rate.h"
52 
53 #include "antioch/eigen_utils.h"
54 #include "antioch/metaphysicl_utils.h"
55 #include "antioch/valarray_utils.h"
56 
57 #include <cmath>
58 #include <limits>
59 
60 #ifdef ANTIOCH_HAVE_GRVY
61 #include "grvy.h"
62 
63 GRVY::GRVY_Timer_Class gt;
64 #endif
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[2*tuple] << " 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[2*tuple+1] << " 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[2*tuple] << " 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[2*tuple+1] << " 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 
134  typedef typename Antioch::value_type<PairScalars>::type Scalar;
135 
136  const Scalar Cf = 1.4L;
137 
138  Antioch::ConstantRate<Scalar> constant_rate(Cf);
139 
140  // Construct from example to avoid resizing issues
141  PairScalars T = example;
142  PairScalars rate_exact = example;
143  PairScalars derive_exact = Antioch::zero_clone(example);
144  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
145  {
146  T[2*tuple] = 1500.1L;
147  T[2*tuple+1] = 1600.1L;
148  rate_exact[2*tuple] = Cf;
149  rate_exact[2*tuple+1] = Cf;
150  }
151  Antioch::KineticsConditions<PairScalars> cond(T);
152 
153  int return_flag = 0;
154 
155 #ifdef ANTIOCH_HAVE_GRVY
156  gt.BeginTimer(testname);
157 #endif
158 
159 // KineticsConditions
160  PairScalars rate = constant_rate(cond);//Antioch::zero_clone(T.T());
161  PairScalars derive = constant_rate.derivative(cond);//Antioch::zero_clone(T.T());
162 
163 #ifdef ANTIOCH_HAVE_GRVY
164  gt.EndTimer(testname);
165 #endif
166 
167  return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag;
168 
169 #ifdef ANTIOCH_HAVE_GRVY
170  gt.BeginTimer(testname);
171 #endif
172 
173  constant_rate.rate_and_derivative(cond,rate,derive);
174 
175 #ifdef ANTIOCH_HAVE_GRVY
176  gt.EndTimer(testname);
177 #endif
178 
179  return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag;
180 
181 // T
182 
183 #ifdef ANTIOCH_HAVE_GRVY
184  gt.BeginTimer(testname);
185 #endif
186 
187  rate = constant_rate(T);
188  derive = constant_rate.derivative(T);
189 
190 #ifdef ANTIOCH_HAVE_GRVY
191  gt.EndTimer(testname);
192 #endif
193 
194  return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag;
195 
196 #ifdef ANTIOCH_HAVE_GRVY
197  gt.BeginTimer(testname);
198 #endif
199 
200  constant_rate.rate_and_derivative(T,rate,derive);
201 
202 #ifdef ANTIOCH_HAVE_GRVY
203  gt.EndTimer(testname);
204 #endif
205 
206  return_flag = check_rate_and_derivative(rate_exact, derive_exact, rate, derive,T) || return_flag;
207 
208  return return_flag;
209 }
210 
211 
212 int main()
213 {
214  int returnval = 0;
215 
216  returnval = returnval ||
217  vectester (std::valarray<float>(2*ANTIOCH_N_TUPLES), "valarray<float>");
218  returnval = returnval ||
219  vectester (std::valarray<double>(2*ANTIOCH_N_TUPLES), "valarray<double>");
220  returnval = returnval ||
221  vectester (std::valarray<long double>(2*ANTIOCH_N_TUPLES), "valarray<ld>");
222 #ifdef ANTIOCH_HAVE_EIGEN
223  returnval = returnval ||
224  vectester (Eigen::Array<float, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXf");
225  returnval = returnval ||
226  vectester (Eigen::Array<double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXd");
227  returnval = returnval ||
228  vectester (Eigen::Array<long double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXld");
229 #endif
230 #ifdef ANTIOCH_HAVE_METAPHYSICL
231  returnval = returnval ||
232  vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray<float>");
233  returnval = returnval ||
234  vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray<double>");
235  returnval = returnval ||
236  vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray<ld>");
237 #endif
238 #ifdef ANTIOCH_HAVE_VEXCL
239  vex::Context ctx_f (vex::Filter::All);
240  if (!ctx_f.empty())
241  returnval = returnval ||
242  vectester (vex::vector<float> (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector<float>");
243 
244  vex::Context ctx_d (vex::Filter::DoublePrecision);
245  if (!ctx_d.empty())
246  returnval = returnval ||
247  vectester (vex::vector<double> (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector<double>");
248 #endif
249 
250 #ifdef ANTIOCH_HAVE_GRVY
251  gt.Finalize();
252  gt.Summarize();
253 #endif
254 
255  return returnval;
256 }
vectester
int vectester(const PairScalars &example, const std::string &testname)
Definition: constant_rate_vec_unit.C:130
Antioch::ConstantRate
Constant rate equation.
Definition: constant_rate.h:55
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: constant_rate_vec_unit.C:67
main
int main()
Definition: constant_rate_vec_unit.C:212
Antioch::ConstantRate::derivative
VectorStateType VectorStateType derivative(const KineticsConditions< StateType, VectorStateType > &cond) const ANTIOCH_AUTOFUNC(StateType
Antioch::ConstantRate::rate_and_derivative
VectorStateType VectorStateType VectorStateType void rate_and_derivative(const KineticsConditions< StateType, VectorStateType > &cond, StateType &rate, StateType &drate_dT) const
Definition: constant_rate.h:211
Antioch::zero_clone
_Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > zero_clone(const _Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > &ex)
Definition: eigen_utils.h:145
Antioch::KineticsConditions
This class contains the conditions of the chemistry.
Definition: kinetics_conditions.h:64
Generated on Thu Jul 7 2016 11:09:46 for antioch-0.4.0 by  doxygen 1.8.8