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 #ifndef ANTIOCH_KINETICS_PARSING_H
28 #define ANTIOCH_KINETICS_PARSING_H
29 
30 //Antioch
31 #include "antioch/antioch_asserts.h"
32 #include "antioch/metaprogramming_decl.h"
33 #include "antioch/physical_constants.h"
34 #include "antioch/kinetics_type.h"
35 #include "antioch/constant_rate.h"
36 #include "antioch/hercourtessen_rate.h"
37 #include "antioch/berthelot_rate.h"
38 #include "antioch/arrhenius_rate.h"
39 #include "antioch/berthelothercourtessen_rate.h"
40 #include "antioch/kooij_rate.h"
41 #include "antioch/vanthoff_rate.h"
42 #include "antioch/photochemical_rate.h"
43 #include "antioch/reaction_enum.h"
44 
45 namespace Antioch
46 {
47 /*! cross-road to kinetics model
48  */
49 
50   template<typename CoeffType, typename VectorCoeffType>
51   KineticsType<CoeffType, VectorCoeffType>* build_rate( const VectorCoeffType &data,
52                                                         const KineticsModel::KineticsModel &kin );
53 
54   template<typename CoeffType, typename VectorCoeffType, typename VectorType>
55   void reset_rate( KineticsType<CoeffType,VectorCoeffType> & kin, const VectorType & coefs);
56 
57   /*!
58    * The rate constant models derived from the Arrhenius model have
59    * an activation energy which is stored and used, for efficiency reasons,
60    * in the reduced \f$\frac{E_a}{\mathrm{R}}\f$ form in K.
61    * This calls for a subtle management of the get/set function.
62    * This is done by a unit management, by default an activation energy is
63    * an energy by quantity, so the default is to consider the SI unit
64    * for an energy, (J/mol). Therefore if the provided value is in Kelvin,
65    * it should be explicitely provided.
66    */
67   template <typename CoeffType, typename VectorCoeffType>
68   void reset_parameter_of_rate(KineticsType<CoeffType,VectorCoeffType> & rate,
69                                KineticsModel::Parameters parameter,
70                                const CoeffType new_value, const std::string & unit = "SI");
71 
72   // vectorized parameter
73   template <typename CoeffType, typename VectorCoeffType>
74   void reset_parameter_of_rate(KineticsType<CoeffType,VectorCoeffType> & rate,
75                                KineticsModel::Parameters parameter,
76                                const CoeffType new_value, int l, const std::string & unit = "SI");
77 
78 
79 //----------------------------------------
80 
81   //! We take here the parameters as:
82   //  - \f$A\f$, \f$\beta\f$, \f$E_a\f$, \f$D\f$, \f$\mathrm{T_{ref}}\f$, \f$\mathrm{scale}\f$.
83   //  The \f$\mathrm{scale}\f$ parameter is the factor of \f$E_a\f$ from its unit
84   //  to K.
85   //
86   //  We check the number of parameters then initialize.
87   template<typename CoeffType, typename VectorCoeffType>
88   inline
build_rate(const VectorCoeffType & data,const KineticsModel::KineticsModel & kin)89   KineticsType<CoeffType, VectorCoeffType>* build_rate( const VectorCoeffType &data,
90                                                         const KineticsModel::KineticsModel &kin )
91   {
92     using std::pow;
93 
94     KineticsType<CoeffType,VectorCoeffType>* rate = NULL;
95 
96     switch(kin)
97       {
98       case(KineticsModel::CONSTANT):
99         {
100           antioch_assert_equal_to(1,data.size());
101           rate = new ConstantRate<CoeffType>(data[0]);//Cf
102         }
103         break;
104 
105       case(KineticsModel::HERCOURT_ESSEN):
106         {
107           antioch_assert_equal_to(3,data.size());
108           rate = new HercourtEssenRate<CoeffType>(data[0],data[1],data[2]);//Cf,eta,Tref
109         }
110         break;
111 
112       case(KineticsModel::BERTHELOT):
113         {
114           antioch_assert_equal_to(2,data.size());
115           rate = new BerthelotRate<CoeffType>(data[0],data[1]);// Cf, D
116         }
117         break;
118 
119       case(KineticsModel::ARRHENIUS):
120         {
121           antioch_assert_equal_to(3,data.size());
122           rate = new ArrheniusRate<CoeffType>(data[0],data[1],data[2]);//Cf,Ea,scale
123         }
124         break;
125 
126       case(KineticsModel::BHE):
127         {
128           antioch_assert_equal_to(4,data.size());
129           rate = new BerthelotHercourtEssenRate<CoeffType>(data[0],data[1],data[2],data[3]);//Cf,eta,D,Tref
130         }
131         break;
132 
133       case(KineticsModel::KOOIJ):
134         {
135           antioch_assert_equal_to(5,data.size());
136           rate = new KooijRate<CoeffType>(data[0],data[1],data[2],data[3],data[4]);//Cf,eta,Ea,Tref,scale
137         }
138         break;
139 
140       case(KineticsModel::VANTHOFF):
141         {
142           antioch_assert_equal_to(6,data.size());
143           rate = new VantHoffRate<CoeffType>(data[0],data[1],data[2],data[3],data[4],data[5]);//Cf,eta,Ea,D,Tref,scale
144         }
145         break;
146 
147       case(KineticsModel::PHOTOCHEM):
148         {
149           antioch_assert_equal_to(0,data.size()%2); //two vectors in one: lambda then cross-section
150           VectorCoeffType cs;
151           VectorCoeffType lambda;
152           for(unsigned int i = 0; i < data.size()/2; i++)
153           {
154              lambda.push_back(data[i]);
155              cs.push_back(data[i + data.size()/2]);
156           }
157           rate = new PhotochemicalRate<CoeffType,VectorCoeffType>(cs,lambda);//cs,lambda
158         }
159         break;
160 
161       default:
162         {
163           antioch_error();
164         }
165 
166       } // switch(kin)
167 
168     return rate;
169   }
170 
171   template<typename CoeffType, typename VectorCoeffType, typename VectorType>
reset_rate(KineticsType<CoeffType,VectorCoeffType> & kin,const VectorType & coefs)172   void reset_rate( KineticsType<CoeffType,VectorCoeffType> & kin, const VectorType & coefs)
173   {
174 
175     switch(kin.type())
176       {
177       case(KineticsModel::CONSTANT):
178         {
179           static_cast<ConstantRate<CoeffType>*>(&kin)->reset_coefs(coefs);
180         }
181         break;
182 
183       case(KineticsModel::HERCOURT_ESSEN):
184         {
185           static_cast< HercourtEssenRate<CoeffType>*>(&kin)->reset_coefs(coefs);
186         }
187         break;
188 
189       case(KineticsModel::BERTHELOT):
190         {
191           static_cast< BerthelotRate<CoeffType>*>(&kin)->reset_coefs(coefs);
192         }
193         break;
194 
195       case(KineticsModel::ARRHENIUS):
196         {
197           static_cast< ArrheniusRate<CoeffType>*>(&kin)->reset_coefs(coefs);
198         }
199         break;
200 
201       case(KineticsModel::BHE):
202         {
203           static_cast< BerthelotHercourtEssenRate<CoeffType>*>(&kin)->reset_coefs(coefs);
204         }
205         break;
206 
207       case(KineticsModel::KOOIJ):
208         {
209           static_cast< KooijRate<CoeffType>*>(&kin)->reset_coefs(coefs);
210         }
211         break;
212 
213       case(KineticsModel::VANTHOFF):
214         {
215           static_cast< VantHoffRate<CoeffType>*>(&kin)->reset_coefs(coefs);
216         }
217         break;
218 
219       case(KineticsModel::PHOTOCHEM):
220         {
221           static_cast<PhotochemicalRate<CoeffType,VectorCoeffType>*>(&kin)->reset_coefs(coefs);
222         }
223         break;
224 
225       default:
226         {
227           antioch_error();
228         }
229 
230       } // switch(kin.type())
231   }
232 
233 
234   template <typename CoeffType, typename VectorCoeffType>
reset_parameter_of_rate(KineticsType<CoeffType,VectorCoeffType> & rate,KineticsModel::Parameters parameter,const CoeffType new_value,const std::string & unit)235   void reset_parameter_of_rate(KineticsType<CoeffType,VectorCoeffType> & rate,
236                                KineticsModel::Parameters parameter,
237                                const CoeffType new_value, const std::string & unit)
238   {
239 
240 // this is crude at the moment, no test
241 // this will be replaced by a unit manager
242 // at some point, to be able to have an explicit
243 // custom internal unit system with appropriate testing
244     CoeffType new_coef = (unit == "SI")?new_value:
245                                         new_value * Units<typename value_type<CoeffType>::type>(unit).get_SI_factor();
246 
247 // Ea management, we want K, two possibilities now
248 // 1 - Ea is already in K
249 // 2 - Ea is in J.mol-1
250    if(parameter == KineticsModel::Parameters::E)
251    {
252       if(unit != "K")
253       {
254          new_coef = new_coef / Constants::R_universal<typename value_type<CoeffType>::type>();
255       }
256    }
257 
258 
259     switch(rate.type())
260       {
261       case(KineticsModel::CONSTANT):
262         {
263           static_cast<ConstantRate<CoeffType>*>(&rate)->set_parameter(parameter,new_coef);
264         }
265         break;
266 
267       case(KineticsModel::HERCOURT_ESSEN):
268         {
269           static_cast< HercourtEssenRate<CoeffType>*>(&rate)->set_parameter(parameter,new_coef);
270         }
271         break;
272 
273       case(KineticsModel::BERTHELOT):
274         {
275           static_cast< BerthelotRate<CoeffType>*>(&rate)->set_parameter(parameter,new_coef);
276         }
277         break;
278 
279       case(KineticsModel::ARRHENIUS):
280         {
281           static_cast< ArrheniusRate<CoeffType>*>(&rate)->set_parameter(parameter,new_coef);
282         }
283         break;
284 
285       case(KineticsModel::BHE):
286         {
287           static_cast< BerthelotHercourtEssenRate<CoeffType>*>(&rate)->set_parameter(parameter,new_coef);
288         }
289         break;
290 
291       case(KineticsModel::KOOIJ):
292         {
293           static_cast< KooijRate<CoeffType>*>(&rate)->set_parameter(parameter,new_coef);
294         }
295         break;
296 
297       case(KineticsModel::VANTHOFF):
298         {
299           static_cast< VantHoffRate<CoeffType>*>(&rate)->set_parameter(parameter,new_coef);
300         }
301         break;
302 
303       default:
304         {
305           antioch_error();
306         }
307 
308       } // switch(kin.type())
309   }
310 
311   template <typename CoeffType, typename VectorCoeffType>
reset_parameter_of_rate(KineticsType<CoeffType,VectorCoeffType> & rate,KineticsModel::Parameters parameter,const CoeffType new_value,int l,const std::string & unit)312   void reset_parameter_of_rate(KineticsType<CoeffType,VectorCoeffType> & rate,
313                                KineticsModel::Parameters parameter,
314                                const CoeffType new_value, int l, const std::string & unit)
315   {
316 
317 // this is crude at the moment, no test
318 // this will be replaced by a unit manager
319 // at some point, to be able to have an explicit
320 // custom internal unit system with appropriate testing
321     CoeffType new_coef = (unit == "SI")?new_value:
322                                         new_value * Units<typename value_type<CoeffType>::type>(unit).get_SI_factor();
323 
324     switch(rate.type())
325     {
326       case(KineticsModel::PHOTOCHEM):
327         {
328           static_cast<PhotochemicalRate<CoeffType,VectorCoeffType>*>(&rate)->set_parameter(parameter,l,new_coef);
329         }
330         break;
331 
332       default:
333         {
334           antioch_error();
335         }
336 
337       } // switch(kin.type())
338   }
339 
340 
341 
342 } // end namespace Antioch
343 
344 #endif // ANTIOCH_REACTION_PARSING_H
345