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>
test_diff(const Element & dij,const ElementOrScalar & dij_exact,const Scalar & tol,const std::string & words)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>
vectester(const PairScalars & example,const std::string & testname)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
main()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 }
261