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 #include "antioch/read_reaction_set_data.h" 28 29 // Antioch 30 #include "antioch/read_reaction_set_data_instantiate_macro.h" 31 #include "antioch/string_utils.h" 32 #include "antioch/units.h" 33 #include "antioch/reaction_set.h" 34 #include "antioch/reaction_parsing.h" 35 #include "antioch/kinetics_parsing.h" 36 #include "antioch/physical_constants.h" 37 #include "antioch/ascii_parser.h" 38 #include "antioch/chemkin_parser.h" 39 #include "antioch/xml_parser.h" 40 41 namespace Antioch 42 { 43 template <typename NumericType> verify_unit_of_parameter(Units<NumericType> & default_unit,const std::string & provided_unit,const std::vector<std::string> & accepted_unit,const std::string & equation,const std::string & parameter_name)44 void verify_unit_of_parameter(Units<NumericType> & default_unit, 45 const std::string & provided_unit, 46 const std::vector<std::string> & accepted_unit, 47 const std::string & equation, 48 const std::string & parameter_name) 49 { 50 if(provided_unit.empty() && default_unit.is_united()) //no unit provided and there is one 51 { 52 antioch_unit_required(parameter_name,default_unit.get_symbol()); 53 }else // test 54 { 55 bool found(false); 56 for(unsigned int i = 0; i < accepted_unit.size(); i++) 57 { 58 if(Units<NumericType>(accepted_unit[i]).is_homogeneous(provided_unit)) 59 { 60 found = true; 61 break; 62 } 63 } 64 if(!found) 65 { 66 std::string errorstring("Error in reaction " + equation); 67 errorstring += "\n" + parameter_name + " parameter's unit should be homogeneous to \"" + accepted_unit[0] + "\""; 68 for(unsigned int i = 1; i < accepted_unit.size(); i++)errorstring += ", or \"" + accepted_unit[i] + "\""; 69 errorstring += " and you provided \"" + provided_unit + "\""; 70 antioch_unit_error(errorstring); 71 } 72 default_unit.set_unit(provided_unit); 73 } 74 } 75 76 template<typename NumericType> read_reaction_set_data(const std::string & filename,const bool verbose,ReactionSet<NumericType> & reaction_set,ParsingType type)77 void read_reaction_set_data( const std::string& filename, 78 const bool verbose, 79 ReactionSet<NumericType>& reaction_set, 80 ParsingType type ) 81 { 82 ParserBase<NumericType> * parser(NULL); 83 switch(type) 84 { 85 case ASCII: 86 parser = new ASCIIParser<NumericType>(filename,verbose); 87 break; 88 case CHEMKIN: 89 parser = new ChemKinParser<NumericType>(filename,verbose); 90 break; 91 case XML: 92 parser = new XMLParser<NumericType>(filename,verbose); 93 break; 94 default: 95 antioch_parsing_error("unknown type"); 96 } 97 98 //error or no reaction data 99 if(!parser->initialize())return; 100 101 const ChemicalMixture<NumericType>& chem_mixture = reaction_set.chemical_mixture(); 102 unsigned int n_species = chem_mixture.n_species(); 103 // Sanity Check on species 104 /* 105 tinyxml2::XMLElement* species = element->FirstChildElement("phase"); 106 species = species->FirstChildElement("speciesArray"); 107 108 std::vector<std::string> species_names; 109 110 std::cout << species->GetText() << std::endl; 111 112 SplitString(std::string(species->GetText()), 113 " ", 114 species_names, 115 false); 116 117 118 if( n_species != species_names.size() ) 119 { 120 std::cerr << "Error: Mismatch in n_species and the number of species specified in" << std::endl 121 << " the kinetics XML file " << filename << std::endl 122 << " Found species: " << species->GetText() << std::endl; 123 antioch_error(); 124 } 125 */ 126 127 std::map<std::string,KineticsModel::KineticsModel> kin_keyword; 128 kin_keyword["Constant"] = KineticsModel::CONSTANT; 129 kin_keyword["HercourtEssen"] = KineticsModel::HERCOURT_ESSEN; 130 kin_keyword["Berthelot"] = KineticsModel::BERTHELOT; 131 kin_keyword["Arrhenius"] = KineticsModel::ARRHENIUS; 132 kin_keyword["BerthelotHercourtEssen"] = KineticsModel::BHE; 133 kin_keyword["Kooij"] = KineticsModel::KOOIJ; 134 kin_keyword["ModifiedArrhenius"] = KineticsModel::KOOIJ; //for Arrhenius fans 135 kin_keyword["VantHoff"] = KineticsModel::VANTHOFF; 136 kin_keyword["photochemistry"] = KineticsModel::PHOTOCHEM; 137 138 std::map<KineticsModel::KineticsModel,unsigned int> kinetics_model_map; 139 kinetics_model_map[KineticsModel::CONSTANT] = 0; 140 kinetics_model_map[KineticsModel::HERCOURT_ESSEN] = 1; 141 kinetics_model_map[KineticsModel::BERTHELOT] = 2; 142 kinetics_model_map[KineticsModel::ARRHENIUS] = 3; 143 kinetics_model_map[KineticsModel::BHE] = 4; 144 kinetics_model_map[KineticsModel::KOOIJ] = 5; 145 kinetics_model_map[KineticsModel::VANTHOFF] = 7; 146 kinetics_model_map[KineticsModel::PHOTOCHEM] = 8; 147 148 std::vector<std::string> models; 149 models.push_back("Constant"); 150 models.push_back("HercourtEssen"); 151 models.push_back("Berthelot"); 152 models.push_back("Arrhenius"); 153 models.push_back("BerthelotHercourtEssen"); 154 models.push_back("Kooij"); 155 models.push_back("ModifiedArrhenius"); 156 models.push_back("VantHoff"); 157 models.push_back("photochemistry"); 158 159 std::map<std::string,ReactionType::ReactionType> proc_keyword; 160 proc_keyword["Elementary"] = ReactionType::ELEMENTARY; 161 proc_keyword["Duplicate"] = ReactionType::DUPLICATE; 162 proc_keyword["ThreeBody"] = ReactionType::THREE_BODY; 163 proc_keyword["threeBody"] = ReactionType::THREE_BODY; // Cantera/backward compatiblity 164 proc_keyword["LindemannFalloff"] = ReactionType::LINDEMANN_FALLOFF; 165 proc_keyword["TroeFalloff"] = ReactionType::TROE_FALLOFF; 166 proc_keyword["LindemannFalloffThreeBody"] = ReactionType::LINDEMANN_FALLOFF_THREE_BODY; 167 proc_keyword["TroeFalloffThreeBody"] = ReactionType::TROE_FALLOFF_THREE_BODY; 168 169 while (parser->reaction()) 170 { 171 if (verbose) std::cout << "Reaction \"" << parser->reaction_id() << "\":\n" 172 << " eqn: " << parser->reaction_equation() 173 << std::endl; 174 175 ReactionType::ReactionType typeReaction(ReactionType::ELEMENTARY); 176 KineticsModel::KineticsModel kineticsModel(KineticsModel::HERCOURT_ESSEN); // = 0 177 178 if (!parser->reaction_chemical_process().empty()) 179 { 180 if (verbose) std::cout << " type: " << parser->reaction_chemical_process() << std::endl; 181 if(!proc_keyword.count(parser->reaction_chemical_process())) 182 { 183 std::cerr << "The type of chemical process you provided (" << parser->reaction_chemical_process() << ")" 184 << " does not correspond to any Antioch knows.\n" 185 << "Implemented chemical processes are:\n" 186 << " Elementary (default)\n" 187 << " Duplicate\n" 188 << " ThreeBody\n" 189 << " LindemannFalloff\n" 190 << " TroeFalloff\n" 191 << " LindemannFalloffThreeBody\n" 192 << " TroeFalloffThreeBody\n" 193 << "See Antioch documentation for more details." 194 << std::endl; 195 antioch_not_implemented(); 196 } 197 typeReaction = proc_keyword[parser->reaction_chemical_process()]; 198 } 199 200 bool reversible(parser->reaction_reversible()); 201 if (verbose) std::cout << "reversible: " << reversible << std::endl; 202 203 kineticsModel = kin_keyword[parser->reaction_kinetics_model(models)]; 204 const std::string reading_kinetics_model = parser->reaction_kinetics_model(models); 205 206 // usually Kooij is called Arrhenius, check here 207 if(kineticsModel == KineticsModel::ARRHENIUS) 208 { 209 if(parser->verify_Kooij_in_place_of_Arrhenius()) 210 { 211 kineticsModel = KineticsModel::KOOIJ; 212 antioch_do_once( 213 std::cout << "In reaction(s) including " << parser->reaction_id() << "\n" 214 << "An equation of the form \"A * (T/Tref)^beta * exp(-Ea/(R*T))\" is a Kooij equation,\n" 215 << "I guess a modified Arrhenius could be a name too. Whatever, the correct label is\n" 216 << "\"Kooij\", or, << à la limite >> \"ModifiedArrhenius\". Please use those terms instead,\n" 217 << "thanks and a good day to you, user." << std::endl; 218 ); // antioch_do_once 219 } 220 } 221 222 // construct a Reaction object 223 Reaction<NumericType>* my_rxn = build_reaction<NumericType>(n_species, parser->reaction_equation(), 224 reversible,typeReaction,kineticsModel); 225 my_rxn->set_id(parser->reaction_id()); 226 227 // We will add the reaction, unless we do not have a 228 // reactant or product 229 bool relevant_reaction = true; 230 NumericType order_reaction(0); 231 std::vector<std::pair<std::string,int> > molecules_pairs; 232 233 if(parser->reactants_pairs(molecules_pairs)) 234 { 235 236 std::map<std::string,NumericType> orders = parser->reactants_orders(); 237 if (verbose) 238 { 239 std::cout << "\n reactants: "; 240 for(unsigned int ir = 0; ir < molecules_pairs.size(); ir++) 241 { 242 NumericType order = (orders.count(molecules_pairs[ir].first))?orders.at(molecules_pairs[ir].first):static_cast<NumericType>(molecules_pairs[ir].second); 243 std::cout << molecules_pairs[ir].first << ":" << molecules_pairs[ir].second << "," << order << ", "; 244 } 245 } 246 247 for( unsigned int p=0; p < molecules_pairs.size(); p++ ) 248 { 249 if(molecules_pairs[p].first == "e-") molecules_pairs[p].first = "e"; 250 251 if(verbose) std::cout << "\n " << molecules_pairs[p].first << " " << molecules_pairs[p].second; 252 253 if( !chem_mixture.species_name_map().count( molecules_pairs[p].first ) ) 254 { 255 relevant_reaction = false; 256 if (verbose) std::cout << "\n -> skipping this reaction (no reactant " << molecules_pairs[p].first << ")"; 257 } 258 else 259 { 260 NumericType order = (orders.count(molecules_pairs[p].first))?orders.at(molecules_pairs[p].first):static_cast<NumericType>(molecules_pairs[p].second); 261 my_rxn->add_reactant( molecules_pairs[p].first, 262 chem_mixture.species_name_map().find( molecules_pairs[p].first )->second, 263 molecules_pairs[p].second, order ); 264 order_reaction += order; 265 } 266 } 267 } 268 269 molecules_pairs.clear(); 270 if(parser->products_pairs(molecules_pairs)) 271 { 272 273 std::map<std::string,NumericType> orders = parser->products_orders(); 274 if(verbose) std::cout << "\n products: "; 275 for(unsigned int ir = 0; ir < molecules_pairs.size(); ir++) 276 { 277 NumericType order = (orders.count(molecules_pairs[ir].first))?orders.at(molecules_pairs[ir].first):static_cast<NumericType>(molecules_pairs[ir].second); 278 std::cout << molecules_pairs[ir].first << ":" << molecules_pairs[ir].second << "," << order << ", "; 279 } 280 281 for (unsigned int p=0; p < molecules_pairs.size(); p++) 282 { 283 if(molecules_pairs[p].first == "e-") molecules_pairs[p].first = "e"; 284 285 if(verbose) std::cout << "\n " << molecules_pairs[p].first << " " << molecules_pairs[p].second; 286 287 if( !chem_mixture.species_name_map().count( molecules_pairs[p].first ) ) 288 { 289 relevant_reaction = false; 290 if (verbose) std::cout << "\n -> skipping this reaction (no product " << molecules_pairs[p].first << ")"; 291 } 292 else 293 { 294 NumericType order = (orders.count(molecules_pairs[p].first))?orders.at(molecules_pairs[p].first):static_cast<NumericType>(molecules_pairs[p].second); 295 my_rxn->add_product( molecules_pairs[p].first, 296 chem_mixture.species_name_map().find( molecules_pairs[p].first )->second, 297 molecules_pairs[p].second, order ); 298 } 299 } 300 if(verbose) std::cout << std::endl; 301 } 302 303 if(!relevant_reaction) 304 { 305 if(verbose) std::cout << "skipped reaction\n\n"; 306 delete my_rxn; 307 continue; 308 } 309 310 while(parser->rate_constant(reading_kinetics_model)) //for duplicate and falloff models, several kinetics rate to load, no mixing allowed 311 { 312 313 /* Any data is formatted by the parser method. 314 * For any data required, parser sends back: 315 * - true/false if data is found 316 * - value of data 317 * - unit of data if found, empty string else 318 * - default unit of data 319 * 320 * The parser defines the defaults, note the special case 321 * of the pre-exponential parameters: 322 * its unit is [quantity-1]^(order - 1)/s, thus the parser 323 * defines only the [quantity-1] unit (SI unit is m^3/mol) 324 * as default. 325 */ 326 327 std::vector<NumericType> data; // for rate constant 328 329 Units<NumericType> def_unit; 330 int pow_unit(order_reaction - 1); 331 bool is_falloff(false); 332 bool is_k0(false); 333 //threebody always 334 if(my_rxn->type() == ReactionType::THREE_BODY)pow_unit++; 335 //falloff for k0 336 if(my_rxn->type() == ReactionType::LINDEMANN_FALLOFF || 337 my_rxn->type() == ReactionType::TROE_FALLOFF || 338 my_rxn->type() == ReactionType::LINDEMANN_FALLOFF_THREE_BODY || 339 my_rxn->type() == ReactionType::TROE_FALLOFF_THREE_BODY) 340 { 341 is_falloff = verbose; 342 //k0 is either determined by an explicit name, or is the first of unnamed rate constants 343 if(parser->is_k0(my_rxn->n_rate_constants(),reading_kinetics_model)) 344 { 345 pow_unit++; 346 is_k0 = verbose; 347 } 348 } 349 350 NumericType par_value(-1.); 351 std::vector<NumericType> par_values; 352 std::string par_unit; 353 std::string default_unit; 354 std::vector<std::string> accepted_unit; 355 356 // verbose as we read along 357 if(verbose)std::cout << " rate: " << models[kinetics_model_map[kineticsModel]] << " model\n"; 358 if(is_falloff){ 359 (is_k0)?std::cout << " Low pressure limit rate constant\n":std::cout << " High pressure limit rate constant\n"; 360 } 361 362 // pre-exponential, everyone 363 if(parser->rate_constant_preexponential_parameter(par_value, par_unit, default_unit)) 364 { 365 // using Units object to build accepted_unit 366 accepted_unit.clear(); 367 def_unit.set_unit("m3/mol"); 368 //to the m-1 power 369 if(pow_unit != 0) 370 { 371 def_unit *= pow_unit; 372 }else 373 { 374 def_unit.clear(); 375 } 376 def_unit.substract("s"); // per second 377 accepted_unit.push_back(def_unit.get_symbol()); 378 379 def_unit.set_unit(default_unit); 380 //to the m-1 power 381 if(pow_unit != 0) 382 { 383 def_unit *= pow_unit; 384 }else 385 { 386 def_unit.clear(); 387 } 388 def_unit.substract("s"); // per second 389 verify_unit_of_parameter(def_unit, par_unit, accepted_unit, my_rxn->equation(), "A"); 390 if(verbose) 391 { 392 std::cout << " A: " << par_value 393 << " " << def_unit.get_symbol() << std::endl; 394 } 395 data.push_back(par_value * def_unit.get_SI_factor()); 396 } 397 398 // beta, not everyone 399 if(( kineticsModel == KineticsModel::HERCOURT_ESSEN || 400 kineticsModel == KineticsModel::BHE || 401 kineticsModel == KineticsModel::KOOIJ || 402 kineticsModel == KineticsModel::VANTHOFF ) && 403 parser->rate_constant_power_parameter(par_value,par_unit,default_unit)) 404 { 405 accepted_unit.clear(); 406 accepted_unit.push_back(""); 407 def_unit.set_unit(default_unit); 408 verify_unit_of_parameter(def_unit, par_unit, accepted_unit, my_rxn->equation(), "beta"); 409 if(par_value == 0. && //if ARRHENIUS parameterized as KOOIJ, bad test, need to rethink it 410 // if a falloff, maybe the other reaction is a Kooij, we keep all the parameters, just in case 411 !(my_rxn->type() == ReactionType::LINDEMANN_FALLOFF || 412 my_rxn->type() == ReactionType::TROE_FALLOFF || 413 my_rxn->type() == ReactionType::LINDEMANN_FALLOFF_THREE_BODY || 414 my_rxn->type() == ReactionType::TROE_FALLOFF_THREE_BODY) 415 ) 416 { 417 418 std::cerr << "In reaction " << parser->reaction_id() << "\n" 419 << "An equation of the form \"A * exp(-Ea/(R*T))\" is an Arrhenius equation,\n" 420 << "and most certainly not a Kooij one\n" 421 << "it has been corrected, but please, change that in your file.\n" 422 << "Thanks and a good day to you, user." << std::endl; 423 kineticsModel = KineticsModel::ARRHENIUS; 424 }else 425 { 426 data.push_back(par_value * def_unit.get_SI_factor()); 427 } 428 if(verbose) 429 { 430 std::cout << " b: " << par_value << std::endl; 431 } 432 } 433 434 // activation energy, not everyone 435 if(( kineticsModel == KineticsModel::ARRHENIUS || 436 kineticsModel == KineticsModel::KOOIJ || 437 kineticsModel == KineticsModel::VANTHOFF ) && 438 parser->rate_constant_activation_energy_parameter(par_value,par_unit,default_unit)) 439 { 440 accepted_unit.clear(); 441 accepted_unit.push_back("J/mol"); 442 accepted_unit.push_back("K"); 443 def_unit.set_unit(default_unit); 444 verify_unit_of_parameter(def_unit, par_unit, accepted_unit, my_rxn->equation(), "Ea"); 445 data.push_back(par_value * def_unit.get_SI_factor()); 446 if(verbose) 447 { 448 std::cout << " E: " << par_value 449 << " " << def_unit.get_symbol() << std::endl; 450 } 451 } 452 453 454 // Berthelot coefficient (D), not everyone 455 if(( kineticsModel == KineticsModel::BERTHELOT || 456 kineticsModel == KineticsModel::BHE || 457 kineticsModel == KineticsModel::VANTHOFF ) && 458 parser->rate_constant_Berthelot_coefficient_parameter(par_value,par_unit,default_unit)) 459 { 460 accepted_unit.clear(); 461 accepted_unit.push_back("K"); 462 def_unit.set_unit(default_unit); 463 verify_unit_of_parameter(def_unit, par_unit, accepted_unit, my_rxn->equation(), "D"); 464 data.push_back(par_value * def_unit.get_SI_factor()); 465 if(verbose) 466 { 467 std::cout << " D: " << par_value 468 << " " << def_unit.get_symbol() << std::endl; 469 } 470 } 471 472 // Tref, not for everyone 473 if(kineticsModel == KineticsModel::HERCOURT_ESSEN || 474 kineticsModel == KineticsModel::BHE || 475 kineticsModel == KineticsModel::KOOIJ || 476 kineticsModel == KineticsModel::VANTHOFF) 477 { 478 par_value = 1.; 479 if(parser->rate_constant_Tref_parameter(par_value,par_unit,default_unit)) 480 { 481 accepted_unit.clear(); 482 accepted_unit.push_back("K"); 483 def_unit.set_unit(default_unit); 484 verify_unit_of_parameter(def_unit, par_unit, accepted_unit, my_rxn->equation(), "Tref"); 485 }else 486 { 487 antioch_parameter_required("Tref","1 K"); 488 } 489 data.push_back(par_value); 490 } 491 492 // scale E -> E/R 493 if(kineticsModel == KineticsModel::ARRHENIUS || 494 kineticsModel == KineticsModel::KOOIJ || 495 kineticsModel == KineticsModel::VANTHOFF) 496 { 497 parser->rate_constant_activation_energy_parameter(par_value,par_unit,default_unit); 498 (par_unit.empty())?def_unit.set_unit(default_unit):def_unit.set_unit(par_unit); 499 // now finding R unit: [Ea] / [K] 500 def_unit.substract("K"); 501 502 par_value = (def_unit.is_united())? 503 Antioch::Constants::R_universal<NumericType>() // Ea already tranformed in SI 504 :1.L; // no unit, so Ea already in K 505 data.push_back(par_value); 506 } 507 508 //photochemistry 509 // lambda is either a length (def nm) or cm-1 510 // cross-section has several possibilities if given 511 // * cm2 per bin: 512 // - length (typically nm) or cm-1 513 // * cm2 no bin given: 514 // - if given, lambda unit 515 // - if not, nm 516 517 // starting with lambda (for bin unit in cross-section) 518 // lambda is not in SI (m is really to violent), it will be nm 519 if(parser->rate_constant_lambda_parameter(par_values,par_unit,default_unit)) 520 { 521 antioch_assert_equal_to(kineticsModel,KineticsModel::PHOTOCHEM); 522 accepted_unit.clear(); 523 accepted_unit.push_back("nm"); 524 accepted_unit.push_back("cm-1"); 525 def_unit.set_unit(default_unit); 526 verify_unit_of_parameter(def_unit, par_unit, accepted_unit, my_rxn->equation(), "lambda"); 527 528 data.clear(); 529 530 if(def_unit.is_homogeneous("nm"))// it's a length 531 { 532 for(unsigned int il = 0; il < par_values.size(); il++) 533 { 534 data.push_back(par_values[il] * def_unit.factor_to_some_unit("nm")); 535 } 536 }else // it's the inverse of a length 537 { 538 for(unsigned int il = 0; il < par_values.size(); il++) 539 { 540 data.push_back(1.L/(par_values[il] * def_unit.factor_to_some_unit("nm-1"))); 541 } 542 } 543 544 //now the cross-section 545 546 NumericType bin_coefficient = (def_unit.is_homogeneous("nm"))?def_unit.factor_to_some_unit("nm"): 547 def_unit.factor_to_some_unit("nm-1"); 548 if(!parser->rate_constant_cross_section_parameter(par_values,par_unit,default_unit)) 549 { 550 std::cerr << "Where is the cross-section? In what universe have you photochemistry with a wavelength grid and no cross-section on it?" << std::endl; 551 antioch_error(); 552 } 553 //test length 554 if(par_values.size() != data.size()) 555 { 556 std::cerr << "Your cross-section vector and your lambda vector don't have the same size!\n" 557 << "What am I supposed to do with that?" 558 << std::endl; 559 antioch_error(); 560 } 561 562 /* here we will use two def unit: 563 * cs_unit, cm2 by default 564 * bin_unit, nm by default. 565 * 566 * strict rigorous unit is 567 * - (cs_unit - bin_unit): cm2/nm 568 * correct unit is 569 * - cs_unit: cm2 570 * 571 * so we need to test against those two possibilities. 572 * Now the funny part is that we test homogeneity, not 573 * equality, for generality purposes, so in case of strict 574 * rigorous unit, we need to decompose the read_unit into 575 * cross_section and bin units, so we can make the appropriate change. 576 * 577 * !TODO make the decomposition instead of strict equality 578 */ 579 580 accepted_unit.clear(); 581 accepted_unit.push_back("cm2"); // only cross-section 582 accepted_unit.push_back("cm2/nm"); // per bin, bin is length-like 583 accepted_unit.push_back("cm2/nm-1"); // per bin, bin is inverse length-like 584 verify_unit_of_parameter(def_unit, par_unit, accepted_unit, my_rxn->equation(), "cross-section"); 585 586 587 if(def_unit.is_homogeneous("cm2")) // gotta use the bin unit, real unit is [cs]/[provided bin unit] 588 { 589 for(unsigned int ics = 0; ics < par_values.size(); ics++) 590 { 591 data.push_back(par_values[ics] * def_unit.get_SI_factor() / bin_coefficient); //cs in SI, bin in nm or nm-1 592 } 593 }else // bin unit is provided 594 { 595 std::string target_unit = (def_unit.is_homogeneous("cm2/nm"))?"m2/nm":"m2/nm-1"; 596 for(unsigned int ics = 0; ics < par_values.size(); ics++) 597 { 598 data.push_back(par_values[ics] * def_unit.factor_to_some_unit(target_unit)); 599 } 600 } 601 602 } //end photochemistry 603 604 if(data.empty()) //replace the old "if no A parameters" as A is not required anymore 605 { 606 std::cerr << "Somehow, I have a bad feeling about a chemical reaction without any data parameters...\n" 607 << "This is too sad, I give up...\n" 608 << "Please, check the reaction " << my_rxn->equation() << " before coming back to me." << std::endl; 609 antioch_error(); //HEY!!! 610 } 611 612 KineticsType<NumericType>* rate = build_rate<NumericType>(data,kineticsModel); 613 614 my_rxn->add_forward_rate(rate); 615 616 } //end of duplicate/falloff kinetics description loop 617 618 // for falloff, we need a way to know which rate constant is the low pressure limit 619 // and which is the high pressure limit 620 // usually by calling the low pressure limite "k0". If nothing given, by default 621 // the first rate constant encountered is the low limit, 622 // so we need to change something only if the second rate constant has a "name" attribute 623 // of value "k0" 624 if(my_rxn->type() == ReactionType::LINDEMANN_FALLOFF || 625 my_rxn->type() == ReactionType::TROE_FALLOFF || 626 my_rxn->type() == ReactionType::LINDEMANN_FALLOFF_THREE_BODY || 627 my_rxn->type() == ReactionType::TROE_FALLOFF_THREE_BODY) 628 { 629 antioch_assert_equal_to(my_rxn->n_rate_constants(),2); 630 if(parser->where_is_k0(reading_kinetics_model) == 1) // second given is k0 631 { 632 my_rxn->swap_forward_rates(0,1); 633 } 634 } 635 636 637 std::vector<std::pair<std::string,NumericType> > efficiencies; 638 //efficiencies are only for three body reactions 639 if(parser->efficiencies(efficiencies)) 640 { 641 antioch_assert(ReactionType::THREE_BODY == my_rxn->type() || 642 ReactionType::LINDEMANN_FALLOFF_THREE_BODY == my_rxn->type() || 643 ReactionType::TROE_FALLOFF_THREE_BODY == my_rxn->type()); 644 645 for(unsigned int p = 0; p < efficiencies.size(); p++) 646 { 647 if(verbose)std::cout << "\n" << efficiencies[p].first << " " << efficiencies[p].second; 648 649 if(efficiencies[p].first == "e-") efficiencies[p].first = "e"; 650 651 // it is possible that the efficiency is specified for a species we are not 652 // modeling - so only add the efficiency if it is included in our list 653 if( chem_mixture.species_name_map().count( efficiencies[p].first ) ) 654 { 655 my_rxn->set_efficiency( efficiencies[p].first, 656 chem_mixture.species_name_map().find( efficiencies[p].first )->second, 657 efficiencies[p].second ); 658 } 659 } 660 if(verbose)std::cout << std::endl; 661 } 662 663 //F parameters only for Troe falloff 664 if(parser->Troe()) 665 { 666 antioch_assert(ReactionType::TROE_FALLOFF == my_rxn->type() || 667 ReactionType::TROE_FALLOFF_THREE_BODY == my_rxn->type()); 668 669 FalloffReaction<NumericType,TroeFalloff<NumericType> > *my_fall_rxn = 670 static_cast<FalloffReaction<NumericType,TroeFalloff<NumericType> > *> (my_rxn); 671 672 Units<NumericType> def_unit; 673 NumericType par_value(-1.); 674 std::string par_unit; 675 std::string default_unit; 676 std::vector<std::string> accepted_unit; 677 678 // alpha 679 if(!parser->Troe_alpha_parameter(par_value,par_unit,default_unit)) 680 { 681 std::cerr << "alpha parameter of Troe falloff missing!" << std::endl; 682 antioch_error(); 683 } 684 accepted_unit.clear(); 685 accepted_unit.push_back(""); 686 def_unit.set_unit(default_unit); 687 verify_unit_of_parameter(def_unit, par_unit, accepted_unit, my_fall_rxn->equation(), "alpha"); 688 my_fall_rxn->F().set_alpha(par_value * def_unit.get_SI_factor()); 689 if(verbose) 690 { 691 std::cout << " alpha: " << par_value 692 << " " << def_unit.get_symbol() << std::endl; 693 } 694 695 // T*** 696 if(!parser->Troe_T3_parameter(par_value,par_unit,default_unit)) 697 { 698 std::cerr << "T*** parameter of Troe falloff missing!" << std::endl; 699 antioch_error(); 700 } 701 accepted_unit.clear(); 702 accepted_unit.push_back("K"); 703 def_unit.set_unit(default_unit); 704 verify_unit_of_parameter(def_unit, par_unit, accepted_unit, my_fall_rxn->equation(), "T***"); 705 my_fall_rxn->F().set_T3(par_value * def_unit.get_SI_factor()); 706 if(verbose) 707 { 708 std::cout << " T***: " << par_value 709 << " " << def_unit.get_symbol() << std::endl; 710 } 711 712 // T* 713 if(!parser->Troe_T1_parameter(par_value,par_unit,default_unit)) 714 { 715 std::cerr << "T* parameter of Troe falloff missing!" << std::endl; 716 antioch_error(); 717 } 718 // accepted unit is the same 719 def_unit.set_unit(default_unit); 720 verify_unit_of_parameter(def_unit, par_unit, accepted_unit, my_fall_rxn->equation(), "T*"); 721 my_fall_rxn->F().set_T1(par_value * def_unit.get_SI_factor()); 722 if(verbose) 723 { 724 std::cout << " T*: " << par_value 725 << " " << def_unit.get_symbol() << std::endl; 726 } 727 728 // T** is optional 729 if(parser->Troe_T2_parameter(par_value,par_unit,default_unit)) 730 { 731 def_unit.set_unit(default_unit); 732 verify_unit_of_parameter(def_unit, par_unit, accepted_unit, my_fall_rxn->equation(), "T**"); 733 my_fall_rxn->F().set_T2(par_value * def_unit.get_SI_factor()); 734 if(verbose) 735 { 736 std::cout << " T**: " << par_value 737 << " " << def_unit.get_symbol() << std::endl; 738 } 739 } 740 } 741 742 reaction_set.add_reaction(my_rxn); 743 744 if(verbose) std::cout << "\n\n"; 745 } 746 747 if(parser) 748 delete parser; 749 } 750 751 // Instantiate 752 ANTIOCH_READ_REACTION_SET_DATA_INSTANTIATE(); 753 754 } // end namespace Antioch 755