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 // This class 28 #include "antioch/chemkin_parser.h" 29 30 // Antioch 31 #include "antioch/antioch_numeric_type_instantiate_macro.h" 32 #include "antioch/nasa_mixture.h" 33 #include "antioch/cea_curve_fit.h" 34 #include "antioch/nasa7_curve_fit.h" 35 36 // C++ 37 #include <sstream> 38 39 namespace Antioch 40 { 41 template <typename NumericType> ChemKinParser(const std::string & filename,bool verbose)42 ChemKinParser<NumericType>::ChemKinParser(const std::string &filename, bool verbose) 43 : ParserBase<NumericType>("ChemKin",filename,verbose,"!"), 44 _doc(filename.c_str()), 45 _duplicate_process(false), 46 _next_is_reverse(false) 47 { 48 if(!_doc.good()) 49 { 50 std::cerr << "ERROR: unable to load ChemKin file " << filename << std::endl; 51 antioch_error(); 52 } 53 54 if(this->verbose())std::cout << "Having opened file " << filename << std::endl; 55 56 _map[ParsingKey::SPECIES_SET] = "SPECIES"; 57 _map[ParsingKey::THERMO] = "THERMO"; 58 _map[ParsingKey::REACTION_DATA] = "REAC"; //REACTIONS || REAC 59 _map[ParsingKey::FALLOFF_LOW_NAME] = "LOW"; 60 _map[ParsingKey::TROE_FALLOFF] = "TROE"; 61 _map[ParsingKey::FORWARD_ORDER] = "FORD"; 62 _map[ParsingKey::BACKWARD_ORDER] = "RORD"; 63 64 // typically chemkin files list 65 // pre-exponential parameters in (cm3/mol)^(m-1)/s 66 // activation energy in cal/mol, but we want it in K. 67 // power parameter without unit 68 // if falloff, we need to know who's k0 and kinfty 69 // if photochemistry, we have a cross-section on a lambda grid 70 // cross-section typically in cm2/nm (cross-section on a resolution bin, 71 // if bin unit not given, it is lambda unit (supposed to anyway), and a warning message) 72 // lambda typically in nm, sometimes in ang, default considered here is nm 73 // you can also have cm-1, conversion is done with 74 // formulae nm = cm-1 * / * adapted factor 75 _default_unit[ParsingKey::PREEXP] = "cm3/mol"; 76 _default_unit[ParsingKey::POWER] = ""; 77 _default_unit[ParsingKey::ACTIVATION_ENERGY] = "cal/mol"; 78 _default_unit[ParsingKey::BERTHELOT_COEFFICIENT] = "K-1"; 79 _default_unit[ParsingKey::TREF] = "K"; 80 _default_unit[ParsingKey::HV_LAMBDA] = "nm"; 81 _default_unit[ParsingKey::HV_CROSS_SECTION] = "cm2/nm"; 82 _default_unit[ParsingKey::EFFICIENCY] = ""; 83 _default_unit[ParsingKey::TROE_F_ALPHA] = ""; 84 _default_unit[ParsingKey::TROE_F_TS] = "K"; 85 _default_unit[ParsingKey::TROE_F_TSS] = "K"; 86 _default_unit[ParsingKey::TROE_F_TSSS] = "K"; 87 88 _unit_custom_ea["CAL/MOL"] = "cal/mol"; 89 _unit_custom_ea["KCAL/MOL"] = "kcal/mol"; 90 _unit_custom_ea["JOULES/MOL"] = "J/mol"; 91 _unit_custom_ea["KELVINS"] = "K"; 92 _unit_custom_ea["EVOLTS"] = "eV"; 93 94 _unit_custom_A["MOLES"] = "cm3/mol"; 95 _unit_custom_A["MOLECULES"] = "cm3/molecule"; 96 } 97 98 template <typename NumericType> ~ChemKinParser()99 ChemKinParser<NumericType>::~ChemKinParser() 100 { 101 _doc.close(); 102 } 103 104 template <typename NumericType> change_file(const std::string & filename)105 void ChemKinParser<NumericType>::change_file(const std::string & filename) 106 { 107 _doc.close(); 108 _doc.open(filename.c_str()); 109 ParserBase<NumericType>::_file = filename; 110 if(!_doc.good()) 111 { 112 std::cerr << "ERROR: unable to load ChemKin file " << filename << std::endl; 113 antioch_error(); 114 } 115 116 if(this->verbose())std::cout << "Having opened file " << filename << std::endl; 117 } 118 119 template <typename NumericType> initialize()120 bool ChemKinParser<NumericType>::initialize() 121 { 122 std::string line; 123 ascii_getline(_doc,line); 124 bool init(false); 125 126 while(line.find(_map.at(ParsingKey::REACTION_DATA)) == std::string::npos) 127 { 128 if(!ascii_getline(_doc,line) || _doc.eof())break; 129 } 130 131 init = !(line.find(_map.at(ParsingKey::REACTION_DATA)) == std::string::npos); 132 133 if(init) 134 { 135 std::vector<std::string> keywords; 136 int nw = SplitString(line," ",keywords,false); 137 if(nw == 0)keywords.push_back(line); 138 for(unsigned int k = 1; k < keywords.size(); k++)// first word is REACTION 139 { 140 if(_unit_custom_ea.count(keywords[k])) 141 { 142 _default_unit[ParsingKey::PREEXP] = _unit_custom_ea.at(keywords[k]); 143 }else if(_unit_custom_A.count(keywords[k])) 144 { 145 _default_unit[ParsingKey::ACTIVATION_ENERGY] = _unit_custom_A.at(keywords[k]); 146 }else 147 { 148 antioch_parsing_error("ChemKin parser: I don't have ChemKin supporting this word as a unit:\n" + keywords[k]); 149 } 150 } 151 } 152 153 return init; 154 } 155 156 template <typename NumericType> species_list()157 const std::vector<std::string> ChemKinParser<NumericType>::species_list() 158 { 159 std::string line; 160 ascii_getline(_doc,line); 161 162 while(line.find(_map.at(ParsingKey::SPECIES_SET)) == std::string::npos) 163 { 164 if(!ascii_getline(_doc,line) || _doc.eof())break; 165 } 166 if(line.find(_map.at(ParsingKey::SPECIES_SET)) == std::string::npos) 167 { 168 antioch_parsing_error("ChemKin parser: haven't found SPECIES block"); 169 } 170 171 std::vector<std::string> species; 172 ascii_getline(_doc,line); 173 while(line.find(_spec.end_tag()) == std::string::npos) 174 { 175 std::vector<std::string> tmp; 176 SplitString(line," ",tmp,false); 177 if(tmp.empty()) 178 { 179 species.push_back(line); 180 }else 181 { 182 for(unsigned int i = 0; i < tmp.size(); i++) 183 { 184 species.push_back(tmp[i]); 185 } 186 } 187 ascii_getline(_doc,line); 188 } 189 190 return species; 191 } 192 193 template <typename NumericType> reaction()194 bool ChemKinParser<NumericType>::reaction() 195 { 196 197 /*default*/ 198 if(!_next_is_reverse) 199 { 200 _reversible = true; 201 _nrates = 0; 202 _crates = 0; 203 _duplicate_process = false; 204 205 _pow_unit = 0; 206 207 _kinetics_model = "Kooij"; //always Kooij, sometimes Arrhenius, should be corrected within the parser 208 _chemical_process = "Elementary"; 209 _equation.clear(); 210 211 _efficiencies.clear(); 212 213 _reactants.clear(); 214 _products.clear(); 215 _reactants_orders.clear(); 216 _products_orders.clear(); 217 }else 218 { 219 _next_is_reverse = false; 220 _reversible = false; 221 _nrates = 1; 222 _crates = 0; 223 _duplicate_process = false; 224 225 _pow_unit = 0; 226 227 _kinetics_model = "Kooij"; //always Kooij, sometimes Arrhenius, should be corrected within the parser 228 _chemical_process = "Elementary"; 229 230 231 _efficiencies.clear(); 232 233 // reverse products and reactants 234 std::vector<std::pair<std::string,NumericType> > tmp(_reactants); 235 std::vector<std::pair<std::string,NumericType> > tmp_o(_reactants_orders); 236 _reactants = _products; 237 _reactants_orders = _products_orders; 238 _products = tmp; 239 _products_orders = tmp_o; 240 241 // reverse the reaction 242 typename ChemKinDefinitions::Delim delim = _spec.equation_delimiter(_equation); 243 std::string str_delim = _spec.delim().at(delim); 244 std::size_t lim = _equation.find(str_delim); 245 std::string reac_to_prod = _equation.substr(0,lim); 246 std::string prod_to_reac = _equation.substr(lim + str_delim.size(), std::string::npos); 247 _equation = prod_to_reac + str_delim + reac_to_prod; 248 } 249 250 _A.clear(); 251 _b.clear(); 252 _Ea.clear(); 253 _D.clear(); 254 255 _Tref = 1.; 256 _Troe_alpha = -1.; 257 _Troe_T1 = -1.; 258 _Troe_T2 = -1.; 259 _Troe_T3 = -1.; 260 261 /* reaction */ 262 bool reac = false; 263 std::string line; 264 if(_cached_line.empty()) 265 { 266 reac = this->next_meaningful_line(line); 267 _cached_line = line; 268 }else 269 { 270 line = _cached_line; // already meaningful 271 reac = true; 272 } 273 274 //reaction found 275 while(!this->next_reaction(line)) 276 { 277 if(line.find(_spec.end_tag()) != std::string::npos || _doc.eof()) // equivalent 278 { 279 reac = false; 280 break; 281 } 282 // comment out 283 if(line.find(_spec.comment()) != std::string::npos)line.erase(line.find(_spec.comment()),std::string::npos); 284 // parsing 285 this->parse_a_line(line); 286 // getting on to the next line 287 this->next_meaningful_line(line); 288 } 289 _cached_line = line; 290 291 return reac; 292 } 293 294 template <typename NumericType> rate_constant(const std::string &)295 bool ChemKinParser<NumericType>::rate_constant(const std::string& /* kinetics_model */) 296 { 297 _crates++; 298 return (_crates <= _nrates); 299 } 300 301 template <typename NumericType> reactants_pairs(std::vector<std::pair<std::string,int>> & reactants_pair)302 bool ChemKinParser<NumericType>::reactants_pairs(std::vector<std::pair<std::string,int> >& reactants_pair) const 303 { 304 reactants_pair.clear(); 305 reactants_pair.resize(_reactants.size()); 306 307 for(unsigned int i = 0; i < _reactants.size(); i++) 308 { 309 reactants_pair[i] = std::make_pair(_reactants[i].first,(int)_reactants[i].second); 310 } 311 312 return !(_reactants.empty()); 313 } 314 315 template <typename NumericType> products_pairs(std::vector<std::pair<std::string,int>> & products_pair)316 bool ChemKinParser<NumericType>::products_pairs(std::vector<std::pair<std::string,int> > & products_pair) const 317 { 318 products_pair.clear(); 319 products_pair.resize(_products.size()); 320 for(unsigned int i = 0; i < _products.size(); i++) 321 { 322 products_pair[i] = std::make_pair(_products[i].first,(int)_products[i].second); 323 } 324 return !(_products.empty()); 325 } 326 327 template <typename NumericType> reactants_orders()328 const std::map<std::string,NumericType> ChemKinParser<NumericType>::reactants_orders() const 329 { 330 std::map<std::string,NumericType> orders; 331 for(unsigned int s = 0; s < _reactants_orders.size(); s++) 332 { 333 orders.insert(_reactants_orders[s]); 334 } 335 return orders; 336 } 337 338 template <typename NumericType> products_orders()339 const std::map<std::string,NumericType> ChemKinParser<NumericType>::products_orders() const 340 { 341 std::map<std::string,NumericType> orders; 342 for(unsigned int s = 0; s < _products_orders.size(); s++) 343 { 344 orders.insert(_products_orders[s]); 345 } 346 return orders; 347 } 348 349 template <typename NumericType> rate_constant_preexponential_parameter(NumericType & A,std::string & A_unit,std::string & def_unit)350 bool ChemKinParser<NumericType>::rate_constant_preexponential_parameter(NumericType & A, std::string & A_unit, std::string & def_unit) const 351 { 352 if(_crates <= _A.size()) 353 { 354 A = _A[_crates - 1]; 355 // there is no default unit (or always default), anyway, units are 356 // always explicit 357 def_unit = _default_unit.at(ParsingKey::PREEXP); 358 359 Units<NumericType> A_u(def_unit); 360 //to the m-1 power 361 int mult = (_chemical_process.find("Falloff") != std::string::npos && _crates == 1)?_pow_unit+1:_pow_unit; //falloff: +1 for k0 362 if(mult != 0) 363 { 364 A_u *= mult; 365 }else 366 { 367 A_u.clear(); 368 } 369 A_u.substract("s"); // per second 370 A_unit = A_u.get_symbol(); 371 } 372 return (_crates <= _A.size()); 373 } 374 375 template <typename NumericType> rate_constant_power_parameter(NumericType & b,std::string & b_unit,std::string & def_unit)376 bool ChemKinParser<NumericType>::rate_constant_power_parameter(NumericType & b, std::string & b_unit, std::string & def_unit) const 377 { 378 if(_crates <= _b.size()) 379 { 380 b = _b[_crates-1]; 381 // there is no default unit (or always default), anyway, units are 382 // always explicit 383 b_unit = _default_unit.at(ParsingKey::POWER); 384 def_unit = b_unit; 385 } 386 return (_crates <= _b.size()); 387 } 388 389 template <typename NumericType> rate_constant_activation_energy_parameter(NumericType & Ea,std::string & Ea_unit,std::string & def_unit)390 bool ChemKinParser<NumericType>::rate_constant_activation_energy_parameter(NumericType & Ea, std::string & Ea_unit, std::string & def_unit) const 391 { 392 if(_crates <= _Ea.size()) 393 { 394 Ea = _Ea[_crates-1]; 395 // there is no default unit (or always default), anyway, units are 396 // always explicit 397 Ea_unit = _default_unit.at(ParsingKey::ACTIVATION_ENERGY); 398 def_unit = Ea_unit; 399 } 400 return (_crates <= _Ea.size()); 401 } 402 403 template <typename NumericType> rate_constant_Tref_parameter(NumericType & Tref,std::string & Tref_unit,std::string & def_unit)404 bool ChemKinParser<NumericType>::rate_constant_Tref_parameter(NumericType & Tref, std::string & Tref_unit, std::string & def_unit) const 405 { 406 Tref = 1.L; 407 Tref_unit = _default_unit.at(ParsingKey::TREF); 408 def_unit = Tref_unit; 409 return (_crates <= _b.size()); 410 } 411 412 template <typename NumericType> efficiencies(std::vector<std::pair<std::string,NumericType>> & par_values)413 bool ChemKinParser<NumericType>::efficiencies(std::vector<std::pair<std::string,NumericType> > & par_values) const 414 { 415 par_values = _efficiencies; 416 return !_efficiencies.empty(); 417 } 418 419 template <typename NumericType> Troe_alpha_parameter(NumericType & alpha,std::string & alpha_unit,std::string & def_unit)420 bool ChemKinParser<NumericType>::Troe_alpha_parameter(NumericType & alpha, std::string & alpha_unit, std::string & def_unit) const 421 { 422 alpha = _Troe_alpha; 423 alpha_unit = _default_unit.at(ParsingKey::TROE_F_ALPHA); 424 def_unit = alpha_unit; 425 426 return this->Troe(); 427 } 428 429 template <typename NumericType> Troe_T1_parameter(NumericType & T1,std::string & T1_unit,std::string & def_unit)430 bool ChemKinParser<NumericType>::Troe_T1_parameter(NumericType & T1, std::string & T1_unit, std::string & def_unit) const 431 { 432 T1 = _Troe_T1; 433 T1_unit = _default_unit.at(ParsingKey::TROE_F_TS); 434 def_unit = T1_unit; 435 436 return this->Troe(); 437 } 438 439 template <typename NumericType> Troe_T2_parameter(NumericType & T2,std::string & T2_unit,std::string & def_unit)440 bool ChemKinParser<NumericType>::Troe_T2_parameter(NumericType & T2, std::string & T2_unit, std::string & def_unit) const 441 { 442 T2 = _Troe_T2; 443 T2_unit = _default_unit.at(ParsingKey::TROE_F_TSS); 444 def_unit = T2_unit; 445 446 return (_Troe_T2 > 0.); 447 } 448 449 template <typename NumericType> Troe_T3_parameter(NumericType & T3,std::string & T3_unit,std::string & def_unit)450 bool ChemKinParser<NumericType>::Troe_T3_parameter(NumericType & T3, std::string & T3_unit, std::string & def_unit) const 451 { 452 T3 = _Troe_T3; 453 T3_unit = _default_unit.at(ParsingKey::TROE_F_TSSS); 454 def_unit = T3_unit; 455 456 return this->Troe(); 457 } 458 459 template <typename NumericType> parse_a_line(const std::string & line)460 void ChemKinParser<NumericType>::parse_a_line(const std::string & line) 461 { 462 std::string capital_line(line); 463 std::transform(capital_line.begin(),capital_line.end(), capital_line.begin(),::toupper); 464 465 // reaction equation 466 if(line.find(_spec.delim().at(_spec.REVERSIBLE)) != std::string::npos) //equation a beta ea, "=" 467 { 468 this->parse_equation_coef(line); 469 _nrates++; 470 471 // reverse reaction 472 }else if(capital_line.find(_spec.reversible()) != std::string::npos) // reversible reaction is explicit "REV" 473 { 474 _reversible = false; 475 if(_next_is_reverse) 476 { 477 this->parse_reversible_parameters(line); 478 _next_is_reverse = false; 479 }else 480 { 481 _next_is_reverse = true; 482 } 483 484 // duplicate reaction 485 }else if(capital_line.find(_spec.duplicate()) != std::string::npos) // duplicate, "DUPLICATE" or "DUP" , search for "DUP" 486 { 487 _chemical_process = "Duplicate"; 488 _duplicate_process = true; 489 490 // custom forward orders 491 }else if(capital_line.find(_map.at(ParsingKey::FORWARD_ORDER)) != std::string::npos) // forward order, 492 { 493 this->parse_forward_orders(line); 494 // custom backward orders 495 }else if(capital_line.find(_map.at(ParsingKey::BACKWARD_ORDER)) != std::string::npos) // backward order, 496 { 497 this->parse_backward_orders(line); 498 // data about pressure dependence 499 }else if(line.find(_spec.parser()) != std::string::npos) // "/" 500 { 501 if(_chemical_process.find("Falloff") != std::string::npos && 502 _chemical_process.find("ThreeBody") == std::string::npos)_chemical_process += "ThreeBody"; 503 this->parse_coefficients_line(line); 504 }else 505 { 506 antioch_parsing_error("ChemKin parser: Can't parse this line:\n" + line); 507 } 508 } 509 510 template <typename NumericType> parse_equation_coef(const std::string & line)511 void ChemKinParser<NumericType>::parse_equation_coef(const std::string & line) 512 { 513 514 std::vector<std::string> out; 515 SplitString(line," ",out,false); 516 if(out.size() < 4)antioch_parsing_error("ChemKin parser: unrecognized reaction input line:\n" + line); 517 518 // parameters are the three last components 519 _Ea.push_back(std::atof(out[out.size()-1].c_str())); // last 520 _b.push_back( std::atof(out[out.size()-2].c_str())); // last - 1 521 _A.push_back( std::atof(out[out.size()-3].c_str())); // last - 2 522 523 // equation is the rest, can be 1 word as nreac + nprod + 1 524 // so making it 1 whatever happens 525 std::string equation; 526 for(unsigned int i = 0; i < out.size() - 3; i++) 527 { 528 equation += out[i]; 529 } 530 // in case of duplicate, the equation will be 531 // printed several times, we care only for the 532 // first 533 if(_reactants.empty())this->parse_equation(equation); 534 } 535 536 template <typename NumericType> parse_reversible_parameters(const std::string & line)537 void ChemKinParser<NumericType>::parse_reversible_parameters(const std::string & line) 538 { 539 // line looks like "REV / A beta Ea /" 540 541 std::vector<std::string> out; 542 SplitString(line,_spec.parser(),out,false); 543 if(out.size() < 2)antioch_parsing_error("ChemKin parser: unrecognized reversible reaction parameters input line:\n" + line); 544 545 std::vector<std::string> pars; 546 SplitString(out[1]," ",pars,false); 547 if(pars.size() < 3)antioch_parsing_error("ChemKin parser: unrecognized reversible reaction parameters input line:\n" + line); 548 549 // parameters are given 550 _A.push_back( std::atof(pars[0].c_str())); 551 _b.push_back( std::atof(pars[1].c_str())); 552 _Ea.push_back(std::atof(pars[2].c_str())); 553 } 554 555 template <typename NumericType> parse_equation(std::string & equation)556 void ChemKinParser<NumericType>::parse_equation(std::string &equation) 557 { 558 _equation = equation; 559 560 // first we need to treat the (+M) case (falloff) 561 // as it is not compatible with SplitString using ChemKinDefinitions::PLUS (+) 562 // as delimiter 563 if(equation.find(_spec.symbol().at(ChemKinDefinitions::FAL)) != std::string::npos) 564 { 565 // Lindemann by default 566 _chemical_process = "LindemannFalloff"; 567 // supress all occurences 568 while(equation.find(_spec.symbol().at(ChemKinDefinitions::FAL)) != std::string::npos) 569 { 570 equation.erase(equation.find(_spec.symbol().at(ChemKinDefinitions::FAL)),4); 571 } 572 } 573 574 std::vector<std::string> out; 575 576 // now we're singling the equation separator (=, =>, <=>) 577 typename ChemKinDefinitions::Delim delim(_spec.equation_delimiter(equation)); 578 if(delim == ChemKinDefinitions::ERROR)antioch_parsing_error("ChemKin parser: badly written equation:\n" + equation); 579 580 //between ChemKinDefinitions::PLUS 581 equation.insert(equation.find(_spec.delim().at(delim)),_spec.delim().at(ChemKinDefinitions::PLUS)); 582 equation.insert(equation.find(_spec.delim().at(delim)) + _spec.delim().at(delim).size(),_spec.delim().at(ChemKinDefinitions::PLUS)); 583 584 585 SplitString(equation,_spec.delim().at(ChemKinDefinitions::PLUS),out,true); //empties are cations charge, formatted as reac_i equation_separator prod_i 586 if(out.size() < 3)antioch_parsing_error("ChemKin parser: unrecognized reaction equation:\n" + equation); 587 /* cases are: 588 - equation_separator => switch between reac and prod 589 - molecules: 590 * reactant or prod, cation if followed by empty 591 * M if three body 592 */ 593 bool prod(false); 594 for(unsigned int i = 0; i < out.size(); i++) 595 { 596 if(out[i] == _spec.delim().at(delim)) 597 { 598 _reversible = !(delim == ChemKinDefinitions::IRREVERSIBLE); 599 prod = true; 600 // is it a third-body reaction? 601 }else if(out[i] == _spec.symbol().at(ChemKinDefinitions::TB)) 602 { 603 if(_chemical_process.find("Falloff") != std::string::npos) 604 std::cerr << "WARNING: ChemKin parser: it seems you want both a falloff and a three-body reaction in your equation:\n" << equation << std::endl; 605 606 _chemical_process = "ThreeBody"; 607 608 }else // here's a regular molecule 609 { 610 611 // no assumption on the charge, you can put as many '+' as you want 612 // adding empties, then skipping them 613 unsigned int j(1); 614 if(i+j < out.size()) // don't go beyond the last one 615 { 616 while(out[i+j].empty()) 617 { 618 out[i] += "+"; 619 j++; 620 } 621 } 622 j--; // back to non empty 623 (prod)?_products.push_back(this->parse_molecule(out[i])) 624 : 625 _reactants.push_back(this->parse_molecule(out[i])); 626 i += j; // skip empties 627 } 628 } 629 630 // checking for real stoichiometric coeffs 631 bool real(false); 632 for(unsigned int i = 0; i < _reactants.size(); i++) 633 { 634 if(this->after_coma_digits(_reactants[i].second)) 635 { 636 real = true; 637 break; 638 } 639 } 640 if(!real) 641 { 642 for(unsigned int i = 0; i < _products.size(); i++) 643 { 644 if(this->after_coma_digits(_products[i].second)) 645 { 646 real = true; 647 break; 648 } 649 } 650 } 651 652 if(real)this->rescale_stoichiometry(); 653 654 // orders are by default the stoichiometric coeffs 655 _reactants_orders = _reactants; 656 _products_orders = _products; 657 658 // order of reaction 659 for(unsigned int r = 0; r < _reactants.size(); r++) 660 { 661 _pow_unit += (unsigned int)_reactants[r].second; 662 } 663 _pow_unit--; 664 if(_chemical_process == "ThreeBody")_pow_unit++; 665 } 666 667 template <typename NumericType> parse_molecule(const std::string & molecule)668 std::pair<std::string,NumericType> ChemKinParser<NumericType>::parse_molecule(const std::string & molecule) 669 { 670 // as long as we have a number ({[0-9],.}), it is the stoichiometric coefficient 671 unsigned int pos(0); 672 while(this->is_real_number(molecule[pos]))pos++; 673 674 NumericType stoi = (pos == 0)?1.:std::atof(molecule.substr(0,pos+1).c_str()); 675 676 return std::make_pair(molecule.substr(pos,std::string::npos),stoi); 677 } 678 679 template <typename NumericType> is_real_number(const char & c)680 bool ChemKinParser<NumericType>::is_real_number(const char & c) const 681 { 682 return (c == '0' || c == '1' || c == '2' || 683 c == '3' || c == '4' || c == '5' || 684 c == '6' || c == '7' || c == '8' || 685 c == '9' || c == '.'); 686 687 } 688 689 template <typename NumericType> parse_orders(const std::string & line,std::vector<std::pair<std::string,NumericType>> & reaction_orders)690 void ChemKinParser<NumericType>::parse_orders(const std::string & line, std::vector<std::pair<std::string, NumericType> > & reaction_orders) 691 { 692 // 1 - parse the line 693 std::vector<std::string> out; 694 SplitString(line,_spec.parser(),out,false); 695 out.erase(out.begin()); // keyword (RORD or FORD) 696 697 std::vector< std::pair<std::string, NumericType> > orders; 698 for(unsigned int i = 0; i < out.size(); i++) 699 { 700 std::vector<std::string> you; 701 SplitString(out[i]," ",you,false); 702 if(you.size() != 2)antioch_parsing_error("ChemKin parser: I don't recognize this part:\n" + out[i]); 703 NumericType order = std::atof(you[1].c_str()); 704 orders.push_back(std::make_pair(you[0],order)); 705 } 706 707 // 2 - replace if found 708 for(unsigned int s = 0; s < reaction_orders.size(); s++) 709 { 710 // good ol' search & find without any clue about where it can be 711 unsigned int i; 712 for(i = 0; i < orders.size(); i++) 713 { 714 if(orders[i].first == reaction_orders[s].first) 715 { 716 reaction_orders[s] = orders[i]; 717 break; 718 } 719 } 720 } 721 } 722 723 template <typename NumericType> parse_forward_orders(const std::string & line)724 void ChemKinParser<NumericType>::parse_forward_orders(const std::string & line) 725 { 726 this->parse_orders(line, _reactants_orders); 727 } 728 729 template <typename NumericType> parse_backward_orders(const std::string & line)730 void ChemKinParser<NumericType>::parse_backward_orders(const std::string & line) 731 { 732 this->parse_orders(line, _products_orders); 733 } 734 735 template <typename NumericType> parse_coefficients_line(const std::string & line)736 void ChemKinParser<NumericType>::parse_coefficients_line(const std::string &line) 737 { 738 std::vector<std::string> out; 739 SplitString(line,_spec.parser(),out,false); 740 if(out.size() < 2)antioch_parsing_error("ChemKin parser: can't parse this line:\n" + line); 741 // can be LOW, TROE or coefficients, anything else is ignored 742 //accounts for blank spaces 743 if(out.front().find(_map.at(ParsingKey::TROE_FALLOFF)) != std::string::npos) //TROE, alpha, T***, T*, T** 744 { 745 antioch_assert_greater_equal(out.size(),2); 746 std::vector<std::string> troe_par; 747 SplitString(out[1]," ",troe_par,false); 748 if(troe_par.size() < 3)antioch_parsing_error("ChemKin parser: Troe parameters error while reading:\n" + line); 749 750 _Troe_alpha = std::atof(troe_par[0].c_str()); 751 _Troe_T3 = std::atof(troe_par[1].c_str()); 752 _Troe_T1 = std::atof(troe_par[2].c_str()); 753 _Troe_T2 = (troe_par.size() == 4)?std::atof(troe_par[3].c_str()):-1.L; 754 755 _chemical_process = "TroeFalloff"; 756 } 757 else if(out.front().find(_map.at(ParsingKey::FALLOFF_LOW_NAME)) != std::string::npos) // k0 758 { 759 antioch_assert_greater_equal(out.size(),2); 760 std::vector<std::string> k0_par; 761 SplitString(out[1]," ",k0_par,false); 762 if(k0_par.size() != 3)antioch_parsing_error("ChemKin parser: Falloff k0 parameters error while reading:\n" + line); 763 764 _A.insert( _A.begin(),std::atof( k0_par[0].c_str())); 765 _b.insert( _b.begin(),std::atof( k0_par[1].c_str())); 766 _Ea.insert(_Ea.begin(),std::atof(k0_par[2].c_str())); 767 768 _nrates++; 769 770 }else if(_chemical_process.find("ThreeBody") != std::string::npos)// efficiencies 771 // in case it is superfluous (or we need to redesign pressure-dependance) 772 { 773 for(unsigned int i = 0; i < out.size(); i++) 774 { 775 // Trim from the left 776 out[i].erase(out[i].begin(), std::find_if(out[i].begin(), out[i].end(), std::not1(std::ptr_fun<int, int>(std::isspace)))); 777 // Trim from the right 778 out[i].erase(std::find_if(out[i].rbegin(), out[i].rend(), std::not1(std::ptr_fun<int, int>(std::isspace))).base(), out[i].end()); 779 780 if(out[i].empty()) 781 { 782 out.erase(out.begin() + i); 783 i--; 784 } 785 } 786 if(out.size()%2 != 0)antioch_parsing_error("ChemKin parser: efficiencies parsing failed:\n" + line); 787 for(unsigned int c = 0; c < out.size(); c += 2) 788 { 789 _efficiencies.push_back(std::make_pair( out[c], std::atof(out[c+1].c_str()))); // mol / coef 790 } 791 } 792 793 } 794 795 template <typename NumericType> rescale_stoichiometry()796 void ChemKinParser<NumericType>::rescale_stoichiometry() 797 { 798 // we suppose rational number r = p / q with p and q integers, 799 // we look for q kinda stupidly (if this is more than 2 or 3, it 800 // is really ridiculous) 801 std::vector<unsigned int> not_int_factor; 802 for(unsigned int i = 0; i < _reactants.size(); i++) 803 { 804 if(this->after_coma_digits(_reactants[i].second)) 805 { 806 unsigned int fac = this->factor_to_int(_reactants[i].second); 807 unsigned int i(0); 808 for(i = 0; i < not_int_factor.size(); i++) 809 { 810 if(fac == not_int_factor[i])break; 811 } 812 if(i < not_int_factor.size() - 1)not_int_factor.push_back(fac); 813 } 814 } 815 for(unsigned int i = 0; i < _products.size(); i++) 816 { 817 if(this->after_coma_digits(_products[i].second)) 818 { 819 unsigned int fac = this->factor_to_int(_products[i].second); 820 unsigned int i(0); 821 for(i = 0; i < not_int_factor.size(); i++) 822 { 823 if(fac == not_int_factor[i])break; 824 } 825 if(i < not_int_factor.size() - 1)not_int_factor.push_back(fac); 826 } 827 } 828 829 // global factor will be brutal multiplication 830 // don't want to decompose into prime numbers 831 832 unsigned int factor(1); 833 for(unsigned int i = 0; i < not_int_factor.size(); i++) 834 { 835 factor *= not_int_factor[i]; 836 } 837 838 // rescale everyone 839 for(unsigned int i = 0; i < _reactants.size(); i++) 840 { 841 _reactants[i].second *= (NumericType)factor; 842 } 843 for(unsigned int i = 0; i < _products.size(); i++) 844 { 845 _products[i].second *= (NumericType)factor; 846 } 847 } 848 849 template <typename NumericType> factor_to_int(NumericType number)850 unsigned int ChemKinParser<NumericType>::factor_to_int(NumericType number) const 851 { 852 // now find p and q associated with this number 853 unsigned int down(2); 854 const unsigned int limit(150); 855 while(this->after_coma_digits(number * (NumericType) down)) 856 { 857 down++; 858 if(down > limit) 859 { 860 std::stringstream os; 861 os << "real is " << number << " and multiplicative factor limit is " << limit; 862 antioch_parsing_error("ChemKin parser: could not find integer from real\n:" + os.str()); 863 } 864 } 865 866 return down; 867 } 868 869 template <typename NumericType> after_coma_digits(NumericType number)870 bool ChemKinParser<NumericType>::after_coma_digits(NumericType number) const 871 { 872 // smallest number tolerated 1e-3, it is already ridiculous 873 const NumericType eps(1e-3); 874 return (std::abs(number - std::floor(number)) > eps); 875 } 876 877 template <typename NumericType> next_reaction(const std::string & input_line)878 bool ChemKinParser<NumericType>::next_reaction(const std::string & input_line) 879 { 880 bool out(false); 881 if(input_line.find(_spec.delim().at(ChemKinDefinitions::REVERSIBLE)) != std::string::npos || 882 input_line.find(_spec.end_tag()) != std::string::npos) out = true; 883 884 if(_next_is_reverse) out = true; // reversible given, get out 885 886 if(input_line == _cached_line || _cached_line.empty()) out = false; // current reaction 887 888 if(_duplicate_process && input_line.find(_spec.delim().at(ChemKinDefinitions::REVERSIBLE)) != std::string::npos)// if we find a reaction and it is the same than the cached one 889 { 890 out = false; // we suppose in duplicate reaction 891 std::vector<std::string> inputs; 892 SplitString(input_line," ",inputs,false); 893 if(inputs.size() < 4)antioch_parsing_error("ChemKin parser: unrecognized reaction input line:\n" + input_line); 894 // getting the equation 895 std::string test_eq; 896 for(unsigned int i = 0; i < inputs.size() - 3; i++) 897 { 898 test_eq += inputs[i]; 899 } 900 901 // the reactants and products should appear in the equation 902 std::vector<unsigned int> index_reac(_reactants.size(),0); 903 for(unsigned int ir = 0; ir < _reactants.size(); ir++) 904 { 905 if(input_line.find(_reactants[ir].first) == std::string::npos) // not the same reaction 906 { 907 out = true; 908 break; 909 } 910 index_reac.push_back(input_line.find(_reactants[ir].first)); 911 } 912 913 if(!out) 914 { 915 for(unsigned int ip = 0; ip < _products.size(); ip++) 916 { 917 if(input_line.find(_products[ip].first) == std::string::npos) // not the same reaction 918 { 919 out = true; 920 break; 921 } 922 for(unsigned int i = 0; i < index_reac.size(); i++) // products appear after reactants 923 { 924 // if the product appears before the reactant 925 if(input_line.find(_products[ip].first) < index_reac[i]) // not the same reaction 926 { 927 out = true; 928 // now in case someone wants the same molecule as reactant and product 929 for(unsigned int jr = 0; jr < _reactants.size(); jr++) 930 { 931 // find because H2O triggers H2O2... 932 if(_reactants[jr].first.find(_products[ip].first) != std::string::npos)out = false; 933 } 934 break; 935 } 936 } 937 if(out)break; 938 } 939 } 940 } 941 942 return out; 943 } 944 945 template <typename NumericType> next_meaningful_line(std::string & line)946 bool ChemKinParser<NumericType>::next_meaningful_line(std::string & line) 947 { 948 if(_next_is_reverse)return false; 949 bool out(true); 950 // skip empty lines 951 ascii_getline(_doc,line); 952 while(line.empty() || _spec.is_comment(line[0])) // fully commented alone line 953 { 954 if(!ascii_getline(_doc,line) || // getline 955 line.find(_spec.end_tag()) != std::string::npos || _doc.eof() // end of file 956 ) 957 { 958 out = false; 959 break; 960 } 961 } 962 963 return out; 964 } 965 966 template <typename NumericType> 967 template <typename CurveType> read_thermodynamic_data_root(NASAThermoMixture<NumericType,CurveType> & thermo)968 void ChemKinParser<NumericType>::read_thermodynamic_data_root(NASAThermoMixture<NumericType, CurveType >& thermo) 969 { 970 971 // finding thermo 972 std::string line; 973 ascii_getline(_doc,line); 974 975 while(line.find(_map.at(ParsingKey::THERMO)) == std::string::npos) 976 { 977 if(!ascii_getline(_doc,line) || _doc.eof())break; 978 } 979 980 if(!_doc.good()) 981 { 982 std::cerr << "Thermodynamics description not found" << std::endl; 983 antioch_error(); 984 } 985 986 const ChemicalMixture<NumericType>& chem_mixture = thermo.chemical_mixture(); 987 988 std::string name; 989 std::vector<NumericType> coeffs; 990 std::vector<NumericType> temps(3,0.); 991 992 this->skip_comments(_doc); 993 994 // chemkin classic 995 // only two intervals 996 // \todo: the parser should allow custom 997 // intervals definition 998 while (_doc.good()) 999 { 1000 std::stringstream tmp; 1001 this->skip_comments(_doc); // comments in middle 1002 1003 if(!ascii_getline(_doc,line))break; 1004 1005 if(line.find(_spec.end_tag()) != std::string::npos)break; //END 1006 1007 // question is: will we have temps or NASA coeffs 1008 // line is 80 character long for coeffs, so if not, it's temps 1009 if(line.size() != 80) 1010 { 1011 std::vector<std::string> temp_tmp; 1012 SplitString(line," ",temp_tmp,false); 1013 if(temp_tmp.size() == 0) 1014 { 1015 std::cerr << "This line is not understandable by the chemkin parser:\n" << line << std::endl; 1016 antioch_error(); 1017 } 1018 temps.resize(temp_tmp.size(),0.); 1019 for(unsigned int t = 0; t < temp_tmp.size(); t++) 1020 { 1021 temps[t] = std::atof(temp_tmp[t].c_str()); 1022 } 1023 continue; 1024 } 1025 1026 tmp << line.substr(0,18); //name 1027 1028 // this is ChemKin doc, 1st: temps E10.0 1029 tmp << line.substr(45,10); // low 1030 tmp << " " << line.substr(55,10); // high 1031 tmp << " " << line.substr(65,10); // inter 1032 1033 // get rid of last character 1034 for(unsigned int n = 0; n < 3; n++) 1035 { 1036 if(!ascii_getline(_doc,line))// we have a problem 1037 { 1038 std::cerr << "NASA input file error" << std::endl; 1039 antioch_error(); 1040 } 1041 //2nd: coeffs E15.8 1042 tmp << line.substr(0,15) << " " 1043 << line.substr(15,15) << " " 1044 << line.substr(30,15) << " " 1045 << line.substr(45,15) << " " 1046 << line.substr(60,15) << " "; 1047 } 1048 1049 coeffs.clear(); 1050 tmp >> name // Species Name 1051 >> temps[0] 1052 >> temps[2] 1053 >> temps[1]; 1054 coeffs.resize(14,0); 1055 for(unsigned int i = 7; i < 21; i++) 1056 { 1057 NumericType a; 1058 tmp >> a; 1059 coeffs[i%14] = a; 1060 } 1061 1062 // If we are still good, we have a valid set of thermodynamic 1063 // data for this species. Otherwise, we read past end-of-file 1064 // in the section above 1065 if (_doc.good()) 1066 { 1067 // Check if this is a species we want. 1068 if( chem_mixture.species_name_map().find(name) != 1069 chem_mixture.species_name_map().end() ) 1070 { 1071 thermo.add_curve_fit(name, coeffs, temps); 1072 } 1073 } 1074 } // end while 1075 } 1076 1077 } // end namespace Antioch 1078 1079 // Instantiate 1080 ANTIOCH_NUMERIC_TYPE_CLASS_INSTANTIATE(Antioch::ChemKinParser); 1081