antioch-0.4.0
molecular_binary_diffusion_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 
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 #include "antioch/vexcl_utils_decl.h"
51 #include "antioch/vector_utils_decl.h"
52 
53 #include "antioch/molecular_binary_diffusion.h"
54 #include "antioch/gsl_spliner.h"
55 
56 #include "antioch/eigen_utils.h"
57 #include "antioch/metaphysicl_utils.h"
58 #include "antioch/valarray_utils.h"
59 #include "antioch/vexcl_utils.h"
60 #include "antioch/vector_utils.h"
61 
62 #ifdef ANTIOCH_HAVE_GRVY
63 #include "grvy.h"
64 
65 GRVY::GRVY_Timer_Class gt;
66 #endif
67 
68 #include <cmath>
69 #include <limits>
70 
71 #ifdef ANTIOCH_HAVE_GSL
72 
73 template <typename Scalar, typename Element, typename ElementOrScalar>
74 int test_diff( const Element & dij, const ElementOrScalar & dij_exact, const Scalar & tol, const std::string & words )
75 {
76  using std::abs;
77 
78  int return_flag = 0;
79 
80  const double rel_error = abs( (dij - dij_exact)/dij_exact);
81 
82  if( rel_error > tol )
83  {
84  std::cerr << std::setprecision(15) << std::scientific;
85  std::cerr << "Error: Mismatch in bimolecular coefficient of " << words << std::endl
86  << "Dij(T) = " << dij << std::endl
87  << "Dij_exact = " << dij_exact << std::endl
88  << "rel_error = " << rel_error << std::endl
89  << "tol = " << tol << std::endl;
90  return_flag = 1;
91  }
92 
93  return return_flag;
94 }
95 
96 
97 template <typename PairScalars>
98 int vectester(const PairScalars& example, const std::string& testname)
99 {
100  typedef typename Antioch::value_type<PairScalars>::type Scalar;
101 
102 // N2
103  const Scalar N2_LJ_eps(97.530L);
104  const Scalar N2_LJ_depth(3.621L);
105  const Scalar N2_dipole(0.L);
106  const Scalar N2_polar(1.760L);
107  const Scalar N2_Zrot(4.L);
108  const Scalar N2_mass(28.016e-3L);
109 // CH4
110  const Scalar CH4_LJ_eps(141.400L);
111  const Scalar CH4_LJ_depth(3.746L);
112  const Scalar CH4_dipole(0.L);
113  const Scalar CH4_polar(2.6L);
114  const Scalar CH4_Zrot(13.L);
115  const Scalar CH4_mass(16.043e-3L);
116 // H2O
117  const Scalar H2O_LJ_eps(572.4L);
118  const Scalar H2O_LJ_depth(2.605L);
119  const Scalar H2O_dipole(1.844L);
120  const Scalar H2O_polar(0.L);
121  const Scalar H2O_Zrot(4.L);
122  const Scalar H2O_mass(18.016e-3L);
123 
124  Antioch::TransportSpecies<Scalar> N2(0,N2_LJ_eps,N2_LJ_depth,N2_dipole,N2_polar,N2_Zrot,N2_mass),
125  CH4(1,CH4_LJ_eps,CH4_LJ_depth,CH4_dipole,CH4_polar,CH4_Zrot,CH4_mass),
126  H2O(2,H2O_LJ_eps,H2O_LJ_depth,H2O_dipole,H2O_polar,H2O_Zrot,H2O_mass);
127 
128 
129  Antioch::MolecularBinaryDiffusion<Scalar,Antioch::GSLSpliner>
130  D00(N2,N2), D01(N2,CH4), D02(N2,H2O),
131  D10(CH4,N2), D11(CH4,CH4), D12(CH4,H2O),
132  D20(H2O,N2), D21(H2O,CH4), D22(H2O,H2O);
133 
134 
135  // Construct from example to avoid resizing issues
136  PairScalars T = example;
137  PairScalars cTot = example;
138  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
139  {
140  T[2*tuple] = 1500.1;
141  T[2*tuple+1] = 1600.1;
142  cTot[2*tuple] = 5e-7;
143  cTot[2*tuple+1] = 5e-7;
144  }
145  // bc gives
146  const Scalar D00_exact0 = 3.575629059282712203350417538824313603e3;
147  const Scalar D11_exact0 = 4.415035849326582722446637772820322872e3;
148  const Scalar D22_exact0 = 8.615250909281137767894964155009068175e3;
149  const Scalar D10_exact0 = 4.048999169845892614961423528844221310e3;
150  const Scalar D20_exact0 = 5.468107000169991297723144486054211050e3;
151  const Scalar D12_exact0 = 5.973341001459783642751059656941311432e3;
152 
153  const Scalar D00_exact1 = 3.692886119994007408894312228191265622e3;
154  const Scalar D11_exact1 = 4.559819918938905357528221639097498410e3;
155  const Scalar D22_exact1 = 8.897774342826358536851687124754748934e3;
156  const Scalar D10_exact1 = 4.181779649478152213931731689698388959e3;
157  const Scalar D20_exact1 = 5.647424861129375458731724279131726416e3;
158  const Scalar D12_exact1 = 6.169227206892386749039436637869885641e3;
159 
160 #ifdef ANTIOCH_HAVE_GRVY
161  gt.BeginTimer(testname);
162 #endif
163 
164  PairScalars d00 = D00(T,cTot) * D00.Stockmayer(T);
165  PairScalars d11 = D11(T,cTot) * D11.Stockmayer(T);
166  PairScalars d22 = D22(T,cTot) * D22.Stockmayer(T);
167 
168  PairScalars d10 = D10(T,cTot) * D10.Stockmayer(T);
169  PairScalars d01 = D01(T,cTot) * D01.Stockmayer(T);
170  PairScalars d20 = D20(T,cTot) * D20.Stockmayer(T);
171  PairScalars d02 = D02(T,cTot) * D02.Stockmayer(T);
172  PairScalars d21 = D21(T,cTot) * D21.Stockmayer(T);
173  PairScalars d12 = D12(T,cTot) * D12.Stockmayer(T);
174 
175  int return_flag = 0;
176 
177 #ifdef ANTIOCH_HAVE_GRVY
178  gt.EndTimer(testname);
179 #endif
180 
181  const Scalar tol = (std::numeric_limits<Scalar>::epsilon()*10 < 5e-16)?5e-16:
182  std::numeric_limits<Scalar>::epsilon()*10;
183 
184  for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
185  {
186  return_flag = test_diff( d10[2*tuple] , d01[2*tuple], tol, "N2 - CH4 symetry" ) || return_flag;
187  return_flag = test_diff( d20[2*tuple] , d02[2*tuple], tol, "N2 - H2O symetry" ) || return_flag;
188  return_flag = test_diff( d12[2*tuple] , d21[2*tuple], tol, "H2O - CH4 symetry" ) || return_flag;
189 //
190  return_flag = test_diff( d00[2*tuple] , D00_exact0, tol, "N2 - N2" ) || return_flag;
191  return_flag = test_diff( d11[2*tuple] , D11_exact0, tol, "CH4 - CH4" ) || return_flag;
192  return_flag = test_diff( d22[2*tuple] , D22_exact0, tol, "H2O - H2O" ) || return_flag;
193  return_flag = test_diff( d10[2*tuple] , D10_exact0, tol, "N2 - CH4" ) || return_flag;
194  return_flag = test_diff( d12[2*tuple] , D12_exact0, tol, "CH4 - H2O" ) || return_flag;
195  return_flag = test_diff( d20[2*tuple] , D20_exact0, tol, "N2 - H2O" ) || return_flag;
196 //
197  return_flag = test_diff( d00[2*tuple+1] , D00_exact1, tol, "N2 - N2" ) || return_flag;
198  return_flag = test_diff( d11[2*tuple+1] , D11_exact1, tol, "CH4 - CH4" ) || return_flag;
199  return_flag = test_diff( d22[2*tuple+1] , D22_exact1, tol, "H2O - H2O" ) || return_flag;
200  return_flag = test_diff( d10[2*tuple+1] , D10_exact1, tol, "N2 - CH4" ) || return_flag;
201  return_flag = test_diff( d12[2*tuple+1] , D12_exact1, tol, "CH4 - H2O" ) || return_flag;
202  return_flag = test_diff( d20[2*tuple+1] , D20_exact1, tol, "N2 - H2O" ) || return_flag;
203  }
204 
205  return return_flag;
206 }
207 #endif // ANTIOCH_HAVE_GSL
208 
209 int main()
210 {
211  int returnval = 0;
212 
213 #ifdef ANTIOCH_HAVE_GSL
214 
215  returnval = returnval ||
216  vectester (std::valarray<float>(2*ANTIOCH_N_TUPLES), "valarray<float>");
217  returnval = returnval ||
218  vectester (std::valarray<double>(2*ANTIOCH_N_TUPLES), "valarray<double>");
219  returnval = returnval ||
220  vectester (std::valarray<long double>(2*ANTIOCH_N_TUPLES), "valarray<ld>");
221 #ifdef ANTIOCH_HAVE_EIGEN
222  returnval = returnval ||
223  vectester (Eigen::Array<float, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXf");
224  returnval = returnval ||
225  vectester (Eigen::Array<double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXd");
226  returnval = returnval ||
227  vectester (Eigen::Array<long double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXld");
228 #endif
229 #ifdef ANTIOCH_HAVE_METAPHYSICL
230  returnval = returnval ||
231  vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray<float>");
232  returnval = returnval ||
233  vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray<double>");
234  returnval = returnval ||
235  vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray<ld>");
236 #endif
237 #ifdef ANTIOCH_HAVE_VEXCL
238  vex::Context ctx_f (vex::Filter::All);
239  if (!ctx_f.empty())
240  returnval = returnval ||
241  vectester (vex::vector<float> (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector<float>");
242 
243  vex::Context ctx_d (vex::Filter::DoublePrecision);
244  if (!ctx_d.empty())
245  returnval = returnval ||
246  vectester (vex::vector<double> (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector<double>");
247 #endif
248 
249 #ifdef ANTIOCH_HAVE_GRVY
250  gt.Finalize();
251  gt.Summarize();
252 #endif
253 
254 #else // ANTIOCH_HAVE_GSL
255  // 77 return code tells Automake we skipped this.
256  returnval = 77;
257 #endif
258 
259  return returnval;
260 }
Antioch::TransportSpecies
Class to encapsulate data relevant for transport for each chemical species.
Definition: transport_species.h:58
vectester
int vectester(const PairScalars &example, const std::string &testname)
Definition: berthelot_rate_vec_unit.C:129
Generated on Thu Jul 7 2016 11:09:46 for antioch-0.4.0 by  doxygen 1.8.8