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

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