antioch-0.4.0
sigma_bin_converter_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 #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/vector_utils_decl.h"
48 #include "antioch/eigen_utils_decl.h"
49 #include "antioch/metaphysicl_utils_decl.h"
50 #include "antioch/valarray_utils_decl.h"
51 #include "antioch/vexcl_utils_decl.h"
52 
53 #include "antioch/sigma_bin_converter.h"
54 #include "antioch/antioch_asserts.h"
55 #include "antioch/physical_constants.h"
56 
57 #include "antioch/vector_utils.h"
58 #include "antioch/eigen_utils.h"
59 #include "antioch/metaphysicl_utils.h"
60 #include "antioch/valarray_utils.h"
61 #include "antioch/vexcl_utils.h"
62 
63 #ifdef ANTIOCH_HAVE_GRVY
64 #include "grvy.h"
65 
66 GRVY::GRVY_Timer_Class gt;
67 #endif
68 
69 // C++
70 #include <cmath>
71 #include <limits>
72 #include <string>
73 #include <iostream>
74 #include <iomanip>
75 #include <fstream>
76 
77 /*
78  x[0] = 1.; y[0] = 1.;
79  x[1] = 2.; y[1] = 2.;
80  x[2] = 3.; y[2] = 5.;
81  x[3] = 4.; y[3] = 8.;
82  x[4] = 5.; y[4] = 6.;
83  x[5] = 6.; y[5] = 10.;
84  x[6] = 7.; y[6] = 7.;
85  x[7] = 8.; y[7] = 4.;
86  x[8] = 9.; y[8] = 2.;
87  x[9] = 10.; y[9] = 0.3; // right stairs, this value is useless
88 */
89 template <typename PairScalars, typename VectorPairScalar>
90 void make_custom(unsigned int choice, const PairScalars & ex, VectorPairScalar & x, VectorPairScalar & y)
91 {
92  // sum_{bin} Delta x * y
93  switch(choice)
94  {
95  case(0): // 9 bin contained in ref -> [2.5;8.5] within [1;10]
96  {
97  x.resize(9,ex); y.resize(9,ex);
98  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
99  {
100  x[0][2*tuple] = 2.50L; y[0][2*tuple] = 2.25L/0.75L; // 0.50][2 + 0.25][5
101  x[1][2*tuple] = 3.25L; y[1][2*tuple] = 3.75L/0.75L; // 0.75][5
102  x[2][2*tuple] = 4.00L; y[2][2*tuple] = 6.00L/0.75L; // 0.75][8
103  x[3][2*tuple] = 4.75L; y[3][2*tuple] = 5.00L/0.75L; // 0.25][8 + 0.5 ][6
104  x[4][2*tuple] = 5.50L; y[4][2*tuple] = 5.50L/0.75L; // 0.25][6 + 0.25][10
105  x[5][2*tuple] = 6.25L; y[5][2*tuple] = 7.50L/0.75L; // 0.75][10
106  x[6][2*tuple] = 7.00L; y[6][2*tuple] = 5.25L/0.75L; // 0.75][7
107  x[7][2*tuple] = 7.75L; y[7][2*tuple] = 3.75L/0.75L; // 0.25][7 + 0.5][4
108  x[8][2*tuple] = 8.50L; y[8][2*tuple] = 0.00L/0.75L; // 0. (right stairs)
109 
110  // 9 bin outside ref -> [0;12] containing [1;10]
111 
112  x[0][2*tuple + 1] = 0.00L; y[0][2*tuple + 1] = 0.50L/1.50L; // 0.5][1
113  x[1][2*tuple + 1] = 1.50L; y[1][2*tuple + 1] = 2.50L/1.50L; // 0.5][1 + 2.0][1
114  x[2][2*tuple + 1] = 3.00L; y[2][2*tuple + 1] = 9.00L/1.50L; // 1.0][5 + 0.5][8
115  x[3][2*tuple + 1] = 4.50L; y[3][2*tuple + 1] = 10.0L/1.50L; // 0.5][8 + 1.0][6
116  x[4][2*tuple + 1] = 6.00L; y[4][2*tuple + 1] = 13.5L/1.50L; // 1.0][10 + 0.5][7
117  x[5][2*tuple + 1] = 7.50L; y[5][2*tuple + 1] = 7.50L/1.50L; // 0.50][7 + 1.0][4
118  x[6][2*tuple + 1] = 9.00L; y[6][2*tuple + 1] = 2.00L/1.50L; // 1.0][2
119  x[7][2*tuple + 1] = 10.5L; y[7][2*tuple + 1] = 0.00L/1.50L; // 0
120  x[8][2*tuple + 1] = 12.0L; y[8][2*tuple + 1] = 0.00L/1.50L; // 0. (right stairs)
121  }
122  break;
123  }
124  case(1): // 5 bins beyond only min -> [-1;3.8] in [0;10]
125  {
126  x.resize(5,ex); y.resize(5,ex);
127  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
128  {
129  x[0][2*tuple] = -1.00L; y[0][2*tuple] = 0.00L/1.20L; // 0
130  x[1][2*tuple] = 0.20L; y[1][2*tuple] = 0.40L/1.20L; // 0.4][1
131  x[2][2*tuple] = 1.40L; y[2][2*tuple] = 1.80L/1.20L; // 0.6][1 + 0.6][2
132  x[3][2*tuple] = 2.60L; y[3][2*tuple] = 4.80L/1.20L; // 0.4][2 + 0.8][5
133  x[4][2*tuple] = 3.80L; y[4][2*tuple] = 00.0L/1.20L; // 0. (right stairs)
134 
135  // 6 bins beyond only max -> [2;10.75] in [0;10]
136  x[0][2*tuple + 1] = 2.0000L; y[0][2*tuple + 1] = 8.50L/2.1875L; // 1.0000][2 + 1.0000][5 + 0.1875][8
137  x[1][2*tuple + 1] = 4.1875L; y[1][2*tuple + 1] = 16.25L/2.1875L; // 0.8125][8 + 1.0000][6 + 0.3750][10
138  x[2][2*tuple + 1] = 6.3750L; y[2][2*tuple + 1] = 15.50L/2.1875L; // 0.6250][10 + 1.0000][7 + 0.5625][4
139  x[3][2*tuple + 1] = 8.5625L; y[3][2*tuple + 1] = 3.75L/2.1875L; // 0.4375][4 + 1.0000][2 + 0
140  x[4][2*tuple + 1] = 10.7500L; y[4][2*tuple + 1] = 0.00L/2.1875L; // 0. (right stairs)
141  }
142  break;
143  }
144  }
145 
146 }
147 
148 template <typename VectorScalar>
149 void make_reference(VectorScalar & x, VectorScalar & y)
150 {
151  x.resize(10,0);
152  y.resize(10,0);
153  for(unsigned int i = 0; i < 10; i++)
154  {
155  x[i] = ((typename Antioch::value_type<VectorScalar>::type)(i+1));
156  }
157  y[0] = 1.L;
158  y[1] = 2.L;
159  y[2] = 5.L;
160  y[3] = 8.L;
161  y[4] = 6.L;
162  y[5] = 10.L;
163  y[6] = 7.L;
164  y[7] = 4.L;
165  y[8] = 2.L;
166  y[9] = 0.3L;
167 
168 }
169 
170 
171 template <typename PairScalars>
172 int vectester(const PairScalars& example, const std::string& testname)
173 {
174  typedef typename Antioch::value_type<PairScalars>::type Scalar;
175 
176  std::vector<Scalar> bin_ref_x,bin_ref_y;
177 
178  make_reference(bin_ref_x,bin_ref_y);
179 
180  Antioch::SigmaBinConverter<std::vector<Scalar> > bin;
181 
182  const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 10;
183 
184  int return_flag = 0;
185 
186 #ifdef ANTIOCH_HAVE_GRVY
187  gt.BeginTimer(testname);
188 #endif
189 
190  // 2 * 2 cases:
191  // - custom inside ref, ref inside custom
192  // - custom beyond only min ref, custom beyond only max ref
193  for(unsigned int i = 0; i < 2; i++)
194  {
195  std::vector<PairScalars> bin_custom_x, exact_sol_y;
196  make_custom(i,example,bin_custom_x,exact_sol_y);
197  std::vector<PairScalars> bin_custom_y(bin_custom_x.size(),Antioch::zero_clone(bin_custom_x[0]));
198 
199  bin.y_on_custom_grid(bin_ref_x,bin_ref_y,
200  bin_custom_x,bin_custom_y);
201 
202 
203  for(unsigned int il = 0; il < bin_custom_x.size() - 1; il++)
204  {
205  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
206  {
207 //tuple
208 
209  Scalar dist = Antioch::if_else(exact_sol_y[il][2*tuple] < tol,
210  std::abs(bin_custom_y[il][2*tuple] - exact_sol_y[il][2*tuple]),
211  std::abs(bin_custom_y[il][2*tuple] - exact_sol_y[il][2*tuple])/exact_sol_y[il][2*tuple]);
212  if( dist > tol )
213  {
214  std::cout << std::scientific << std::setprecision(16)
215  << "Error: Mismatch in bin values for " << testname << std::endl
216  << "case (" << bin_custom_x[il][2*tuple] << ";" << bin_custom_x[il + 1][2*tuple] << ")" << std::endl
217  << "bin = " << bin_custom_y[il][2*tuple] << std::endl
218  << "bin_exact = " << exact_sol_y[il][2*tuple] << std::endl
219  << "relative error = " << dist << std::endl
220  << "tolerance = " << tol << std::endl;
221 
222  return_flag = 1;
223  }
224 
225 //tuple + 1
226  dist = Antioch::if_else(exact_sol_y[il][2*tuple + 1] < tol,
227  std::abs(bin_custom_y[il][2*tuple + 1] - exact_sol_y[il][2*tuple + 1]),
228  std::abs(bin_custom_y[il][2*tuple + 1] - exact_sol_y[il][2*tuple + 1])/exact_sol_y[il][2*tuple + 1]);
229  if( dist > tol )
230  {
231  std::cout << std::scientific << std::setprecision(16)
232  << "Error: Mismatch in bin values for " << testname << std::endl
233  << "case (" << bin_custom_x[il][2*tuple + 1] << ";" << bin_custom_x[il+1][2*tuple + 1] << ")" << std::endl
234  << "bin = " << bin_custom_y[il][2*tuple + 1] << std::endl
235  << "bin_exact = " << exact_sol_y[il][2*tuple + 1] << std::endl
236  << "relative error = " << dist << std::endl
237  << "tolerance = " << tol << std::endl;
238 
239  return_flag = 1;
240  }
241 
242 
243  }
244  }
245  }
246 
247 #ifdef ANTIOCH_HAVE_GRVY
248  gt.EndTimer(testname);
249 #endif
250 
251  return return_flag;
252 }
253 
254 int main()
255 {
256 
257  int returnval(0);
258 
259  returnval = returnval ||
260  vectester (std::valarray<float>(2*ANTIOCH_N_TUPLES), "valarray<float>");
261  returnval = returnval ||
262  vectester (std::valarray<double>(2*ANTIOCH_N_TUPLES), "valarray<double>");
263  returnval = returnval ||
264  vectester (std::valarray<long double>(2*ANTIOCH_N_TUPLES), "valarray<ld>");
265 
266 #ifdef ANTIOCH_HAVE_EIGEN
267  returnval = returnval ||
268  vectester (Eigen::Array<float, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXf");
269  returnval = returnval ||
270  vectester (Eigen::Array<double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXd");
271  returnval = returnval ||
272  vectester (Eigen::Array<long double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXld");
273 #endif
274 #ifdef ANTIOCH_HAVE_METAPHYSICL
275  returnval = returnval ||
276  vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray<float>");
277  returnval = returnval ||
278  vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray<double>");
279  returnval = returnval ||
280  vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray<ld>");
281 #endif
282 #ifdef ANTIOCH_HAVE_VEXCL
283  vex::Context ctx_f (vex::Filter::All);
284  if (!ctx_f.empty())
285  returnval = returnval ||
286  vectester (vex::vector<float> (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector<float>");
287 
288  vex::Context ctx_d (vex::Filter::DoublePrecision);
289  if (!ctx_d.empty())
290  returnval = returnval ||
291  vectester (vex::vector<double> (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector<double>");
292 #endif
293 
294 #ifdef ANTIOCH_HAVE_GRVY
295  gt.Finalize();
296  gt.Summarize();
297 #endif
298 
299  return returnval;
300 
301 }
make_reference
void make_reference(VectorScalar &x, VectorScalar &y)
Definition: sigma_bin_converter_vec_unit.C:149
vectester
int vectester(const PairScalars &example, const std::string &testname)
Definition: sigma_bin_converter_vec_unit.C:172
Antioch::if_else
enable_if_c< is_eigen< T1 >::value &&is_eigen< T2 >::value, typename state_type< T1 >::type >::type if_else(const Condition &condition, const T1 &if_true, const T2 &if_false)
Definition: eigen_utils.h:250
make_custom
void make_custom(unsigned int choice, const PairScalars &ex, VectorPairScalar &x, VectorPairScalar &y)
Definition: sigma_bin_converter_vec_unit.C:90
Antioch::SigmaBinConverter::y_on_custom_grid
void y_on_custom_grid(const VectorCoeffType &x_old, const VectorCoeffType &y_old, const VectorStateType &x_new, VectorStateType &y_new) const
Definition: sigma_bin_converter.h:79
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
Generated on Thu Jul 7 2016 11:09:47 for antioch-0.4.0 by  doxygen 1.8.8