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 // C++
28 #include <limits>
29 #include <vector>
30 // Antioch
31 #include "antioch/berthelot_rate.h"
32
33 template <typename Scalar>
check_rate_and_derivative(const Scalar & rate_exact,const Scalar & derive_exact,const Scalar & rate,const Scalar & derive,const Scalar & T)34 int check_rate_and_derivative(const Scalar & rate_exact, const Scalar & derive_exact,
35 const Scalar & rate, const Scalar & derive, const Scalar & T)
36 {
37 const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 2;
38 int return_flag(0);
39
40 using std::abs;
41
42 if( abs( (rate - rate_exact)/rate_exact ) > tol )
43 {
44 std::cout << std::scientific << std::setprecision(16)
45 << "Error: Mismatch in rate values." << std::endl
46 << "T = " << T << " K" << std::endl
47 << "rate(T) = " << rate << std::endl
48 << "rate_exact = " << rate_exact << std::endl
49 << "relative difference = " << abs( (rate - rate_exact)/rate_exact ) << std::endl
50 << "tolerance = " << tol << std::endl;
51
52 return_flag = 1;
53 }
54 if( abs( (derive - derive_exact)/derive_exact ) > tol )
55 {
56 std::cout << std::scientific << std::setprecision(16)
57 << "Error: Mismatch in rate derivative values." << std::endl
58 << "T = " << T << " K" << std::endl
59 << "drate_dT(T) = " << derive << std::endl
60 << "derive_exact = " << derive_exact << std::endl
61 << "relative difference = " << abs( (derive - derive_exact)/derive_exact ) << std::endl
62 << "tolerance = " << tol << std::endl;
63
64 return_flag = 1;
65 }
66
67 return return_flag;
68 }
69
70 template <typename Scalar>
test_values(const Scalar & Cf,const Scalar & D,const Antioch::BerthelotRate<Scalar> & berthelot_rate)71 int test_values(const Scalar & Cf, const Scalar & D, const Antioch::BerthelotRate<Scalar> & berthelot_rate)
72 {
73 using std::abs;
74 using std::exp;
75 int return_flag = 0;
76
77 for(Scalar T = 300.1; T <= 2500.1; T += 10.)
78 {
79 const Scalar rate_exact = Cf*exp(D*T);
80 const Scalar derive_exact = D * Cf * exp(D*T);
81 Antioch::KineticsConditions<Scalar> cond(T);
82
83 //KineticsConditions
84 Scalar rate = berthelot_rate(cond);
85 Scalar deriveRate = berthelot_rate.derivative(cond);
86
87 return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag;
88
89 berthelot_rate.rate_and_derivative(cond,rate,deriveRate);
90
91 return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag;
92
93 // T
94 rate = berthelot_rate(T);
95 deriveRate = berthelot_rate.derivative(T);
96
97 return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag;
98
99 berthelot_rate.rate_and_derivative(T,rate,deriveRate);
100
101 return_flag = check_rate_and_derivative(rate_exact,derive_exact,rate,deriveRate,T) || return_flag;
102
103 }
104 return return_flag;
105 }
106
107 template <typename Scalar>
tester()108 int tester()
109 {
110 Scalar Cf = 1.4;
111 Scalar D = -5.0;
112
113 Antioch::BerthelotRate<Scalar> berthelot_rate(Cf,D);
114
115 int return_flag = test_values(Cf,D,berthelot_rate);
116
117 Cf = 1e-7;
118 D = 1e-3;
119 berthelot_rate.set_Cf(Cf);
120 berthelot_rate.set_D(D);
121 return_flag = test_values(Cf,D,berthelot_rate) || return_flag;
122
123
124 Cf = 2.1e-11;
125 D = -1.2;
126 std::vector<Scalar> values(2);
127 values[0] = Cf;
128 values[1] = D;
129 berthelot_rate.reset_coefs(values);
130 return_flag = test_values(Cf,D,berthelot_rate) || return_flag;
131
132 return return_flag;
133 }
134
main()135 int main()
136 {
137 return (tester<double>() ||
138 tester<long double>() ||
139 tester<float>());
140 }
141