FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

chemkinReader.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include <reactionThermophysicalModels/chemkinReader.H>
00027 #include <fstream>
00028 #include <specie/atomicWeights.H>
00029 #include <specie/IrreversibleReaction.H>
00030 #include <specie/ReversibleReaction.H>
00031 #include <specie/NonEquilibriumReversibleReaction.H>
00032 #include <specie/ArrheniusReactionRate.H>
00033 #include <specie/thirdBodyArrheniusReactionRate.H>
00034 #include <specie/FallOffReactionRate.H>
00035 #include <specie/ChemicallyActivatedReactionRate.H>
00036 #include <specie/LindemannFallOffFunction.H>
00037 #include <specie/TroeFallOffFunction.H>
00038 #include <specie/SRIFallOffFunction.H>
00039 #include <specie/LandauTellerReactionRate.H>
00040 #include <specie/JanevReactionRate.H>
00041 #include <specie/powerSeriesReactionRate.H>
00042 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00043 
00044 /* * * * * * * * * * * * * * * * * Static data * * * * * * * * * * * * * * * */
00045 
00046 namespace Foam
00047 {
00048     addChemistryReaderType(chemkinReader, gasThermoPhysics);
00049 };
00050 
00051 /* * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * */
00052 
00053 const char* Foam::chemkinReader::reactionTypeNames[4] =
00054 {
00055     "irreversible",
00056     "reversible",
00057     "nonEquilibriumReversible",
00058     "unknownReactionType"
00059 };
00060 
00061 const char* Foam::chemkinReader::reactionRateTypeNames[8] =
00062 {
00063     "Arrhenius",
00064     "thirdBodyArrhenius",
00065     "unimolecularFallOff",
00066     "chemicallyActivatedBimolecular",
00067     "LandauTeller",
00068     "Janev",
00069     "powerSeries",
00070     "unknownReactionRateType"
00071 };
00072 
00073 const char* Foam::chemkinReader::fallOffFunctionNames[4] =
00074 {
00075     "Lindemann",
00076     "Troe",
00077     "SRI",
00078     "unknownFallOffFunctionType"
00079 };
00080 
00081 void Foam::chemkinReader::initReactionKeywordTable()
00082 {
00083     reactionKeywordTable_.insert("M", thirdBodyReactionType);
00084     reactionKeywordTable_.insert("LOW", unimolecularFallOffReactionType);
00085     reactionKeywordTable_.insert("HIGH", chemicallyActivatedBimolecularReactionType);
00086     reactionKeywordTable_.insert("TROE", TroeReactionType);
00087     reactionKeywordTable_.insert("SRI", SRIReactionType);
00088     reactionKeywordTable_.insert("LT", LandauTellerReactionType);
00089     reactionKeywordTable_.insert("RLT", reverseLandauTellerReactionType);
00090     reactionKeywordTable_.insert("JAN", JanevReactionType);
00091     reactionKeywordTable_.insert("FIT1", powerSeriesReactionRateType);
00092     reactionKeywordTable_.insert("HV", radiationActivatedReactionType);
00093     reactionKeywordTable_.insert("TDEP", speciesTempReactionType);
00094     reactionKeywordTable_.insert("EXCI", energyLossReactionType);
00095     reactionKeywordTable_.insert("MOME", plasmaMomentumTransfer);
00096     reactionKeywordTable_.insert("XSMI", collisionCrossSection);
00097     reactionKeywordTable_.insert("REV", nonEquilibriumReversibleReactionType);
00098     reactionKeywordTable_.insert("DUPLICATE", duplicateReactionType);
00099     reactionKeywordTable_.insert("DUP", duplicateReactionType);
00100     reactionKeywordTable_.insert("FORD", speciesOrderForward);
00101     reactionKeywordTable_.insert("RORD", speciesOrderReverse);
00102     reactionKeywordTable_.insert("UNITS", UnitsOfReaction);
00103     reactionKeywordTable_.insert("END", end);
00104 }
00105 
00106 
00107 Foam::scalar Foam::chemkinReader::molecularWeight
00108 (
00109     const List<specieElement>& specieComposition
00110 ) const
00111 {
00112     scalar molWt = 0.0;
00113 
00114     forAll (specieComposition, i)
00115     {
00116         label nAtoms = specieComposition[i].nAtoms;
00117         const word& elementName = specieComposition[i].elementName;
00118 
00119         if (isotopeAtomicWts_.found(elementName))
00120         {
00121             molWt += nAtoms*isotopeAtomicWts_[elementName];
00122         }
00123         else if (atomicWeights.found(elementName))
00124         {
00125             molWt += nAtoms*atomicWeights[elementName];
00126         }
00127         else
00128         {
00129             FatalErrorIn("chemkinReader::lex()")
00130                 << "Unknown element " << elementName
00131                 << " on line " << lineNo_-1 << nl
00132                 << "    specieComposition: " << specieComposition
00133                 << exit(FatalError);
00134         }
00135     }
00136 
00137     return molWt;
00138 }
00139 
00140 
00141 void Foam::chemkinReader::checkCoeffs
00142 (
00143     const scalarList& reactionCoeffs,
00144     const char* reactionRateName,
00145     const label nCoeffs
00146 ) const
00147 {
00148     if (reactionCoeffs.size() != nCoeffs)
00149     {
00150         FatalErrorIn("chemkinReader::checkCoeffs")
00151             << "Wrong number of coefficients for the " << reactionRateName
00152             << " rate expression on line "
00153             << lineNo_-1 << ", should be "
00154             << nCoeffs << " but " << reactionCoeffs.size() << " supplied." << nl
00155             << "Coefficients are "
00156             << reactionCoeffs << nl
00157             << exit(FatalError);
00158     }
00159 }
00160 
00161 template<class ReactionRateType>
00162 void Foam::chemkinReader::addReactionType
00163 (
00164     const reactionType rType,
00165     DynamicList<gasReaction::specieCoeffs>& lhs,
00166     DynamicList<gasReaction::specieCoeffs>& rhs,
00167     const ReactionRateType& rr
00168 )
00169 {
00170     switch (rType)
00171     {
00172         case irreversible:
00173         {
00174             reactions_.append
00175             (
00176                 new IrreversibleReaction<gasThermoPhysics, ReactionRateType>
00177                 (
00178                     Reaction<gasThermoPhysics>
00179                     (
00180                         speciesTable_,
00181                         lhs.shrink(),
00182                         rhs.shrink(),
00183                         speciesThermo_
00184                     ),
00185                     rr
00186                 )
00187             );
00188         }
00189         break;
00190 
00191         case reversible:
00192         {
00193             reactions_.append
00194             (
00195                 new ReversibleReaction<gasThermoPhysics, ReactionRateType>
00196                 (
00197                     Reaction<gasThermoPhysics>
00198                     (
00199                         speciesTable_,
00200                         lhs.shrink(),
00201                         rhs.shrink(),
00202                         speciesThermo_
00203                     ),
00204                     rr
00205                 )
00206             );
00207         }
00208         break;
00209 
00210         default:
00211 
00212             if (rType < 3)
00213             {
00214                 FatalErrorIn("chemkinReader::addReactionType")
00215                     << "Reaction type " << reactionTypeNames[rType]
00216                     << " on line " << lineNo_-1
00217                     << " not handled by this function"
00218                     << exit(FatalError);
00219             }
00220             else
00221             {
00222                 FatalErrorIn("chemkinReader::addReactionType")
00223                     << "Unknown reaction type " << rType
00224                     << " on line " << lineNo_-1
00225                     << exit(FatalError);
00226             }
00227     }
00228 }
00229 
00230 template<template<class, class> class PressureDependencyType>
00231 void Foam::chemkinReader::addPressureDependentReaction
00232 (
00233     const reactionType rType,
00234     const fallOffFunctionType fofType,
00235     DynamicList<gasReaction::specieCoeffs>& lhs,
00236     DynamicList<gasReaction::specieCoeffs>& rhs,
00237     const scalarList& efficiencies,
00238     const scalarList& k0Coeffs,
00239     const scalarList& kInfCoeffs,
00240     const HashTable<scalarList>& reactionCoeffsTable,
00241     const scalar Afactor0,
00242     const scalar AfactorInf,
00243     const scalar RR
00244 )
00245 {
00246     checkCoeffs(k0Coeffs, "k0", 3);
00247     checkCoeffs(kInfCoeffs, "kInf", 3);
00248 
00249     switch (fofType)
00250     {
00251         case Lindemann:
00252         {
00253             addReactionType
00254             (
00255                 rType,
00256                 lhs, rhs,
00257                 PressureDependencyType
00258                     <ArrheniusReactionRate, LindemannFallOffFunction>
00259                 (
00260                     ArrheniusReactionRate
00261                     (
00262                         Afactor0*k0Coeffs[0],
00263                         k0Coeffs[1],
00264                         k0Coeffs[2]/RR
00265                     ),
00266                     ArrheniusReactionRate
00267                     (
00268                         AfactorInf*kInfCoeffs[0],
00269                         kInfCoeffs[1],
00270                         kInfCoeffs[2]/RR
00271                     ),
00272                     LindemannFallOffFunction(),
00273                     thirdBodyEfficiencies(speciesTable_, efficiencies)
00274                 )
00275             );
00276         }
00277         break;
00278 
00279         case Troe:
00280         {
00281             scalarList TroeCoeffs
00282             (
00283                 reactionCoeffsTable[fallOffFunctionNames[fofType]]
00284             );
00285 
00286             if (TroeCoeffs.size() != 4 && TroeCoeffs.size() != 3)
00287             {
00288                 FatalErrorIn("chemkinReader::addPressureDependentReaction")
00289                     << "Wrong number of coefficients for Troe rate expression"
00290                        " on line " << lineNo_-1 << ", should be 3 or 4 but "
00291                     << TroeCoeffs.size() << " supplied." << nl
00292                     << "Coefficients are "
00293                     << TroeCoeffs << nl
00294                     << exit(FatalError);
00295             }
00296 
00297             if (TroeCoeffs.size() == 3)
00298             {
00299                 TroeCoeffs.setSize(4);
00300                 TroeCoeffs[3] = GREAT;
00301             }
00302 
00303             addReactionType
00304             (
00305                 rType,
00306                 lhs, rhs,
00307                 PressureDependencyType
00308                     <ArrheniusReactionRate, TroeFallOffFunction>
00309                 (
00310                     ArrheniusReactionRate
00311                     (
00312                         Afactor0*k0Coeffs[0],
00313                         k0Coeffs[1],
00314                         k0Coeffs[2]/RR
00315                     ),
00316                     ArrheniusReactionRate
00317                     (
00318                         AfactorInf*kInfCoeffs[0],
00319                         kInfCoeffs[1],
00320                         kInfCoeffs[2]/RR
00321                     ),
00322                     TroeFallOffFunction
00323                     (
00324                         TroeCoeffs[0],
00325                         TroeCoeffs[1],
00326                         TroeCoeffs[2],
00327                         TroeCoeffs[3]
00328                     ),
00329                     thirdBodyEfficiencies(speciesTable_, efficiencies)
00330                 )
00331             );
00332         }
00333         break;
00334 
00335         case SRI:
00336         {
00337             scalarList SRICoeffs
00338             (
00339                 reactionCoeffsTable[fallOffFunctionNames[fofType]]
00340             );
00341 
00342             if (SRICoeffs.size() != 5 && SRICoeffs.size() != 3)
00343             {
00344                 FatalErrorIn("chemkinReader::addPressureDependentReaction")
00345                     << "Wrong number of coefficients for SRI rate expression"
00346                        " on line " << lineNo_-1 << ", should be 3 or 5 but "
00347                     << SRICoeffs.size() << " supplied." << nl
00348                     << "Coefficients are "
00349                     << SRICoeffs << nl
00350                     << exit(FatalError);
00351             }
00352 
00353             if (SRICoeffs.size() == 3)
00354             {
00355                 SRICoeffs.setSize(5);
00356                 SRICoeffs[3] = 1.0;
00357                 SRICoeffs[4] = 0.0;
00358             }
00359 
00360             addReactionType
00361             (
00362                 rType,
00363                 lhs, rhs,
00364                 PressureDependencyType
00365                     <ArrheniusReactionRate, SRIFallOffFunction>
00366                 (
00367                     ArrheniusReactionRate
00368                     (
00369                         Afactor0*k0Coeffs[0],
00370                         k0Coeffs[1],
00371                         k0Coeffs[2]/RR
00372                     ),
00373                     ArrheniusReactionRate
00374                     (
00375                         AfactorInf*kInfCoeffs[0],
00376                         kInfCoeffs[1],
00377                         kInfCoeffs[2]/RR
00378                     ),
00379                     SRIFallOffFunction
00380                     (
00381                         SRICoeffs[0],
00382                         SRICoeffs[1],
00383                         SRICoeffs[2],
00384                         SRICoeffs[3],
00385                         SRICoeffs[4]
00386                     ),
00387                     thirdBodyEfficiencies(speciesTable_, efficiencies)
00388                 )
00389             );
00390         }
00391         break;
00392 
00393         default:
00394         {
00395             if (fofType < 4)
00396             {
00397                 FatalErrorIn("chemkinReader::addPressureDependentReaction")
00398                     << "Fall-off function type "
00399                     << fallOffFunctionNames[fofType]
00400                     << " on line " << lineNo_-1
00401                     << " not implemented"
00402                     << exit(FatalError);
00403             }
00404             else
00405             {
00406                 FatalErrorIn("chemkinReader::addPressureDependentReaction")
00407                     << "Unknown fall-off function type " << fofType
00408                     << " on line " << lineNo_-1
00409                     << exit(FatalError);
00410             }
00411         }
00412     }
00413 }
00414 
00415 
00416 void Foam::chemkinReader::addReaction
00417 (
00418     DynamicList<gasReaction::specieCoeffs>& lhs,
00419     DynamicList<gasReaction::specieCoeffs>& rhs,
00420     const scalarList& efficiencies,
00421     const reactionType rType,
00422     const reactionRateType rrType,
00423     const fallOffFunctionType fofType,
00424     const scalarList& ArrheniusCoeffs,
00425     HashTable<scalarList>& reactionCoeffsTable,
00426     const scalar RR
00427 )
00428 {
00429     checkCoeffs(ArrheniusCoeffs, "Arrhenius", 3);
00430 
00431     scalarList nAtoms(elementNames_.size(), 0.0);
00432 
00433     forAll (lhs, i)
00434     {
00435         const List<specieElement>& specieComposition =
00436             specieComposition_[speciesTable_[lhs[i].index]];
00437 
00438         forAll (specieComposition, j)
00439         {
00440             label elementi = elementIndices_[specieComposition[j].elementName];
00441             nAtoms[elementi] += lhs[i].stoichCoeff*specieComposition[j].nAtoms;
00442         }
00443     }
00444 
00445     forAll (rhs, i)
00446     {
00447         const List<specieElement>& specieComposition =
00448             specieComposition_[speciesTable_[rhs[i].index]];
00449 
00450         forAll (specieComposition, j)
00451         {
00452             label elementi = elementIndices_[specieComposition[j].elementName];
00453             nAtoms[elementi] -= rhs[i].stoichCoeff*specieComposition[j].nAtoms;
00454         }
00455     }
00456 
00457 
00458     // Calculate the unit conversion factor for the A coefficient
00459     // for the change from mol/cm^3 to kmol/m^3 concentraction units
00460     const scalar concFactor = 0.001;
00461     scalar sumExp = 0.0;
00462     forAll (lhs, i)
00463     {
00464         sumExp += lhs[i].exponent;
00465     }
00466     scalar Afactor = pow(concFactor, sumExp - 1.0);
00467 
00468     scalar AfactorRev = Afactor;
00469 
00470     if (rType == nonEquilibriumReversible)
00471     {
00472         sumExp = 0.0;
00473         forAll (rhs, i)
00474         {
00475             sumExp += rhs[i].exponent;
00476         }
00477         AfactorRev = pow(concFactor, sumExp - 1.0);
00478     }
00479 
00480     switch (rrType)
00481     {
00482         case Arrhenius:
00483         {
00484             if (rType == nonEquilibriumReversible)
00485             {
00486                 const scalarList& reverseArrheniusCoeffs =
00487                     reactionCoeffsTable[reactionTypeNames[rType]];
00488 
00489                 checkCoeffs(reverseArrheniusCoeffs, "reverse Arrhenius", 3);
00490 
00491                 reactions_.append
00492                 (
00493                     new NonEquilibriumReversibleReaction
00494                         <gasThermoPhysics, ArrheniusReactionRate>
00495                     (
00496                         Reaction<gasThermoPhysics>
00497                         (
00498                             speciesTable_,
00499                             lhs.shrink(),
00500                             rhs.shrink(),
00501                             speciesThermo_
00502                         ),
00503                         ArrheniusReactionRate
00504                         (
00505                             Afactor*ArrheniusCoeffs[0],
00506                             ArrheniusCoeffs[1],
00507                             ArrheniusCoeffs[2]/RR
00508                         ),
00509                         ArrheniusReactionRate
00510                         (
00511                             AfactorRev*reverseArrheniusCoeffs[0],
00512                             reverseArrheniusCoeffs[1],
00513                             reverseArrheniusCoeffs[2]/RR
00514                         )
00515                     )
00516                 );
00517             }
00518             else
00519             {
00520                 addReactionType
00521                 (
00522                     rType,
00523                     lhs, rhs,
00524                     ArrheniusReactionRate
00525                     (
00526                         Afactor*ArrheniusCoeffs[0],
00527                         ArrheniusCoeffs[1],
00528                         ArrheniusCoeffs[2]/RR
00529                     )
00530                 );
00531             }
00532         }
00533         break;
00534 
00535         case thirdBodyArrhenius:
00536         {
00537             if (rType == nonEquilibriumReversible)
00538             {
00539                 const scalarList& reverseArrheniusCoeffs =
00540                     reactionCoeffsTable[reactionTypeNames[rType]];
00541 
00542                 checkCoeffs(reverseArrheniusCoeffs, "reverse Arrhenius", 3);
00543 
00544                 reactions_.append
00545                 (
00546                     new NonEquilibriumReversibleReaction
00547                         <gasThermoPhysics, thirdBodyArrheniusReactionRate>
00548                     (
00549                         Reaction<gasThermoPhysics>
00550                         (
00551                             speciesTable_,
00552                             lhs.shrink(),
00553                             rhs.shrink(),
00554                             speciesThermo_
00555                         ),
00556                         thirdBodyArrheniusReactionRate
00557                         (
00558                             Afactor*concFactor*ArrheniusCoeffs[0],
00559                             ArrheniusCoeffs[1],
00560                             ArrheniusCoeffs[2]/RR,
00561                             thirdBodyEfficiencies(speciesTable_, efficiencies)
00562                         ),
00563                         thirdBodyArrheniusReactionRate
00564                         (
00565                             AfactorRev*concFactor*reverseArrheniusCoeffs[0],
00566                             reverseArrheniusCoeffs[1],
00567                             reverseArrheniusCoeffs[2]/RR,
00568                             thirdBodyEfficiencies(speciesTable_, efficiencies)
00569                         )
00570                     )
00571                 );
00572             }
00573             else
00574             {
00575                 addReactionType
00576                 (
00577                     rType,
00578                     lhs, rhs,
00579                     thirdBodyArrheniusReactionRate
00580                     (
00581                         Afactor*concFactor*ArrheniusCoeffs[0],
00582                         ArrheniusCoeffs[1],
00583                         ArrheniusCoeffs[2]/RR,
00584                         thirdBodyEfficiencies(speciesTable_, efficiencies)
00585                     )
00586                 );
00587             }
00588         }
00589         break;
00590 
00591         case unimolecularFallOff:
00592         {
00593             addPressureDependentReaction<FallOffReactionRate>
00594             (
00595                 rType,
00596                 fofType,
00597                 lhs,
00598                 rhs,
00599                 efficiencies,
00600                 reactionCoeffsTable[reactionRateTypeNames[rrType]],
00601                 ArrheniusCoeffs,
00602                 reactionCoeffsTable,
00603                 concFactor*Afactor,
00604                 Afactor,
00605                 RR
00606             );
00607         }
00608         break;
00609 
00610         case chemicallyActivatedBimolecular:
00611         {
00612             addPressureDependentReaction<ChemicallyActivatedReactionRate>
00613             (
00614                 rType,
00615                 fofType,
00616                 lhs,
00617                 rhs,
00618                 efficiencies,
00619                 ArrheniusCoeffs,
00620                 reactionCoeffsTable[reactionRateTypeNames[rrType]],
00621                 reactionCoeffsTable,
00622                 Afactor,
00623                 Afactor/concFactor,
00624                 RR
00625             );
00626         }
00627         break;
00628 
00629         case LandauTeller:
00630         {
00631             const scalarList& LandauTellerCoeffs =
00632                 reactionCoeffsTable[reactionRateTypeNames[rrType]];
00633             checkCoeffs(LandauTellerCoeffs, "Landau-Teller", 2);
00634 
00635             if (rType == nonEquilibriumReversible)
00636             {
00637                 const scalarList& reverseArrheniusCoeffs =
00638                     reactionCoeffsTable[reactionTypeNames[rType]];
00639                 checkCoeffs(reverseArrheniusCoeffs, "reverse Arrhenius", 3);
00640 
00641                 const scalarList& reverseLandauTellerCoeffs =
00642                     reactionCoeffsTable
00643                     [
00644                         word(reactionTypeNames[rType])
00645                       + reactionRateTypeNames[rrType]
00646                     ];
00647                 checkCoeffs(LandauTellerCoeffs, "reverse Landau-Teller", 2);
00648 
00649                 reactions_.append
00650                 (
00651                     new NonEquilibriumReversibleReaction
00652                         <gasThermoPhysics, LandauTellerReactionRate>
00653                     (
00654                         Reaction<gasThermoPhysics>
00655                         (
00656                             speciesTable_,
00657                             lhs.shrink(),
00658                             rhs.shrink(),
00659                             speciesThermo_
00660                         ),
00661                         LandauTellerReactionRate
00662                         (
00663                             Afactor*ArrheniusCoeffs[0],
00664                             ArrheniusCoeffs[1],
00665                             ArrheniusCoeffs[2]/RR,
00666                             LandauTellerCoeffs[0],
00667                             LandauTellerCoeffs[1]
00668                         ),
00669                         LandauTellerReactionRate
00670                         (
00671                             AfactorRev*reverseArrheniusCoeffs[0],
00672                             reverseArrheniusCoeffs[1],
00673                             reverseArrheniusCoeffs[2]/RR,
00674                             reverseLandauTellerCoeffs[0],
00675                             reverseLandauTellerCoeffs[1]
00676                         )
00677                     )
00678                 );
00679             }
00680             else
00681             {
00682                 addReactionType
00683                 (
00684                     rType,
00685                     lhs, rhs,
00686                     LandauTellerReactionRate
00687                     (
00688                         Afactor*ArrheniusCoeffs[0],
00689                         ArrheniusCoeffs[1],
00690                         ArrheniusCoeffs[2]/RR,
00691                         LandauTellerCoeffs[0],
00692                         LandauTellerCoeffs[1]
00693                     )
00694                 );
00695             }
00696         }
00697         break;
00698 
00699         case Janev:
00700         {
00701             const scalarList& JanevCoeffs =
00702                 reactionCoeffsTable[reactionRateTypeNames[rrType]];
00703 
00704             checkCoeffs(JanevCoeffs, "Janev", 9);
00705 
00706             addReactionType
00707             (
00708                 rType,
00709                 lhs, rhs,
00710                 JanevReactionRate
00711                 (
00712                     Afactor*ArrheniusCoeffs[0],
00713                     ArrheniusCoeffs[1],
00714                     ArrheniusCoeffs[2]/RR,
00715                     JanevCoeffs.begin()
00716                 )
00717             );
00718         }
00719         break;
00720 
00721         case powerSeries:
00722         {
00723             const scalarList& powerSeriesCoeffs =
00724                 reactionCoeffsTable[reactionRateTypeNames[rrType]];
00725 
00726             checkCoeffs(powerSeriesCoeffs, "power-series", 4);
00727 
00728             addReactionType
00729             (
00730                 rType,
00731                 lhs, rhs,
00732                 powerSeriesReactionRate
00733                 (
00734                     Afactor*ArrheniusCoeffs[0],
00735                     ArrheniusCoeffs[1],
00736                     ArrheniusCoeffs[2]/RR,
00737                     powerSeriesCoeffs.begin()
00738                 )
00739             );
00740         }
00741         break;
00742 
00743         case unknownReactionRateType:
00744         {
00745             FatalErrorIn("chemkinReader::addReaction")
00746                 << "Internal error on line " << lineNo_-1
00747                 << ": reaction rate type has not been set"
00748                 << exit(FatalError);
00749         }
00750         break;
00751 
00752         default:
00753         {
00754             if (rrType < 9)
00755             {
00756                 FatalErrorIn("chemkinReader::addReaction")
00757                     << "Reaction rate type " << reactionRateTypeNames[rrType]
00758                     << " on line " << lineNo_-1
00759                     << " not implemented"
00760                     << exit(FatalError);
00761             }
00762             else
00763             {
00764                 FatalErrorIn("chemkinReader::addReaction")
00765                     << "Unknown reaction rate type " << rrType
00766                     << " on line " << lineNo_-1
00767                     << exit(FatalError);
00768             }
00769         }
00770     }
00771 
00772 
00773     forAll (nAtoms, i)
00774     {
00775         if (mag(nAtoms[i]) > SMALL)
00776         {
00777             FatalErrorIn("chemkinReader::addReaction")
00778                 << "Elemental imbalance in " << elementNames_[i]
00779                 << " in reaction" << nl
00780                 << reactions_.last() << nl
00781                 << " on line " << lineNo_-1
00782                 << exit(FatalError);
00783         }
00784     }
00785 
00786     lhs.clear();
00787     rhs.clear();
00788     reactionCoeffsTable.clear();
00789 }
00790 
00791 
00792 void Foam::chemkinReader::read
00793 (
00794     const fileName& CHEMKINFileName,
00795     const fileName& thermoFileName
00796 )
00797 {
00798     if (thermoFileName != fileName::null)
00799     {
00800         std::ifstream thermoStream(thermoFileName.c_str());
00801 
00802         if (!thermoStream)
00803         {
00804             FatalErrorIn
00805             (
00806                 "chemkin::chemkin(const fileName& CHEMKINFileName, "
00807                 "const fileName& thermoFileName)"
00808             )   << "file " << thermoFileName << " not found"
00809                 << exit(FatalError);
00810         }
00811 
00812         yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
00813         yy_switch_to_buffer(bufferPtr);
00814 
00815         while(lex() != 0)
00816         {}
00817 
00818         yy_delete_buffer(bufferPtr);
00819 
00820         lineNo_ = 1;
00821     }
00822 
00823     std::ifstream CHEMKINStream(CHEMKINFileName.c_str());
00824 
00825     if (!CHEMKINStream)
00826     {
00827         FatalErrorIn
00828         (
00829             "chemkin::chemkin(const fileName& CHEMKINFileName, "
00830             "const fileName& thermoFileName)"
00831         )   << "file " << CHEMKINFileName << " not found"
00832             << exit(FatalError);
00833     }
00834 
00835     yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
00836     yy_switch_to_buffer(bufferPtr);
00837 
00838     initReactionKeywordTable();
00839 
00840     while(lex() != 0)
00841     {}
00842 
00843     yy_delete_buffer(bufferPtr);
00844 }
00845 
00846 
00847 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
00848 
00849 // Construct from components
00850 Foam::chemkinReader::chemkinReader
00851 (
00852     const fileName& CHEMKINFileName,
00853     const fileName& thermoFileName
00854 )
00855 :
00856     lineNo_(1),
00857     specieNames_(10),
00858     speciesTable_(static_cast<const wordList&>(wordList()))
00859 {
00860     read(CHEMKINFileName, thermoFileName);
00861 }
00862 
00863 
00864 // Construct from components
00865 Foam::chemkinReader::chemkinReader(const dictionary& thermoDict)
00866 :
00867     lineNo_(1),
00868     specieNames_(10),
00869     speciesTable_(static_cast<const wordList&>(wordList()))
00870 {
00871     fileName chemkinFile
00872     (
00873         fileName(thermoDict.lookup("CHEMKINFile")).expand()
00874     );
00875 
00876     fileName thermoFile = fileName::null;
00877 
00878     if (thermoDict.found("CHEMKINThermoFile"))
00879     {
00880         thermoFile = fileName(thermoDict.lookup("CHEMKINThermoFile")).expand();
00881     }
00882 
00883     // allow relative file names
00884     fileName relPath = thermoDict.name().path();
00885     if (relPath.size())
00886     {
00887         if (chemkinFile.size() && chemkinFile[0] != '/')
00888         {
00889             chemkinFile = relPath/chemkinFile;
00890         }
00891 
00892         if (thermoFile.size() && thermoFile[0] != '/')
00893         {
00894             thermoFile = relPath/thermoFile;
00895         }
00896     }
00897 
00898     read(chemkinFile, thermoFile);
00899 }
00900 
00901 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines