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 // Antioch
30 #include "antioch/vector_utils.h"
31
32 #include "antioch/reaction.h"
33 #include "antioch/elementary_reaction.h"
34
35 template <typename Scalar>
tester()36 int tester()
37 {
38 using std::abs;
39 using std::exp;
40 using std::pow;
41
42 const Scalar Cf = 1.4L;
43 const Scalar Ea = 5.0L;
44 const Scalar beta = 1.2L;
45 const Scalar D = 2.5e-2L;
46
47 const std::string equation("A + B -> C + D");
48 const unsigned int n_species(4);
49
50 int return_flag = 0;
51 std::vector<Scalar> mol_densities;
52 mol_densities.push_back(1e-2L);
53 mol_densities.push_back(1e-2L);
54 mol_densities.push_back(1e-2L);
55 mol_densities.push_back(1e-2L);
56
57 const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 100;
58
59 for(Scalar T = 300.1; T <= 2500.1; T += 10.)
60 {
61
62 const Antioch::KineticsConditions<Scalar> conditions(T);
63
64 for(unsigned int ikinmod = 0; ikinmod < 6; ikinmod++)
65 {
66
67 Antioch::KineticsType<Scalar> *rate_kinetics(NULL);
68 Scalar rate_exact;
69 Scalar derive_exact;
70 Antioch::KineticsModel::KineticsModel kin_mod;
71
72 switch(ikinmod)
73 {
74 case 0:
75 {
76 kin_mod = Antioch::KineticsModel::HERCOURT_ESSEN;
77 rate_kinetics = new Antioch::HercourtEssenRate<Scalar>(Cf,beta,1.);
78 rate_exact = Cf * pow(T,beta);
79 derive_exact = Cf * pow (T,beta) * beta/T;
80 break;
81 }
82 case 1:
83 {
84 kin_mod = Antioch::KineticsModel::BERTHELOT;
85 rate_kinetics = new Antioch::BerthelotRate<Scalar>(Cf,D);
86 rate_exact = Cf * exp(D*T);
87 derive_exact = Cf * exp(D*T) * D;
88 break;
89 }
90 case 2:
91 {
92 kin_mod = Antioch::KineticsModel::ARRHENIUS;
93 rate_kinetics = new Antioch::ArrheniusRate<Scalar>(Cf,Ea,1.);
94 rate_exact = Cf * exp(-Ea/T);
95 derive_exact = Cf * exp(-Ea/T) * Ea/pow(T,2);
96 break;
97 }
98 case 3:
99 {
100 kin_mod = Antioch::KineticsModel::BHE;
101 rate_kinetics = new Antioch::BerthelotHercourtEssenRate<Scalar>(Cf,beta,D,1.);
102 rate_exact = Cf * pow(T,beta) * exp(D * T);
103 derive_exact = Cf * pow(T,beta) * exp(D * T) * (beta/T + D);
104 break;
105 }
106 case 4:
107 {
108 kin_mod = Antioch::KineticsModel::KOOIJ;
109 rate_kinetics = new Antioch::KooijRate<Scalar>(Cf,beta,Ea,1.,1.);
110 rate_exact = Cf * pow(T,beta) * exp(-Ea/T);
111 derive_exact = Cf * pow(T,beta) * exp(-Ea/T) * (beta/T + Ea/pow(T,2));
112 break;
113 }
114 case 5:
115 {
116 kin_mod = Antioch::KineticsModel::VANTHOFF;
117 rate_kinetics = new Antioch::VantHoffRate<Scalar>(Cf,beta,Ea,D,1.,1.);
118 rate_exact = Cf * pow(T,beta) * exp(-Ea/T + D * T);
119 derive_exact = Cf * pow(T,beta) * exp(-Ea/T + D * T) * (D + beta/T + Ea/pow(T,2));
120 break;
121 }
122 }
123
124 Antioch::ElementaryReaction<Scalar> * elem_reaction = new Antioch::ElementaryReaction<Scalar>(n_species,equation,true,kin_mod);
125 elem_reaction->add_forward_rate(rate_kinetics);
126 Scalar rate1 = elem_reaction->compute_forward_rate_coefficient(mol_densities,conditions);
127 Scalar rate;
128 Scalar drate_dT;
129 std::vector<Scalar> drate_dx;
130 drate_dx.resize(n_species);
131 elem_reaction->compute_forward_rate_coefficient_and_derivatives(mol_densities,conditions,rate,drate_dT,drate_dx);
132
133 if( abs( (rate1 - rate_exact)/rate_exact ) > tol )
134 {
135 std::cout << std::scientific << std::setprecision(16)
136 << "Error: Mismatch in rate values." << std::endl
137 << "Kinetics model (see enum) " << kin_mod << std::endl
138 << "T = " << T << " K" << std::endl
139 << "rate(T) = " << rate1 << std::endl
140 << "rate_exact = " << rate_exact << std::endl;
141
142 return_flag = 1;
143 }
144 if( abs( (rate - rate_exact)/rate_exact ) > tol )
145 {
146 std::cout << std::scientific << std::setprecision(16)
147 << "Error: Mismatch in rate values." << std::endl
148 << "Kinetics model (see enum) " << kin_mod << std::endl
149 << "T = " << T << " K" << std::endl
150 << "rate(T) = " << rate << std::endl
151 << "rate_exact = " << rate_exact << std::endl;
152
153 return_flag = 1;
154 }
155 if( abs( (drate_dT - derive_exact)/derive_exact ) > tol )
156 {
157 std::cout << std::scientific << std::setprecision(16)
158 << "Error: Mismatch in rate derivative values." << std::endl
159 << "Kinetics model (see enum) " << kin_mod << std::endl
160 << "T = " << T << " K" << std::endl
161 << "drate_dT(T) = " << drate_dT << std::endl
162 << "derive_exact = " << derive_exact << std::endl;
163
164 return_flag = 1;
165 }
166 delete elem_reaction;
167 if(return_flag)return return_flag;
168 }
169 }
170
171 return return_flag;
172 }
173
main()174 int main()
175 {
176 return (tester<double>() ||
177 tester<long double>() ||
178 tester<float>());
179 }
180