src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "chemkinReader.H"
  27 #include <fstream>
  28 #include "atomicWeights.H"
  29 #include "IrreversibleReaction.H"
  30 #include "ReversibleReaction.H"
  31 #include "NonEquilibriumReversibleReaction.H"
  32 #include "ArrheniusReactionRate.H"
  33 #include "thirdBodyArrheniusReactionRate.H"
  34 #include "FallOffReactionRate.H"
  35 #include "ChemicallyActivatedReactionRate.H"
  36 #include "LindemannFallOffFunction.H"
  37 #include "TroeFallOffFunction.H"
  38 #include "SRIFallOffFunction.H"
  39 #include "LandauTellerReactionRate.H"
  40 #include "JanevReactionRate.H"
  41 #include "powerSeriesReactionRate.H"
  42 #include "addToRunTimeSelectionTable.H"
  43
  44 /* * * * * * * * * * * * * * * * * Static data * * * * * * * * * * * * * * * */
  45
  46 namespace Foam
  47 {
  48     addChemistryReaderType(chemkinReader, gasThermoPhysics);
  49 };
  50
  51 /* * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * */
  52
  53 const char* Foam::chemkinReader::reactionTypeNames[4] =
  54 {
  55     "irreversible",
  56     "reversible",
  57     "nonEquilibriumReversible",
  58     "unknownReactionType"
  59 };
  60
  61 const char* Foam::chemkinReader::reactionRateTypeNames[8] =
  62 {
  63     "Arrhenius",
  64     "thirdBodyArrhenius",
  65     "unimolecularFallOff",
  66     "chemicallyActivatedBimolecular",
  67     "LandauTeller",
  68     "Janev",
  69     "powerSeries",
  70     "unknownReactionRateType"
  71 };
  72
  73 const char* Foam::chemkinReader::fallOffFunctionNames[4] =
  74 {
  75     "Lindemann",
  76     "Troe",
  77     "SRI",
  78     "unknownFallOffFunctionType"
  79 };
  80
  81 void Foam::chemkinReader::initReactionKeywordTable()
  82 {
  83     reactionKeywordTable_.insert("M", thirdBodyReactionType);
  84     reactionKeywordTable_.insert("LOW", unimolecularFallOffReactionType);
  85     reactionKeywordTable_.insert("HIGH", chemicallyActivatedBimolecularReactionType);
  86     reactionKeywordTable_.insert("TROE", TroeReactionType);
  87     reactionKeywordTable_.insert("SRI", SRIReactionType);
  88     reactionKeywordTable_.insert("LT", LandauTellerReactionType);
  89     reactionKeywordTable_.insert("RLT", reverseLandauTellerReactionType);
  90     reactionKeywordTable_.insert("JAN", JanevReactionType);
  91     reactionKeywordTable_.insert("FIT1", powerSeriesReactionRateType);
  92     reactionKeywordTable_.insert("HV", radiationActivatedReactionType);
  93     reactionKeywordTable_.insert("TDEP", speciesTempReactionType);
  94     reactionKeywordTable_.insert("EXCI", energyLossReactionType);
  95     reactionKeywordTable_.insert("MOME", plasmaMomentumTransfer);
  96     reactionKeywordTable_.insert("XSMI", collisionCrossSection);
  97     reactionKeywordTable_.insert("REV", nonEquilibriumReversibleReactionType);
  98     reactionKeywordTable_.insert("DUPLICATE", duplicateReactionType);
  99     reactionKeywordTable_.insert("DUP", duplicateReactionType);
 100     reactionKeywordTable_.insert("FORD", speciesOrderForward);
 101     reactionKeywordTable_.insert("RORD", speciesOrderReverse);
 102     reactionKeywordTable_.insert("UNITS", UnitsOfReaction);
 103     reactionKeywordTable_.insert("END", end);
 104 }
 105
 106
 107 Foam::scalar Foam::chemkinReader::molecularWeight
 108 (
 109     const List<specieElement>& specieComposition
 110 ) const
 111 {
 112     scalar molWt = 0.0;
 113
 114     forAll (specieComposition, i)
 115     {
 116         label nAtoms = specieComposition[i].nAtoms;
 117         const word& elementName = specieComposition[i].elementName;
 118
 119         if (isotopeAtomicWts_.found(elementName))
 120         {
 121             molWt += nAtoms*isotopeAtomicWts_[elementName];
 122         }
 123         else if (atomicWeights.found(elementName))
 124         {
 125             molWt += nAtoms*atomicWeights[elementName];
 126         }
 127         else
 128         {
 129             FatalErrorIn("chemkinReader::lex()")
 130                 << "Unknown element " << elementName
 131                 << " on line " << lineNo_-1 << nl
 132                 << "    specieComposition: " << specieComposition
 133                 << exit(FatalError);
 134         }
 135     }
 136
 137     return molWt;
 138 }
 139
 140
 141 void Foam::chemkinReader::checkCoeffs
 142 (
 143     const scalarList& reactionCoeffs,
 144     const char* reactionRateName,
 145     const label nCoeffs
 146 ) const
 147 {
 148     if (reactionCoeffs.size() != nCoeffs)
 149     {
 150         FatalErrorIn("chemkinReader::checkCoeffs")
 151             << "Wrong number of coefficients for the " << reactionRateName
 152             << " rate expression on line "
 153             << lineNo_-1 << ", should be "
 154             << nCoeffs << " but " << reactionCoeffs.size() << " supplied." << nl
 155             << "Coefficients are "
 156             << reactionCoeffs << nl
 157             << exit(FatalError);
 158     }
 159 }
 160
 161 template<class ReactionRateType>
 162 void Foam::chemkinReader::addReactionType
 163 (
 164     const reactionType rType,
 165     DynamicList<gasReaction::specieCoeffs>& lhs,
 166     DynamicList<gasReaction::specieCoeffs>& rhs,
 167     const ReactionRateType& rr
 168 )
 169 {
 170     switch (rType)
 171     {
 172         case irreversible:
 173         {
 174             reactions_.append
 175             (
 176                 new IrreversibleReaction<gasThermoPhysics, ReactionRateType>
 177                 (
 178                     Reaction<gasThermoPhysics>
 179                     (
 180                         speciesTable_,
 181                         lhs.shrink(),
 182                         rhs.shrink(),
 183                         speciesThermo_
 184                     ),
 185                     rr
 186                 )
 187             );
 188         }
 189         break;
 190
 191         case reversible:
 192         {
 193             reactions_.append
 194             (
 195                 new ReversibleReaction<gasThermoPhysics, ReactionRateType>
 196                 (
 197                     Reaction<gasThermoPhysics>
 198                     (
 199                         speciesTable_,
 200                         lhs.shrink(),
 201                         rhs.shrink(),
 202                         speciesThermo_
 203                     ),
 204                     rr
 205                 )
 206             );
 207         }
 208         break;
 209
 210         default:
 211
 212             if (rType < 3)
 213             {
 214                 FatalErrorIn("chemkinReader::addReactionType")
 215                     << "Reaction type " << reactionTypeNames[rType]
 216                     << " on line " << lineNo_-1
 217                     << " not handled by this function"
 218                     << exit(FatalError);
 219             }
 220             else
 221             {
 222                 FatalErrorIn("chemkinReader::addReactionType")
 223                     << "Unknown reaction type " << rType
 224                     << " on line " << lineNo_-1
 225                     << exit(FatalError);
 226             }
 227     }
 228 }
 229
 230 template<template<class, class> class PressureDependencyType>
 231 void Foam::chemkinReader::addPressureDependentReaction
 232 (
 233     const reactionType rType,
 234     const fallOffFunctionType fofType,
 235     DynamicList<gasReaction::specieCoeffs>& lhs,
 236     DynamicList<gasReaction::specieCoeffs>& rhs,
 237     const scalarList& efficiencies,
 238     const scalarList& k0Coeffs,
 239     const scalarList& kInfCoeffs,
 240     const HashTable<scalarList>& reactionCoeffsTable,
 241     const scalar Afactor0,
 242     const scalar AfactorInf,
 243     const scalar RR
 244 )
 245 {
 246     checkCoeffs(k0Coeffs, "k0", 3);
 247     checkCoeffs(kInfCoeffs, "kInf", 3);
 248
 249     switch (fofType)
 250     {
 251         case Lindemann:
 252         {
 253             addReactionType
 254             (
 255                 rType,
 256                 lhs, rhs,
 257                 PressureDependencyType
 258                     <ArrheniusReactionRate, LindemannFallOffFunction>
 259                 (
 260                     ArrheniusReactionRate
 261                     (
 262                         Afactor0*k0Coeffs[0],
 263                         k0Coeffs[1],
 264                         k0Coeffs[2]/RR
 265                     ),
 266                     ArrheniusReactionRate
 267                     (
 268                         AfactorInf*kInfCoeffs[0],
 269                         kInfCoeffs[1],
 270                         kInfCoeffs[2]/RR
 271                     ),
 272                     LindemannFallOffFunction(),
 273                     thirdBodyEfficiencies(speciesTable_, efficiencies)
 274                 )
 275             );
 276         }
 277         break;
 278
 279         case Troe:
 280         {
 281             scalarList TroeCoeffs
 282             (
 283                 reactionCoeffsTable[fallOffFunctionNames[fofType]]
 284             );
 285
 286             if (TroeCoeffs.size() != 4 && TroeCoeffs.size() != 3)
 287             {
 288                 FatalErrorIn("chemkinReader::addPressureDependentReaction")
 289                     << "Wrong number of coefficients for Troe rate expression"
 290                        " on line " << lineNo_-1 << ", should be 3 or 4 but "
 291                     << TroeCoeffs.size() << " supplied." << nl
 292                     << "Coefficients are "
 293                     << TroeCoeffs << nl
 294                     << exit(FatalError);
 295             }
 296
 297             if (TroeCoeffs.size() == 3)
 298             {
 299                 TroeCoeffs.setSize(4);
 300                 TroeCoeffs[3] = GREAT;
 301             }
 302
 303             addReactionType
 304             (
 305                 rType,
 306                 lhs, rhs,
 307                 PressureDependencyType
 308                     <ArrheniusReactionRate, TroeFallOffFunction>
 309                 (
 310                     ArrheniusReactionRate
 311                     (
 312                         Afactor0*k0Coeffs[0],
 313                         k0Coeffs[1],
 314                         k0Coeffs[2]/RR
 315                     ),
 316                     ArrheniusReactionRate
 317                     (
 318                         AfactorInf*kInfCoeffs[0],
 319                         kInfCoeffs[1],
 320                         kInfCoeffs[2]/RR
 321                     ),
 322                     TroeFallOffFunction
 323                     (
 324                         TroeCoeffs[0],
 325                         TroeCoeffs[1],
 326                         TroeCoeffs[2],
 327                         TroeCoeffs[3]
 328                     ),
 329                     thirdBodyEfficiencies(speciesTable_, efficiencies)
 330                 )
 331             );
 332         }
 333         break;
 334
 335         case SRI:
 336         {
 337             scalarList SRICoeffs
 338             (
 339                 reactionCoeffsTable[fallOffFunctionNames[fofType]]
 340             );
 341
 342             if (SRICoeffs.size() != 5 && SRICoeffs.size() != 3)
 343             {
 344                 FatalErrorIn("chemkinReader::addPressureDependentReaction")
 345                     << "Wrong number of coefficients for SRI rate expression"
 346                        " on line " << lineNo_-1 << ", should be 3 or 5 but "
 347                     << SRICoeffs.size() << " supplied." << nl
 348                     << "Coefficients are "
 349                     << SRICoeffs << nl
 350                     << exit(FatalError);
 351             }
 352
 353             if (SRICoeffs.size() == 3)
 354             {
 355                 SRICoeffs.setSize(5);
 356                 SRICoeffs[3] = 1.0;
 357                 SRICoeffs[4] = 0.0;
 358             }
 359
 360             addReactionType
 361             (
 362                 rType,
 363                 lhs, rhs,
 364                 PressureDependencyType
 365                     <ArrheniusReactionRate, SRIFallOffFunction>
 366                 (
 367                     ArrheniusReactionRate
 368                     (
 369                         Afactor0*k0Coeffs[0],
 370                         k0Coeffs[1],
 371                         k0Coeffs[2]/RR
 372                     ),
 373                     ArrheniusReactionRate
 374                     (
 375                         AfactorInf*kInfCoeffs[0],
 376                         kInfCoeffs[1],
 377                         kInfCoeffs[2]/RR
 378                     ),
 379                     SRIFallOffFunction
 380                     (
 381                         SRICoeffs[0],
 382                         SRICoeffs[1],
 383                         SRICoeffs[2],
 384                         SRICoeffs[3],
 385                         SRICoeffs[4]
 386                     ),
 387                     thirdBodyEfficiencies(speciesTable_, efficiencies)
 388                 )
 389             );
 390         }
 391         break;
 392
 393         default:
 394         {
 395             if (fofType < 4)
 396             {
 397                 FatalErrorIn("chemkinReader::addPressureDependentReaction")
 398                     << "Fall-off function type "
 399                     << fallOffFunctionNames[fofType]
 400                     << " on line " << lineNo_-1
 401                     << " not implemented"
 402                     << exit(FatalError);
 403             }
 404             else
 405             {
 406                 FatalErrorIn("chemkinReader::addPressureDependentReaction")
 407                     << "Unknown fall-off function type " << fofType
 408                     << " on line " << lineNo_-1
 409                     << exit(FatalError);
 410             }
 411         }
 412     }
 413 }
 414
 415
 416 void Foam::chemkinReader::addReaction
 417 (
 418     DynamicList<gasReaction::specieCoeffs>& lhs,
 419     DynamicList<gasReaction::specieCoeffs>& rhs,
 420     const scalarList& efficiencies,
 421     const reactionType rType,
 422     const reactionRateType rrType,
 423     const fallOffFunctionType fofType,
 424     const scalarList& ArrheniusCoeffs,
 425     HashTable<scalarList>& reactionCoeffsTable,
 426     const scalar RR
 427 )
 428 {
 429     checkCoeffs(ArrheniusCoeffs, "Arrhenius", 3);
 430
 431     scalarList nAtoms(elementNames_.size(), 0.0);
 432
 433     forAll (lhs, i)
 434     {
 435         const List<specieElement>& specieComposition =
 436             specieComposition_[speciesTable_[lhs[i].index]];
 437
 438         forAll (specieComposition, j)
 439         {
 440             label elementi = elementIndices_[specieComposition[j].elementName];
 441             nAtoms[elementi] += lhs[i].stoichCoeff*specieComposition[j].nAtoms;
 442         }
 443     }
 444
 445     forAll (rhs, i)
 446     {
 447         const List<specieElement>& specieComposition =
 448             specieComposition_[speciesTable_[rhs[i].index]];
 449
 450         forAll (specieComposition, j)
 451         {
 452             label elementi = elementIndices_[specieComposition[j].elementName];
 453             nAtoms[elementi] -= rhs[i].stoichCoeff*specieComposition[j].nAtoms;
 454         }
 455     }
 456
 457
 458     // Calculate the unit conversion factor for the A coefficient
 459     // for the change from mol/cm^3 to kmol/m^3 concentraction units
 460     const scalar concFactor = 0.001;
 461     scalar sumExp = 0.0;
 462     forAll (lhs, i)
 463     {
 464         sumExp += lhs[i].exponent;
 465     }
 466     scalar Afactor = pow(concFactor, sumExp - 1.0);
 467
 468     scalar AfactorRev = Afactor;
 469
 470     if (rType == nonEquilibriumReversible)
 471     {
 472         sumExp = 0.0;
 473         forAll (rhs, i)
 474         {
 475             sumExp += rhs[i].exponent;
 476         }
 477         AfactorRev = pow(concFactor, sumExp - 1.0);
 478     }
 479
 480     switch (rrType)
 481     {
 482         case Arrhenius:
 483         {
 484             if (rType == nonEquilibriumReversible)
 485             {
 486                 const scalarList& reverseArrheniusCoeffs =
 487                     reactionCoeffsTable[reactionTypeNames[rType]];
 488
 489                 checkCoeffs(reverseArrheniusCoeffs, "reverse Arrhenius", 3);
 490
 491                 reactions_.append
 492                 (
 493                     new NonEquilibriumReversibleReaction
 494                         <gasThermoPhysics, ArrheniusReactionRate>
 495                     (
 496                         Reaction<gasThermoPhysics>
 497                         (
 498                             speciesTable_,
 499                             lhs.shrink(),
 500                             rhs.shrink(),
 501                             speciesThermo_
 502                         ),
 503                         ArrheniusReactionRate
 504                         (
 505                             Afactor*ArrheniusCoeffs[0],
 506                             ArrheniusCoeffs[1],
 507                             ArrheniusCoeffs[2]/RR
 508                         ),
 509                         ArrheniusReactionRate
 510                         (
 511                             AfactorRev*reverseArrheniusCoeffs[0],
 512                             reverseArrheniusCoeffs[1],
 513                             reverseArrheniusCoeffs[2]/RR
 514                         )
 515                     )
 516                 );
 517             }
 518             else
 519             {
 520                 addReactionType
 521                 (
 522                     rType,
 523                     lhs, rhs,
 524                     ArrheniusReactionRate
 525                     (
 526                         Afactor*ArrheniusCoeffs[0],
 527                         ArrheniusCoeffs[1],
 528                         ArrheniusCoeffs[2]/RR
 529                     )
 530                 );
 531             }
 532         }
 533         break;
 534
 535         case thirdBodyArrhenius:
 536         {
 537             if (rType == nonEquilibriumReversible)
 538             {
 539                 const scalarList& reverseArrheniusCoeffs =
 540                     reactionCoeffsTable[reactionTypeNames[rType]];
 541
 542                 checkCoeffs(reverseArrheniusCoeffs, "reverse Arrhenius", 3);
 543
 544                 reactions_.append
 545                 (
 546                     new NonEquilibriumReversibleReaction
 547                         <gasThermoPhysics, thirdBodyArrheniusReactionRate>
 548                     (
 549                         Reaction<gasThermoPhysics>
 550                         (
 551                             speciesTable_,
 552                             lhs.shrink(),
 553                             rhs.shrink(),
 554                             speciesThermo_
 555                         ),
 556                         thirdBodyArrheniusReactionRate
 557                         (
 558                             Afactor*concFactor*ArrheniusCoeffs[0],
 559                             ArrheniusCoeffs[1],
 560                             ArrheniusCoeffs[2]/RR,
 561                             thirdBodyEfficiencies(speciesTable_, efficiencies)
 562                         ),
 563                         thirdBodyArrheniusReactionRate
 564                         (
 565                             AfactorRev*concFactor*reverseArrheniusCoeffs[0],
 566                             reverseArrheniusCoeffs[1],
 567                             reverseArrheniusCoeffs[2]/RR,
 568                             thirdBodyEfficiencies(speciesTable_, efficiencies)
 569                         )
 570                     )
 571                 );
 572             }
 573             else
 574             {
 575                 addReactionType
 576                 (
 577                     rType,
 578                     lhs, rhs,
 579                     thirdBodyArrheniusReactionRate
 580                     (
 581                         Afactor*concFactor*ArrheniusCoeffs[0],
 582                         ArrheniusCoeffs[1],
 583                         ArrheniusCoeffs[2]/RR,
 584                         thirdBodyEfficiencies(speciesTable_, efficiencies)
 585                     )
 586                 );
 587             }
 588         }
 589         break;
 590
 591         case unimolecularFallOff:
 592         {
 593             addPressureDependentReaction<FallOffReactionRate>
 594             (
 595                 rType,
 596                 fofType,
 597                 lhs,
 598                 rhs,
 599                 efficiencies,
 600                 reactionCoeffsTable[reactionRateTypeNames[rrType]],
 601                 ArrheniusCoeffs,
 602                 reactionCoeffsTable,
 603                 concFactor*Afactor,
 604                 Afactor,
 605                 RR
 606             );
 607         }
 608         break;
 609
 610         case chemicallyActivatedBimolecular:
 611         {
 612             addPressureDependentReaction<ChemicallyActivatedReactionRate>
 613             (
 614                 rType,
 615                 fofType,
 616                 lhs,
 617                 rhs,
 618                 efficiencies,
 619                 ArrheniusCoeffs,
 620                 reactionCoeffsTable[reactionRateTypeNames[rrType]],
 621                 reactionCoeffsTable,
 622                 Afactor,
 623                 Afactor/concFactor,
 624                 RR
 625             );
 626         }
 627         break;
 628
 629         case LandauTeller:
 630         {
 631             const scalarList& LandauTellerCoeffs =
 632                 reactionCoeffsTable[reactionRateTypeNames[rrType]];
 633             checkCoeffs(LandauTellerCoeffs, "Landau-Teller", 2);
 634
 635             if (rType == nonEquilibriumReversible)
 636             {
 637                 const scalarList& reverseArrheniusCoeffs =
 638                     reactionCoeffsTable[reactionTypeNames[rType]];
 639                 checkCoeffs(reverseArrheniusCoeffs, "reverse Arrhenius", 3);
 640
 641                 const scalarList& reverseLandauTellerCoeffs =
 642                     reactionCoeffsTable
 643                     [
 644                         word(reactionTypeNames[rType])
 645                       + reactionRateTypeNames[rrType]
 646                     ];
 647                 checkCoeffs(LandauTellerCoeffs, "reverse Landau-Teller", 2);
 648
 649                 reactions_.append
 650                 (
 651                     new NonEquilibriumReversibleReaction
 652                         <gasThermoPhysics, LandauTellerReactionRate>
 653                     (
 654                         Reaction<gasThermoPhysics>
 655                         (
 656                             speciesTable_,
 657                             lhs.shrink(),
 658                             rhs.shrink(),
 659                             speciesThermo_
 660                         ),
 661                         LandauTellerReactionRate
 662                         (
 663                             Afactor*ArrheniusCoeffs[0],
 664                             ArrheniusCoeffs[1],
 665                             ArrheniusCoeffs[2]/RR,
 666                             LandauTellerCoeffs[0],
 667                             LandauTellerCoeffs[1]
 668                         ),
 669                         LandauTellerReactionRate
 670                         (
 671                             AfactorRev*reverseArrheniusCoeffs[0],
 672                             reverseArrheniusCoeffs[1],
 673                             reverseArrheniusCoeffs[2]/RR,
 674                             reverseLandauTellerCoeffs[0],
 675                             reverseLandauTellerCoeffs[1]
 676                         )
 677                     )
 678                 );
 679             }
 680             else
 681             {
 682                 addReactionType
 683                 (
 684                     rType,
 685                     lhs, rhs,
 686                     LandauTellerReactionRate
 687                     (
 688                         Afactor*ArrheniusCoeffs[0],
 689                         ArrheniusCoeffs[1],
 690                         ArrheniusCoeffs[2]/RR,
 691                         LandauTellerCoeffs[0],
 692                         LandauTellerCoeffs[1]
 693                     )
 694                 );
 695             }
 696         }
 697         break;
 698
 699         case Janev:
 700         {
 701             const scalarList& JanevCoeffs =
 702                 reactionCoeffsTable[reactionRateTypeNames[rrType]];
 703
 704             checkCoeffs(JanevCoeffs, "Janev", 9);
 705
 706             addReactionType
 707             (
 708                 rType,
 709                 lhs, rhs,
 710                 JanevReactionRate
 711                 (
 712                     Afactor*ArrheniusCoeffs[0],
 713                     ArrheniusCoeffs[1],
 714                     ArrheniusCoeffs[2]/RR,
 715                     JanevCoeffs.begin()
 716                 )
 717             );
 718         }
 719         break;
 720
 721         case powerSeries:
 722         {
 723             const scalarList& powerSeriesCoeffs =
 724                 reactionCoeffsTable[reactionRateTypeNames[rrType]];
 725
 726             checkCoeffs(powerSeriesCoeffs, "power-series", 4);
 727
 728             addReactionType
 729             (
 730                 rType,
 731                 lhs, rhs,
 732                 powerSeriesReactionRate
 733                 (
 734                     Afactor*ArrheniusCoeffs[0],
 735                     ArrheniusCoeffs[1],
 736                     ArrheniusCoeffs[2]/RR,
 737                     powerSeriesCoeffs.begin()
 738                 )
 739             );
 740         }
 741         break;
 742
 743         case unknownReactionRateType:
 744         {
 745             FatalErrorIn("chemkinReader::addReaction")
 746                 << "Internal error on line " << lineNo_-1
 747                 << ": reaction rate type has not been set"
 748                 << exit(FatalError);
 749         }
 750         break;
 751
 752         default:
 753         {
 754             if (rrType < 9)
 755             {
 756                 FatalErrorIn("chemkinReader::addReaction")
 757                     << "Reaction rate type " << reactionRateTypeNames[rrType]
 758                     << " on line " << lineNo_-1
 759                     << " not implemented"
 760                     << exit(FatalError);
 761             }
 762             else
 763             {
 764                 FatalErrorIn("chemkinReader::addReaction")
 765                     << "Unknown reaction rate type " << rrType
 766                     << " on line " << lineNo_-1
 767                     << exit(FatalError);
 768             }
 769         }
 770     }
 771
 772
 773     forAll (nAtoms, i)
 774     {
 775         if (mag(nAtoms[i]) > SMALL)
 776         {
 777             FatalErrorIn("chemkinReader::addReaction")
 778                 << "Elemental imbalance in " << elementNames_[i]
 779                 << " in reaction" << nl
 780                 << reactions_.last() << nl
 781                 << " on line " << lineNo_-1
 782                 << exit(FatalError);
 783         }
 784     }
 785
 786     lhs.clear();
 787     rhs.clear();
 788     reactionCoeffsTable.clear();
 789 }
 790
 791
 792 void Foam::chemkinReader::read
 793 (
 794     const fileName& CHEMKINFileName,
 795     const fileName& thermoFileName
 796 )
 797 {
 798     if (thermoFileName != fileName::null)
 799     {
 800         std::ifstream thermoStream(thermoFileName.c_str());
 801
 802         if (!thermoStream)
 803         {
 804             FatalErrorIn
 805             (
 806                 "chemkin::chemkin(const fileName& CHEMKINFileName, "
 807                 "const fileName& thermoFileName)"
 808             )   << "file " << thermoFileName << " not found"
 809                 << exit(FatalError);
 810         }
 811
 812         yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
 813         yy_switch_to_buffer(bufferPtr);
 814
 815         while(lex() != 0)
 816         {}
 817
 818         yy_delete_buffer(bufferPtr);
 819
 820         lineNo_ = 1;
 821     }
 822
 823     std::ifstream CHEMKINStream(CHEMKINFileName.c_str());
 824
 825     if (!CHEMKINStream)
 826     {
 827         FatalErrorIn
 828         (
 829             "chemkin::chemkin(const fileName& CHEMKINFileName, "
 830             "const fileName& thermoFileName)"
 831         )   << "file " << CHEMKINFileName << " not found"
 832             << exit(FatalError);
 833     }
 834
 835     yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
 836     yy_switch_to_buffer(bufferPtr);
 837
 838     initReactionKeywordTable();
 839
 840     while(lex() != 0)
 841     {}
 842
 843     yy_delete_buffer(bufferPtr);
 844 }
 845
 846
 847 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
 848
 849 // Construct from components
 850 Foam::chemkinReader::chemkinReader
 851 (
 852     const fileName& CHEMKINFileName,
 853     const fileName& thermoFileName
 854 )
 855 :
 856     lineNo_(1),
 857     specieNames_(10),
 858     speciesTable_(static_cast<const wordList&>(wordList()))
 859 {
 860     read(CHEMKINFileName, thermoFileName);
 861 }
 862
 863
 864 // Construct from components
 865 Foam::chemkinReader::chemkinReader(const dictionary& thermoDict)
 866 :
 867     lineNo_(1),
 868     specieNames_(10),
 869     speciesTable_(static_cast<const wordList&>(wordList()))
 870 {
 871     fileName chemkinFile
 872     (
 873         fileName(thermoDict.lookup("CHEMKINFile")).expand()
 874     );
 875
 876     fileName thermoFile = fileName::null;
 877
 878     if (thermoDict.found("CHEMKINThermoFile"))
 879     {
 880         thermoFile = fileName(thermoDict.lookup("CHEMKINThermoFile")).expand();
 881     }
 882
 883     // allow relative file names
 884     fileName relPath = thermoDict.name().path();
 885     if (relPath.size())
 886     {
 887         if (chemkinFile.size() && chemkinFile[0] != '/')
 888         {
 889             chemkinFile = relPath/chemkinFile;
 890         }
 891
 892         if (thermoFile.size() && thermoFile[0] != '/')
 893         {
 894             thermoFile = relPath/thermoFile;
 895         }
 896     }
 897
 898     read(chemkinFile, thermoFile);
 899 }
 900
 901 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //