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