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
32 #include "antioch_config.h"
33
34 #include <valarray>
35
36 #ifdef ANTIOCH_HAVE_EIGEN
37 #include "Eigen/Dense"
38 #endif
39
40 #ifdef ANTIOCH_HAVE_METAPHYSICL
41 #include "metaphysicl/numberarray.h"
42 #endif
43
44 #ifdef ANTIOCH_HAVE_VEXCL
45 #include "vexcl/vexcl.hpp"
46 #endif
47
48 // C++
49 #include <iostream>
50 #include <cmath>
51 #include <limits>
52
53 // Antioch
54
55 // Declare metaprogramming overloads before they're used
56 #include "antioch/eigen_utils_decl.h"
57 #include "antioch/metaphysicl_utils_decl.h"
58 #include "antioch/valarray_utils_decl.h"
59 #include "antioch/vector_utils_decl.h"
60 #include "antioch/vexcl_utils_decl.h"
61
62 #include "antioch/sutherland_viscosity.h"
63
64 #include "antioch/eigen_utils.h"
65 #include "antioch/metaphysicl_utils.h"
66 #include "antioch/valarray_utils.h"
67 #include "antioch/vexcl_utils.h"
68
69 #ifdef ANTIOCH_HAVE_GRVY
70 #include "grvy.h"
71
72 GRVY::GRVY_Timer_Class gt;
73 #endif
74
75 template <typename Scalar, typename PairScalars>
test_viscosity(const PairScalars mu,const PairScalars mu_exact,const Scalar tol)76 int test_viscosity( const PairScalars mu, const PairScalars mu_exact, const Scalar tol )
77 {
78 using std::abs;
79
80 int return_flag = 0;
81
82 const PairScalars rel_error = abs( (mu - mu_exact)/mu_exact);
83
84 if( Antioch::max(rel_error) > tol )
85 {
86 std::cerr << "Error: Mismatch in viscosity" << std::endl
87 << "mu(T) = " << mu << std::endl
88 << "mu_exact = " << mu_exact << std::endl
89 << "rel_error = " << rel_error << std::endl
90 << "tol = " << tol << std::endl;
91 return_flag = 1;
92 }
93
94 return return_flag;
95 }
96
97
98 template <typename PairScalars>
vectester(const PairScalars & example,const std::string & testname)99 int vectester(const PairScalars& example, const std::string& testname)
100 {
101 typedef typename Antioch::value_type<PairScalars>::type Scalar;
102
103 const Scalar mu_ref = 1.0e-3L;
104 const Scalar T_ref = 300.0L;
105
106 Antioch::SutherlandViscosity<Scalar> mu(mu_ref,T_ref);
107
108 std::cout << mu << std::endl;
109
110 PairScalars T = example;
111 PairScalars mu_exact = example;
112 PairScalars mu_exact2 = example;
113
114 for (unsigned int tuple=0; tuple != ANTIOCH_N_TUPLES; ++tuple)
115 {
116 T[2*tuple ] = 1521.2L;
117 T[2*tuple+1] = 1721.2L;
118
119 // bc with scale=40 gives
120 mu_exact[2*tuple ] = .0325778060534850406481862157435995107036L;
121 mu_exact[2*tuple+1] = .0353295183373055195000058747316029365368L;
122
123 mu_exact2[2*tuple ] = .0959985656417205050367745642313443587197L;
124 mu_exact2[2*tuple+1] = .1047500160115581483776648561664869285592L;
125 }
126
127 int return_flag = 0;
128
129 const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 10;
130
131 #ifdef ANTIOCH_HAVE_GRVY
132 gt.BeginTimer(testname);
133 #endif
134
135 PairScalars muT = mu(T);
136
137 #ifdef ANTIOCH_HAVE_GRVY
138 gt.EndTimer(testname);
139 #endif
140
141 return_flag = test_viscosity( muT, mu_exact, tol );
142
143 const Scalar mu_ref2 = 3.14159e-3L;
144 const Scalar T_ref2 = 420.42L;
145
146 mu.reset_coeffs(mu_ref2,T_ref2);
147
148 #ifdef ANTIOCH_HAVE_GRVY
149 gt.BeginTimer(testname);
150 #endif
151
152 muT = mu(T);
153
154 #ifdef ANTIOCH_HAVE_GRVY
155 gt.EndTimer(testname);
156 #endif
157
158 return_flag = test_viscosity( muT, mu_exact2, tol );
159
160 return return_flag;
161 }
162
163
main()164 int main()
165 {
166 int returnval = 0;
167
168 returnval = returnval ||
169 vectester (std::valarray<float>(2*ANTIOCH_N_TUPLES), "valarray<float>");
170 returnval = returnval ||
171 vectester (std::valarray<double>(2*ANTIOCH_N_TUPLES), "valarray<double>");
172 returnval = returnval ||
173 vectester (std::valarray<long double>(2*ANTIOCH_N_TUPLES), "valarray<ld>");
174 #ifdef ANTIOCH_HAVE_EIGEN
175 returnval = returnval ||
176 vectester (Eigen::Array<float, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXf");
177 returnval = returnval ||
178 vectester (Eigen::Array<double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXd");
179 returnval = returnval ||
180 vectester (Eigen::Array<long double, 2*ANTIOCH_N_TUPLES, 1>(), "Eigen::ArrayXld");
181 #endif
182 #ifdef ANTIOCH_HAVE_METAPHYSICL
183 returnval = returnval ||
184 vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, float> (0), "NumberArray<float>");
185 returnval = returnval ||
186 vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, double> (0), "NumberArray<double>");
187 returnval = returnval ||
188 vectester (MetaPhysicL::NumberArray<2*ANTIOCH_N_TUPLES, long double> (0), "NumberArray<ld>");
189 #endif
190 #ifdef ANTIOCH_HAVE_VEXCL
191 vex::Context ctx_f (vex::Filter::All);
192 if (!ctx_f.empty())
193 returnval = returnval ||
194 vectester (vex::vector<float> (ctx_f, 2*ANTIOCH_N_TUPLES), "vex::vector<float>");
195
196 vex::Context ctx_d (vex::Filter::DoublePrecision);
197 if (!ctx_d.empty())
198 returnval = returnval ||
199 vectester (vex::vector<double> (ctx_d, 2*ANTIOCH_N_TUPLES), "vex::vector<double>");
200 #endif
201
202 #ifdef ANTIOCH_HAVE_GRVY
203 gt.Finalize();
204 gt.Summarize();
205 #endif
206
207 return returnval;
208 }
209