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