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/xml_parser.h" 29 30 // Antioch 31 #include "antioch/chemical_mixture.h" 32 #include "antioch/antioch_numeric_type_instantiate_macro.h" 33 #include "antioch/xml_parser_instantiation_macro.h" 34 #include "antioch/nasa_mixture.h" 35 #include "antioch/cea_curve_fit.h" 36 #include "antioch/nasa7_curve_fit.h" 37 38 //XML 39 #include "antioch/tinyxml2_imp.h" 40 41 // C++ 42 #include <sstream> 43 #include <limits> 44 45 namespace Antioch 46 { 47 template <typename NumericType> XMLParser(const std::string & filename,bool verbose)48 XMLParser<NumericType>::XMLParser(const std::string &filename, bool verbose): 49 ParserBase<NumericType>("XML",filename,verbose), 50 _species_block(NULL), 51 _reaction_block(NULL), 52 _reaction(NULL), 53 _rate_constant(NULL), 54 _Troe(NULL) 55 { 56 _doc = new tinyxml2::XMLDocument; 57 if(_doc->LoadFile(filename.c_str())) 58 { 59 std::cerr << "ERROR: unable to load xml file " << filename << std::endl; 60 std::cerr << "Error of tinyxml2 library:\n" 61 << "\tID = " << _doc->ErrorID() << "\n" 62 << "\tError String1 = " << _doc->GetErrorStr1() << "\n" 63 << "\tError String2 = " << _doc->GetErrorStr2() << std::endl; 64 antioch_error(); 65 } 66 67 if(this->verbose())std::cout << "Having opened file " << filename << std::endl; 68 69 70 // XML block/section names 71 _map[ParsingKey::PHASE_BLOCK] = "phase"; 72 _map[ParsingKey::SPECIES_SET] = "speciesArray"; 73 _map[ParsingKey::SPECIES_DATA] = "speciesData"; 74 _map[ParsingKey::SPECIES] = "species"; 75 _map[ParsingKey::THERMO] = "thermo"; // thermo in ,<speciesData> <species> <thermo> <NASA> <floatArray> </floatArray></NASA> </thermo> </species> </speciesData> 76 77 // Thermo parameters 78 _map[ParsingKey::TMIN] = "Tmin"; 79 _map[ParsingKey::TMAX] = "Tmax"; 80 _map[ParsingKey::NASADATA] = "floatArray"; 81 _map[ParsingKey::NASA7] = "NASA"; 82 _map[ParsingKey::NASA9] = "NASA9"; 83 84 // Kinetics parameters 85 _map[ParsingKey::REACTION_DATA] = "reactionData"; 86 _map[ParsingKey::REACTION] = "reaction"; 87 _map[ParsingKey::REVERSIBLE] = "reversible"; 88 _map[ParsingKey::ID] = "id"; 89 _map[ParsingKey::EQUATION] = "equation"; 90 _map[ParsingKey::CHEMICAL_PROCESS] = "type"; 91 _map[ParsingKey::KINETICS_MODEL] = "rateCoeff"; 92 _map[ParsingKey::REACTANTS] = "reactants"; 93 _map[ParsingKey::PRODUCTS] = "products"; 94 _map[ParsingKey::FORWARD_ORDER] = "ford"; 95 _map[ParsingKey::BACKWARD_ORDER] = "rord"; 96 _map[ParsingKey::PREEXP] = "A"; 97 _map[ParsingKey::POWER] = "b"; 98 _map[ParsingKey::ACTIVATION_ENERGY] = "E"; 99 _map[ParsingKey::BERTHELOT_COEFFICIENT] = "D"; 100 _map[ParsingKey::TREF] = "Tref"; 101 _map[ParsingKey::HV_LAMBDA] = "lambda"; 102 _map[ParsingKey::HV_CROSS_SECTION] = "cross_section"; 103 _map[ParsingKey::UNIT] = "units"; 104 _map[ParsingKey::EFFICIENCY] = "efficiencies"; 105 _map[ParsingKey::FALLOFF_LOW] = "k0"; 106 _map[ParsingKey::FALLOFF_LOW_NAME] = "name"; 107 _map[ParsingKey::TROE_FALLOFF] = "Troe"; 108 _map[ParsingKey::TROE_F_ALPHA] = "alpha"; 109 _map[ParsingKey::TROE_F_TS] = "T1"; 110 _map[ParsingKey::TROE_F_TSS] = "T2"; 111 _map[ParsingKey::TROE_F_TSSS] = "T3"; 112 113 // typically Cantera files list 114 // pre-exponential parameters in (m3/kmol)^(m-1)/s 115 // activation energy in cal/mol, but we want it in K. 116 // power parameter without unit 117 // if falloff, we need to know who's k0 and kinfty 118 // if photochemistry, we have a cross-section on a lambda grid 119 // cross-section typically in cm2/nm (cross-section on a resolution bin, 120 // if bin unit not given, it is lambda unit (supposed to anyway), and a warning message) 121 // lambda typically in nm, sometimes in ang, default considered here is nm 122 // you can also have cm-1, conversion is done with 123 // formulae nm = cm-1 * / * adapted factor 124 _default_unit[ParsingKey::PREEXP] = "m3/kmol"; 125 _default_unit[ParsingKey::POWER] = ""; 126 _default_unit[ParsingKey::ACTIVATION_ENERGY] = "cal/mol"; 127 _default_unit[ParsingKey::BERTHELOT_COEFFICIENT] = "K-1"; 128 _default_unit[ParsingKey::TREF] = "K"; 129 _default_unit[ParsingKey::HV_LAMBDA] = "nm"; 130 _default_unit[ParsingKey::HV_CROSS_SECTION] = "cm2/nm"; 131 _default_unit[ParsingKey::EFFICIENCY] = ""; 132 _default_unit[ParsingKey::TROE_F_ALPHA] = ""; 133 _default_unit[ParsingKey::TROE_F_TS] = "K"; 134 _default_unit[ParsingKey::TROE_F_TSS] = "K"; 135 _default_unit[ParsingKey::TROE_F_TSSS] = "K"; 136 137 //gri30 138 _gri_map[GRI30Comp::FALLOFF] = "falloff"; 139 _gri_map[GRI30Comp::FALLOFF_TYPE] = "type"; 140 _gri_map[GRI30Comp::TROE] = "Troe"; 141 142 this->initialize(); 143 } 144 145 template <typename NumericType> ~XMLParser()146 XMLParser<NumericType>::~XMLParser() 147 { 148 delete _doc; 149 } 150 151 template <typename NumericType> change_file(const std::string & filename)152 void XMLParser<NumericType>::change_file(const std::string & filename) 153 { 154 ParserBase<NumericType>::_file = filename; 155 _species_block = NULL; 156 _reaction_block = NULL; 157 _reaction = NULL; 158 _rate_constant = NULL; 159 _Troe = NULL; 160 161 delete _doc; 162 _doc = new tinyxml2::XMLDocument; 163 if(_doc->LoadFile(filename.c_str())) 164 { 165 std::cerr << "ERROR: unable to load xml file " << filename << std::endl; 166 std::cerr << "Error of tinyxml2 library:\n" 167 << "\tID = " << _doc->ErrorID() << "\n" 168 << "\tError String1 = " << _doc->GetErrorStr1() << "\n" 169 << "\tError String2 = " << _doc->GetErrorStr2() << std::endl; 170 antioch_error(); 171 } 172 173 if(this->verbose())std::cout << "Having opened file " << filename << std::endl; 174 175 this->initialize(); 176 } 177 178 template <typename NumericType> initialize()179 bool XMLParser<NumericType>::initialize() 180 { 181 //we start here 182 _reaction_block = _doc->FirstChildElement("ctml"); 183 if (!_reaction_block) 184 { 185 std::cerr << "ERROR: no <ctml> tag found in input file" 186 << std::endl; 187 antioch_error(); 188 } 189 190 _species_block = _reaction_block->FirstChildElement(_map.at(ParsingKey::PHASE_BLOCK).c_str()); 191 if(_species_block) 192 _species_block = _species_block->FirstChildElement(_map.at(ParsingKey::SPECIES_SET).c_str()); 193 194 _thermo_block = _reaction_block->FirstChildElement(_map.at(ParsingKey::SPECIES_DATA).c_str()); 195 196 _reaction_block = _reaction_block->FirstChildElement(_map.at(ParsingKey::REACTION_DATA).c_str()); 197 198 _reaction = NULL; 199 _rate_constant = NULL; 200 201 return _reaction_block; 202 } 203 204 205 206 template <typename NumericType> species_list()207 const std::vector<std::string> XMLParser<NumericType>::species_list() 208 { 209 if(!_species_block) 210 antioch_error_msg("ERROR: Could not find "+_map.at(ParsingKey::SPECIES_SET)+" section in input file!"); 211 212 std::vector<std::string> molecules; 213 214 split_string(std::string(_species_block->GetText())," ",molecules); 215 remove_newline_from_strings(molecules); 216 217 return molecules; 218 } 219 220 template <typename NumericType> reaction()221 bool XMLParser<NumericType>::reaction() 222 { 223 antioch_assert(_reaction_block); 224 _reaction = (_reaction)? 225 _reaction->NextSiblingElement(_map.at(ParsingKey::REACTION).c_str()): 226 _reaction_block->FirstChildElement(_map.at(ParsingKey::REACTION).c_str()); 227 228 _rate_constant = NULL; 229 _Troe = NULL; 230 231 return _reaction; 232 } 233 234 template <typename NumericType> rate_constant(const std::string & kinetics_model)235 bool XMLParser<NumericType>::rate_constant(const std::string & kinetics_model) 236 { 237 // if in a reaction 238 if(_reaction) 239 { 240 // not the first one 241 if(_rate_constant) 242 { 243 _rate_constant = _rate_constant->NextSiblingElement(kinetics_model.c_str()); 244 }else 245 { 246 // first one, we need to set _rate_constant and _Troe, because they contain environments 247 // we suppose that there is a rateCoeff environement 248 // _rate_constant => <rateCoeff> <kin model> </kin model> </rateCoeff> 249 // _Troe => <rateCoeff> <Troe> </Troe> </rateCoeff> || <rateCoeff><falloff type="Troe"></falloff></rateCoef> 250 antioch_assert(_reaction->FirstChildElement(_map.at(ParsingKey::KINETICS_MODEL).c_str())); 251 _rate_constant = _reaction->FirstChildElement(_map.at(ParsingKey::KINETICS_MODEL).c_str())->FirstChildElement(kinetics_model.c_str()); 252 _Troe = _reaction->FirstChildElement(_map.at(ParsingKey::KINETICS_MODEL).c_str())->FirstChildElement(_map.at(ParsingKey::TROE_FALLOFF).c_str()); 253 if(_Troe == NULL) 254 { 255 _Troe = _reaction->FirstChildElement(_map.at(ParsingKey::KINETICS_MODEL).c_str())->FirstChildElement(_gri_map.at(GRI30Comp::FALLOFF).c_str()); 256 if(_Troe) 257 { 258 if(std::string(_Troe->Attribute(_gri_map.at(GRI30Comp::FALLOFF_TYPE).c_str())).compare(_gri_map.at(GRI30Comp::TROE)) != 0) 259 { 260 _Troe = NULL; 261 } 262 } 263 } 264 } 265 }else 266 { 267 _rate_constant = NULL; 268 } 269 270 return _rate_constant; 271 } 272 273 template <typename NumericType> Troe()274 bool XMLParser<NumericType>::Troe() const 275 { 276 return _Troe; 277 } 278 279 template <typename NumericType> reaction_id()280 const std::string XMLParser<NumericType>::reaction_id() const 281 { 282 std::stringstream id; 283 id << _reaction->Attribute(_map.at(ParsingKey::ID).c_str()); 284 return id.str(); 285 } 286 287 template <typename NumericType> reaction_equation()288 const std::string XMLParser<NumericType>::reaction_equation() const 289 { 290 return _reaction->FirstChildElement(_map.at(ParsingKey::EQUATION).c_str())->GetText(); 291 } 292 293 template <typename NumericType> reaction_chemical_process()294 const std::string XMLParser<NumericType>::reaction_chemical_process() const 295 { 296 const char * chem_proc = _reaction->Attribute(_map.at(ParsingKey::CHEMICAL_PROCESS).c_str()); 297 298 if(chem_proc) 299 { 300 // we're GRI 301 if(std::string(chem_proc).compare(_gri_map.at(GRI30Comp::FALLOFF)) == 0) 302 { 303 // find the falloff block 304 // are we threebody? 305 // if equation contains (+M) 306 const std::string & eq = this->reaction_equation(); 307 // <reaction><rateCoeff><falloff type=''/></rateCoeff></reaction> 308 std::string Lind = "LindemannFalloff"; 309 std::string Troe = "TroeFalloff"; 310 311 // Travis CI gives us regex errors so we'll fall back on this: 312 const char * searchfor = "(+M)"; 313 for (std::size_t pos = eq.find('('); pos < eq.size() && *searchfor; ++pos) 314 { 315 if (eq[pos] == *searchfor) 316 { 317 ++searchfor; 318 continue; 319 } 320 if(eq[pos] ==' ') 321 continue; 322 // If not '(+M)' parenthesis, search for the next parenthesis and check 323 searchfor = "(+M)"; 324 pos = eq.find('(',pos+1); 325 // Adjust position before iteration 326 --pos; 327 } 328 329 // If we found everything then *searchfor is NULL 330 if (!*searchfor) 331 { 332 Lind += "ThreeBody"; 333 Troe += "ThreeBody"; 334 } 335 tinyxml2::XMLElement * fall = _reaction->FirstChildElement(_map.at(ParsingKey::KINETICS_MODEL).c_str())->FirstChildElement(_gri_map.at(GRI30Comp::FALLOFF).c_str()); 336 const char * falloff = fall->Attribute(_map.at(ParsingKey::CHEMICAL_PROCESS).c_str()); 337 const char * cp = (Lind.find(falloff) != std::string::npos) ? Lind.c_str() : 338 (Troe.find(falloff) != std::string::npos) ? Troe.c_str() : chem_proc; 339 340 return (cp)? 341 std::string(cp) : 342 std::string(); 343 } 344 345 } 346 347 return (chem_proc)? 348 std::string(chem_proc) : 349 std::string(); 350 } 351 352 template <typename NumericType> reaction_reversible()353 bool XMLParser<NumericType>::reaction_reversible() const 354 { 355 return (_reaction->Attribute(_map.at(ParsingKey::REVERSIBLE).c_str()))? 356 (std::string(_reaction->Attribute(_map.at(ParsingKey::REVERSIBLE).c_str())) == std::string("no"))?false:true //explicit 357 : 358 true; //default 359 } 360 361 template <typename NumericType> reaction_kinetics_model(const std::vector<std::string> & kinetics_models)362 const std::string XMLParser<NumericType>::reaction_kinetics_model(const std::vector<std::string> &kinetics_models) const 363 { 364 unsigned int imod(0); 365 tinyxml2::XMLElement * rate_constant = _reaction->FirstChildElement(_map.at(ParsingKey::KINETICS_MODEL).c_str())->FirstChildElement(kinetics_models[imod].c_str()); 366 while(!rate_constant) 367 { 368 if(imod == kinetics_models.size() - 1) 369 { 370 std::cerr << "Could not find a suitable kinetics model.\n" 371 << "Implemented kinetics models are:\n"; 372 373 for(unsigned int m = 0; m < kinetics_models.size(); m++) 374 std::cerr << " " << kinetics_models[m] << "\n"; 375 376 std::cerr << "See Antioch documentation for more details." 377 << std::endl; 378 antioch_not_implemented(); 379 } 380 imod++; 381 rate_constant = _reaction->FirstChildElement(_map.at(ParsingKey::KINETICS_MODEL).c_str())->FirstChildElement(kinetics_models[imod].c_str()); 382 } 383 384 return kinetics_models[imod]; 385 } 386 387 template <typename NumericType> reactants_pairs(std::vector<std::pair<std::string,int>> & reactants_pair)388 bool XMLParser<NumericType>::reactants_pairs(std::vector<std::pair<std::string,int> > & reactants_pair) const 389 { 390 tinyxml2::XMLElement* reactants = _reaction->FirstChildElement(_map.at(ParsingKey::REACTANTS).c_str()); 391 return this->molecules_pairs(reactants, reactants_pair); 392 } 393 394 template <typename NumericType> products_pairs(std::vector<std::pair<std::string,int>> & products_pair)395 bool XMLParser<NumericType>::products_pairs(std::vector<std::pair<std::string,int> > & products_pair) const 396 { 397 tinyxml2::XMLElement* products = _reaction->FirstChildElement(_map.at(ParsingKey::PRODUCTS).c_str()); 398 return this->molecules_pairs(products, products_pair); 399 } 400 401 template <typename NumericType> reactants_orders()402 const std::map<std::string,NumericType> XMLParser<NumericType>::reactants_orders() const 403 { 404 tinyxml2::XMLElement* orders = _reaction->FirstChildElement(_map.at(ParsingKey::FORWARD_ORDER).c_str()); 405 std::map<std::string,NumericType> map; 406 if(orders){ 407 std::vector<std::pair<std::string,NumericType> > pairs; 408 if(this->molecules_pairs(orders,pairs)) 409 { 410 for(unsigned int s = 0; s < pairs.size(); s++) 411 { 412 map.insert(pairs[s]); 413 } 414 } 415 } 416 417 return map; 418 } 419 420 template <typename NumericType> products_orders()421 const std::map<std::string,NumericType> XMLParser<NumericType>::products_orders() const 422 { 423 tinyxml2::XMLElement* orders = _reaction->FirstChildElement(_map.at(ParsingKey::BACKWARD_ORDER).c_str()); 424 std::map<std::string,NumericType> map; 425 if(orders){ 426 std::vector<std::pair<std::string,NumericType> > pairs; 427 if(this->molecules_pairs(orders,pairs)) 428 { 429 for(unsigned int s = 0; s < pairs.size(); s++) 430 { 431 map.insert(pairs[s]); 432 } 433 } 434 } 435 return map; 436 } 437 438 439 template <typename NumericType> 440 template <typename PairedType> molecules_pairs(tinyxml2::XMLElement * molecules,std::vector<std::pair<std::string,PairedType>> & molecules_pairs)441 bool XMLParser<NumericType>::molecules_pairs(tinyxml2::XMLElement * molecules, std::vector<std::pair<std::string,PairedType> > & molecules_pairs) const 442 { 443 bool out(true); 444 if(molecules) 445 { 446 447 std::vector<std::string> mol_pairs; 448 449 // Split the reactant string on whitespace. If no entries were found, 450 // there is no whitespace - and assume then only one reactant is listed. 451 split_string(std::string(molecules->GetText()), " ", mol_pairs); 452 453 for( unsigned int p=0; p < mol_pairs.size(); p++ ) 454 { 455 std::pair<std::string,PairedType> pair( split_string_on_colon<PairedType>(mol_pairs[p]) ); 456 457 molecules_pairs.push_back(pair); 458 } 459 }else 460 { 461 out = false; 462 } 463 464 return out; 465 } 466 467 template <typename NumericType> is_k0(unsigned int nrc,const std::string & kin_model)468 bool XMLParser<NumericType>::is_k0(unsigned int nrc, const std::string & kin_model) const 469 { 470 bool k0(false); 471 if(_rate_constant->Attribute(_map.at(ParsingKey::FALLOFF_LOW_NAME).c_str())) 472 { 473 if(std::string(_rate_constant->Attribute(_map.at(ParsingKey::FALLOFF_LOW_NAME).c_str())) == _map.at(ParsingKey::FALLOFF_LOW)) 474 { 475 k0 = true; 476 // now verifying the second one 477 if(nrc == 0) // first reaction rate block 478 { 479 antioch_assert(_rate_constant->NextSiblingElement(kin_model.c_str())); 480 if(_rate_constant->NextSiblingElement(kin_model.c_str())->Attribute(_map.at(ParsingKey::FALLOFF_LOW_NAME).c_str())) // HAHA 481 { 482 std::string error = "I can understand the need to put attributes everywhere, really, but in this case, I'm "; 483 error += "afraid that it's not a good idea to have two \'name\' attributes: only the low pressure limit should have it."; 484 antioch_parsing_error(error); 485 } 486 } 487 }else 488 { 489 std::string error = "The keyword associated with the \'name\' attribute within the description of a falloff should be, and only be, "; 490 error += "\'k0\' to specify the low pressure limit. It seems that the one you provided, \'"; 491 error += std::string(_rate_constant->Attribute(_map.at(ParsingKey::FALLOFF_LOW_NAME).c_str())); 492 error += "\' is not this one. Please correct it at reaction"; 493 error += this->reaction_id(); 494 error += ": "; 495 error += this->reaction_equation(); 496 error += "."; 497 antioch_parsing_error(error); 498 } 499 }else if(nrc == 0) // if we're indeed at the first reading 500 { 501 antioch_assert(_rate_constant->NextSiblingElement(kin_model.c_str())); 502 if(!_rate_constant->NextSiblingElement(kin_model.c_str())->Attribute(_map.at(ParsingKey::FALLOFF_LOW_NAME).c_str())) // and the next doesn't have a name 503 { 504 k0 = true; 505 }else 506 { 507 if(std::string(_rate_constant->NextSiblingElement(kin_model.c_str())->Attribute(_map.at(ParsingKey::FALLOFF_LOW_NAME).c_str())) == _map.at(ParsingKey::FALLOFF_LOW))k0 = false; 508 } 509 } 510 return k0; 511 } 512 513 template <typename NumericType> where_is_k0(const std::string & kin_model)514 unsigned int XMLParser<NumericType>::where_is_k0(const std::string & kin_model) const 515 { 516 antioch_assert(!_rate_constant); //should be done exterior to rate constant block 517 antioch_assert(_reaction); //should be done interior to reaction block 518 519 tinyxml2::XMLElement * rate_constant = _reaction->FirstChildElement(_map.at(ParsingKey::KINETICS_MODEL).c_str()); 520 antioch_assert(rate_constant); 521 rate_constant = rate_constant->FirstChildElement(kin_model.c_str()); 522 unsigned int k0(0); 523 if(rate_constant->NextSiblingElement()->Attribute(_map.at(ParsingKey::FALLOFF_LOW_NAME).c_str())) 524 { 525 if(std::string(rate_constant->NextSiblingElement()->Attribute(_map.at(ParsingKey::FALLOFF_LOW_NAME).c_str())) == _map.at(ParsingKey::FALLOFF_LOW))k0=1; 526 } 527 528 return k0; 529 } 530 531 template <typename NumericType> get_parameter(const tinyxml2::XMLElement * ptr,const std::string & par,NumericType & par_value,std::string & par_unit)532 bool XMLParser<NumericType>::get_parameter(const tinyxml2::XMLElement * ptr, const std::string & par, NumericType & par_value, std::string & par_unit) const 533 { 534 antioch_assert(ptr); 535 536 bool out(false); 537 par_unit.clear(); 538 if(ptr->FirstChildElement(par.c_str())) 539 { 540 par_value = std::atof(ptr->FirstChildElement(par.c_str())->GetText()); 541 if(ptr->FirstChildElement(par.c_str())->Attribute(_map.at(ParsingKey::UNIT).c_str())) 542 par_unit = ptr->FirstChildElement(par.c_str())->Attribute(_map.at(ParsingKey::UNIT).c_str()); 543 out = true; 544 } 545 546 return out; 547 } 548 549 template <typename NumericType> get_parameter(const tinyxml2::XMLElement * ptr,const std::string & par,std::vector<NumericType> & par_values,std::string & par_unit)550 bool XMLParser<NumericType>::get_parameter(const tinyxml2::XMLElement * ptr, const std::string & par, std::vector<NumericType> & par_values, std::string & par_unit) const 551 { 552 antioch_assert(ptr); 553 554 bool out(false); 555 par_unit.clear(); 556 if(ptr->FirstChildElement(par.c_str())) 557 { 558 std::vector<std::string> values; 559 split_string(ptr->FirstChildElement(par.c_str())->GetText()," ",values); 560 561 par_values.resize(values.size()); 562 for(unsigned int i = 0; i < values.size(); i++) 563 par_values[i] = string_to_T<NumericType>(values[i].c_str()); 564 565 if(ptr->FirstChildElement(par.c_str())->Attribute(_map.at(ParsingKey::UNIT).c_str())) 566 par_unit = ptr->FirstChildElement(par.c_str())->Attribute(_map.at(ParsingKey::UNIT).c_str()); 567 568 out = true; 569 } 570 571 return out; 572 } 573 574 template <typename NumericType> find_element_with_attribute(const tinyxml2::XMLElement * element,const std::string & elem_name,const std::string & attribute,const std::string & attr_value)575 tinyxml2::XMLElement* XMLParser<NumericType>::find_element_with_attribute( const tinyxml2::XMLElement * element, 576 const std::string& elem_name, 577 const std::string& attribute, 578 const std::string& attr_value ) const 579 { 580 antioch_assert(element); 581 582 const tinyxml2::XMLElement * elem_with_attr = NULL; 583 584 if( !element->Attribute(attribute.c_str()) ) 585 antioch_error_msg("ERROR: Could not find attribute "+attribute+" for current element!"); 586 587 // First check if the first element has the attribute we're looking for 588 if( std::string( element->Attribute(attribute.c_str()) ) == attr_value ) 589 elem_with_attr = element; 590 591 // Otherwise, look at all the siblings 592 else 593 { 594 elem_with_attr = element->NextSiblingElement(elem_name.c_str()); 595 596 while( elem_with_attr ) 597 { 598 std::string curr_attr; 599 if( elem_with_attr->Attribute(attribute.c_str()) ) 600 curr_attr = std::string(elem_with_attr->Attribute(attribute.c_str())); 601 602 if( curr_attr == attr_value ) 603 break; 604 605 elem_with_attr = elem_with_attr->NextSiblingElement(elem_name.c_str()); 606 } 607 608 // Error out if we couldn't find the attribute with the correct value 609 if( !elem_with_attr ) 610 antioch_error_msg("ERROR: Could not find XMLElement with attribute = "+attribute+" whose value is "+attr_value+"!"); 611 612 } 613 614 return const_cast<tinyxml2::XMLElement*>(elem_with_attr); 615 } 616 617 template <typename NumericType> rate_constant_preexponential_parameter(NumericType & A,std::string & A_unit,std::string & def_unit)618 bool XMLParser<NumericType>::rate_constant_preexponential_parameter(NumericType & A, std::string & A_unit, std::string & def_unit) const 619 { 620 def_unit = _default_unit.at(ParsingKey::PREEXP); 621 return this->get_parameter(_rate_constant,_map.at(ParsingKey::PREEXP).c_str(),A,A_unit); 622 } 623 624 625 template <typename NumericType> rate_constant_power_parameter(NumericType & b,std::string & b_unit,std::string & def_unit)626 bool XMLParser<NumericType>::rate_constant_power_parameter(NumericType & b, std::string & b_unit, std::string & def_unit) const 627 { 628 def_unit = _default_unit.at(ParsingKey::POWER); 629 return this->get_parameter(_rate_constant,_map.at(ParsingKey::POWER).c_str(),b,b_unit); 630 } 631 632 template <typename NumericType> rate_constant_activation_energy_parameter(NumericType & Ea,std::string & Ea_unit,std::string & def_unit)633 bool XMLParser<NumericType>::rate_constant_activation_energy_parameter(NumericType & Ea, std::string & Ea_unit, std::string & def_unit) const 634 { 635 def_unit = _default_unit.at(ParsingKey::ACTIVATION_ENERGY); 636 return this->get_parameter(_rate_constant,_map.at(ParsingKey::ACTIVATION_ENERGY).c_str(),Ea,Ea_unit); 637 } 638 639 template <typename NumericType> rate_constant_Berthelot_coefficient_parameter(NumericType & D,std::string & D_unit,std::string & def_unit)640 bool XMLParser<NumericType>::rate_constant_Berthelot_coefficient_parameter(NumericType & D, std::string & D_unit, std::string & def_unit) const 641 { 642 def_unit = _default_unit.at(ParsingKey::BERTHELOT_COEFFICIENT); 643 return this->get_parameter(_rate_constant,_map.at(ParsingKey::BERTHELOT_COEFFICIENT).c_str(),D,D_unit); 644 } 645 646 template <typename NumericType> rate_constant_lambda_parameter(std::vector<NumericType> & lambda,std::string & lambda_unit,std::string & def_unit)647 bool XMLParser<NumericType>::rate_constant_lambda_parameter(std::vector<NumericType> & lambda, std::string & lambda_unit, std::string & def_unit) const 648 { 649 def_unit = _default_unit.at(ParsingKey::HV_LAMBDA); 650 return this->get_parameter(_rate_constant,_map.at(ParsingKey::HV_LAMBDA).c_str(),lambda,lambda_unit); 651 } 652 653 template <typename NumericType> rate_constant_cross_section_parameter(std::vector<NumericType> & sigma,std::string & sigma_unit,std::string & def_unit)654 bool XMLParser<NumericType>::rate_constant_cross_section_parameter(std::vector<NumericType> & sigma, std::string & sigma_unit, std::string & def_unit) const 655 { 656 def_unit = _default_unit.at(ParsingKey::HV_CROSS_SECTION); 657 return this->get_parameter(_rate_constant,_map.at(ParsingKey::HV_CROSS_SECTION).c_str(),sigma,sigma_unit); 658 } 659 660 template <typename NumericType> rate_constant_Tref_parameter(NumericType & Tref,std::string & Tref_unit,std::string & def_unit)661 bool XMLParser<NumericType>::rate_constant_Tref_parameter(NumericType & Tref, std::string & Tref_unit, std::string & def_unit) const 662 { 663 def_unit = _default_unit.at(ParsingKey::TREF); 664 return this->get_parameter(_rate_constant,_map.at(ParsingKey::TREF).c_str(),Tref,Tref_unit); 665 } 666 667 template <typename NumericType> verify_Kooij_in_place_of_Arrhenius()668 bool XMLParser<NumericType>::verify_Kooij_in_place_of_Arrhenius() const 669 { 670 bool out(false); 671 tinyxml2::XMLElement * rate_constant = _reaction->FirstChildElement(_map.at(ParsingKey::KINETICS_MODEL).c_str()); 672 antioch_assert(rate_constant->FirstChildElement("Arrhenius")); 673 rate_constant = rate_constant->FirstChildElement("Arrhenius"); 674 const char * chem_proc = _reaction->Attribute(_map.at(ParsingKey::CHEMICAL_PROCESS).c_str()); 675 bool i_am_a_falloff = (chem_proc && std::string(chem_proc).compare(_gri_map.at(GRI30Comp::FALLOFF)) == 0); 676 if(rate_constant->FirstChildElement(_map.at(ParsingKey::POWER).c_str())) 677 { 678 if(std::atof(rate_constant->FirstChildElement(_map.at(ParsingKey::POWER).c_str())->GetText()) != 0. || //not a very good test 679 i_am_a_falloff) 680 out = true; 681 } 682 683 return out; 684 } 685 686 template <typename NumericType> efficiencies(std::vector<std::pair<std::string,NumericType>> & par_values)687 bool XMLParser<NumericType>::efficiencies(std::vector<std::pair<std::string,NumericType> > & par_values) const 688 { 689 bool out = false; 690 if(_reaction) 691 { 692 tinyxml2::XMLElement * rate_constant = _reaction->FirstChildElement(_map.at(ParsingKey::KINETICS_MODEL).c_str()); 693 if(rate_constant) 694 { 695 if(rate_constant->FirstChildElement(_map.at(ParsingKey::EFFICIENCY).c_str())) 696 { 697 std::vector<std::string> values; 698 std::string value_string = std::string((rate_constant->FirstChildElement(_map.at(ParsingKey::EFFICIENCY).c_str())->GetText())?rate_constant->FirstChildElement(_map.at(ParsingKey::EFFICIENCY).c_str())->GetText():""); 699 700 split_string(value_string, " ", values); 701 702 for(unsigned int i = 0; i < values.size(); i++) 703 par_values.push_back(split_string_on_colon<NumericType>(values[i])); 704 705 out = true; 706 } 707 } 708 } 709 return out; 710 } 711 712 template <typename NumericType> Troe_GRI_parameter(NumericType & parameter,unsigned int index)713 bool XMLParser<NumericType>::Troe_GRI_parameter(NumericType & parameter, unsigned int index) const 714 { 715 // Troe parameters block 716 // [alpha, T***, T*, T**] 717 std::string Troe_block = std::string(_Troe->GetText()); 718 // do we have something? 719 bool gri = Troe_block.size() != 0; 720 if(gri) 721 { 722 std::vector<std::string> values; 723 split_string(Troe_block, " ", values); 724 // T** is optional, we generalize ?? should we ?? 725 if(index < values.size()) 726 { 727 parameter = std::atof(values[index].c_str()); 728 }else{ 729 gri = false; 730 } 731 } 732 return gri; 733 } 734 735 template <typename NumericType> Troe_alpha_parameter(NumericType & alpha,std::string & alpha_unit,std::string & def_unit)736 bool XMLParser<NumericType>::Troe_alpha_parameter(NumericType & alpha, std::string & alpha_unit, std::string & def_unit) const 737 { 738 def_unit = _default_unit.at(ParsingKey::TROE_F_ALPHA); 739 // if defined as Antioch wants it 740 bool antioch = this->get_parameter(_Troe,_map.at(ParsingKey::TROE_F_ALPHA),alpha,alpha_unit); 741 // if not, we try GRI30 742 return antioch ? antioch : this->Troe_GRI_parameter(alpha,0); 743 } 744 745 template <typename NumericType> Troe_T1_parameter(NumericType & T1,std::string & T1_unit,std::string & def_unit)746 bool XMLParser<NumericType>::Troe_T1_parameter(NumericType & T1, std::string & T1_unit, std::string & def_unit) const 747 { 748 def_unit = _default_unit.at(ParsingKey::TROE_F_TS); 749 // if defined as Antioch wants it 750 bool antioch = this->get_parameter(_Troe,_map.at(ParsingKey::TROE_F_TS),T1,T1_unit); 751 // if not, we try GRI30 752 return antioch ? antioch : this->Troe_GRI_parameter(T1,2); 753 } 754 755 template <typename NumericType> Troe_T2_parameter(NumericType & T2,std::string & T2_unit,std::string & def_unit)756 bool XMLParser<NumericType>::Troe_T2_parameter(NumericType & T2, std::string & T2_unit, std::string & def_unit) const 757 { 758 def_unit = _default_unit.at(ParsingKey::TROE_F_TSS); 759 // if defined as Antioch wants it 760 bool antioch = this->get_parameter(_Troe,_map.at(ParsingKey::TROE_F_TSS),T2,T2_unit); 761 // if not, we try GRI30 762 return antioch ? antioch : this->Troe_GRI_parameter(T2,3); 763 } 764 765 template <typename NumericType> Troe_T3_parameter(NumericType & T3,std::string & T3_unit,std::string & def_unit)766 bool XMLParser<NumericType>::Troe_T3_parameter(NumericType & T3, std::string & T3_unit, std::string & def_unit) const 767 { 768 def_unit = _default_unit.at(ParsingKey::TROE_F_TSSS); 769 // if defined as Antioch wants it 770 bool antioch = this->get_parameter(_Troe,_map.at(ParsingKey::TROE_F_TSSS),T3,T3_unit); 771 // if not, we try GRI30 772 return antioch ? antioch : this->Troe_GRI_parameter(T3,1); 773 } 774 775 776 template <typename NumericType> 777 template <typename ThermoType> read_thermodynamic_data_root(ThermoType & thermo)778 void XMLParser<NumericType>::read_thermodynamic_data_root(ThermoType & thermo) 779 { 780 if(!_thermo_block) 781 antioch_error_msg("ERROR: No "+_map.at(ParsingKey::SPECIES_DATA)+" section found! Cannot parse thermo!"); 782 783 const ChemicalMixture<NumericType> & chem_mixture = thermo.chemical_mixture(); 784 const std::vector<ChemicalSpecies<NumericType>*>& chem_species = chem_mixture.chemical_species(); 785 786 // Based on the ThermoType, namely the CurveFit, we deduce what the section name is. 787 std::string nasa_xml_section = this->nasa_xml_section(thermo); 788 789 for(unsigned int s = 0; s < chem_mixture.n_species(); s++) 790 { 791 // Step to first species block 792 tinyxml2::XMLElement * species_block = _thermo_block->FirstChildElement(_map.at(ParsingKey::SPECIES).c_str()); 793 794 if(!species_block) 795 antioch_error_msg("ERROR: No "+_map.at(ParsingKey::SPECIES)+" block found within "+_map.at(ParsingKey::SPECIES_DATA)+" section! Cannot parse thermo!"); 796 797 798 const std::string& name = chem_species[s]->species(); 799 800 tinyxml2::XMLElement * spec = NULL; 801 spec = this->find_element_with_attribute( species_block, 802 _map.at(ParsingKey::SPECIES), 803 "name", 804 name ); 805 806 if(!spec) 807 antioch_error_msg("ERROR: Species "+name+" has not been found in the "+_map.at(ParsingKey::SPECIES_DATA)+" section! Cannot parse thermo!"); 808 809 else 810 { 811 spec = spec->FirstChildElement(_map.at(ParsingKey::THERMO).c_str()); 812 813 if(!spec) 814 antioch_error_msg("ERROR: No "+_map.at(ParsingKey::THERMO)+" block found for species "+name+"! Cannot parse thermo!"); 815 816 // containers for parsing thermo data 817 std::vector<NumericType> temps; 818 std::vector<NumericType> values; 819 tinyxml2::XMLElement * coeffs; 820 std::vector<std::string> coeffs_str; 821 822 // looping for each of the temperature intervals for this species 823 tinyxml2::XMLElement * nasa = spec->FirstChildElement(nasa_xml_section.c_str()); 824 if(!nasa) 825 antioch_error_msg("ERROR: Could not find "+nasa_xml_section+" thermo section!"); 826 827 while(nasa) 828 { 829 if( !(nasa->Attribute(_map.at(ParsingKey::TMIN).c_str())) ) 830 antioch_error_msg("ERROR: Could not find "+_map.at(ParsingKey::TMIN)+" attribute for species "+name+"!"); 831 832 if( !(nasa->Attribute(_map.at(ParsingKey::TMAX).c_str())) ) 833 antioch_error_msg("ERROR: Could not find "+_map.at(ParsingKey::TMAX)+" attribute for species "+name+"!"); 834 // By convention, we put the first TMIN in, and then only the TMAX thereafter 835 // We have a consistency check below to make sure the TMIN's in the input are consistent 836 if( temps.empty() ) 837 temps.push_back(string_to_T<NumericType>(nasa->Attribute(_map.at(ParsingKey::TMIN).c_str()))); 838 839 // temperatures, only Tmax as Tmin is suppose to be last Tmax 840 temps.push_back(string_to_T<NumericType>(nasa->Attribute(_map.at(ParsingKey::TMAX).c_str()))); 841 842 // now coeffs 843 if( !(nasa->FirstChildElement(_map.at(ParsingKey::NASADATA).c_str())) ) 844 antioch_error_msg("ERROR: Could not find "+_map.at(ParsingKey::NASADATA)+" data for species "+name+"!"); 845 846 coeffs = nasa->FirstChildElement(_map.at(ParsingKey::NASADATA).c_str()); 847 split_string(std::string(coeffs->GetText())," ",coeffs_str); 848 remove_newline_from_strings(coeffs_str); 849 850 for(unsigned int d = 0; d < coeffs_str.size(); d++) 851 values.push_back(string_to_T<NumericType>(coeffs_str[d])); 852 853 // If we have more than one interval, make sure the temperature intervals match up 854 if( temps.size() > 1 ) 855 { 856 NumericType prev_Tmax = *(temps.end()-2); 857 NumericType Tmin = string_to_T<NumericType>(nasa->Attribute(_map.at(ParsingKey::TMIN).c_str())); 858 NumericType diff = (Tmin - prev_Tmax)/prev_Tmax; 859 860 const NumericType tol = std::numeric_limits<NumericType>::epsilon() * 10.; 861 862 if(std::abs(diff) > tol) 863 antioch_error_msg("ERROR: Tmax/Tmin mismatch for species "+name+"!"); 864 } 865 866 // This is meant to store only data one interval at a time, so we must clear at each iteration 867 coeffs_str.clear(); 868 869 // Move onto next interval of data 870 nasa = nasa->NextSiblingElement(nasa_xml_section.c_str()); 871 872 } // end while loop 873 874 thermo.add_curve_fit(name, values, temps); 875 } 876 877 } // end species loop 878 } 879 880 // Instantiate 881 ANTIOCH_NUMERIC_TYPE_CLASS_INSTANTIATE(XMLParser); 882 ANTIOCH_XML_PARSER_INSTANTIATE(); 883 884 } // end namespace Antioch 885