foamChemistryReader: Added support for elements and specie composition

Based on a patch contributed by Francesco Contino, Tommaso Lucchini,
Gianluca D’Errico, Hervé Jeanmart, Nicolas Bourgeois and Stéphane
Backaert.

applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C

@@ -88,8 +88,8 @@ int main(int argc, char *argv[])
         {
             elementsDict.add
             (
-                cr.specieComposition()[speciesList[si]][ei].elementName,
+                cr.specieComposition()[speciesList[si]][ei].name(),
-                cr.specieComposition()[speciesList[si]][ei].nAtoms
+                cr.specieComposition()[speciesList[si]][ei].nAtoms()
             );
         }
 

src/thermophysicalModels/reactionThermo/chemistryReaders/chemistryReader/chemistryReader.H

@@ -36,15 +36,19 @@ SourceFiles
 #define chemistryReader_H
 
 #include "typeInfo.H"
-#include "runTimeSelectionTables.H"
+#include "specieElement.H"
 #include "Reaction.H"
 #include "ReactionList.H"
+#include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+typedef HashTable<List<specieElement>> speciesCompositionTable;
+
+
 /*---------------------------------------------------------------------------*\
                       Class chemistryReader Declaration
 \*---------------------------------------------------------------------------*/
@@ -112,6 +116,9 @@ public:
         //- Return access to the list of species
         virtual const speciesTable& species() const = 0;
 
+        //- Table of species composition
+        virtual const speciesCompositionTable& specieComposition() const = 0;
+
         //- Return access to the thermo packages
         virtual const HashPtrTable<ThermoType>& speciesThermo() const = 0;
 

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L

@@ -68,9 +68,6 @@ int yyFlexLexer::yywrap()
 Foam::string foamSpecieString(const char* YYText)
 {
     Foam::string specieString(YYText);
-    // Transforming parentheses is no longer necessary
-    //specieString.replaceAll('(', '<');
-    //specieString.replaceAll(')', '>');
     return specieString;
 }
 
@@ -483,9 +480,9 @@ bool finishReaction = false;
             if (elementName.size() && nAtoms)
             {
                 correctElementName(elementName);
-                currentSpecieComposition[nSpecieElements].elementName =
+                currentSpecieComposition[nSpecieElements].name() =
                     elementName;
-                currentSpecieComposition[nSpecieElements++].nAtoms = nAtoms;
+                currentSpecieComposition[nSpecieElements++].nAtoms() = nAtoms;
             }
         }
 
@@ -540,24 +537,24 @@ bool finishReaction = false;
         )
         {
             correctElementName(elementName);
-            currentSpecieComposition[nSpecieElements].elementName =
+            currentSpecieComposition[nSpecieElements].name() =
                 elementName;
-            currentSpecieComposition[nSpecieElements++].nAtoms = nAtoms;
+            currentSpecieComposition[nSpecieElements++].nAtoms() = nAtoms;
         }
 
         currentSpecieComposition.setSize(nSpecieElements);
 
-        HashTable<List<specieElement>>::iterator specieCompositionIter
+        speciesCompositionTable::iterator specieCompositionIter
         (
-            specieComposition_.find(currentSpecieName)
+            speciesComposition_.find(currentSpecieName)
         );
 
-        if (specieCompositionIter != specieComposition_.end())
+        if (specieCompositionIter != speciesComposition_.end())
         {
-            specieComposition_.erase(specieCompositionIter);
+            speciesComposition_.erase(specieCompositionIter);
         }
 
-        specieComposition_.insert
+        speciesComposition_.insert
         (
             currentSpecieName,
             currentSpecieComposition

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

@@ -119,8 +119,8 @@ Foam::scalar Foam::chemkinReader::molecularWeight
 
     forAll(specieComposition, i)
     {
-        label nAtoms = specieComposition[i].nAtoms;
+        label nAtoms = specieComposition[i].nAtoms();
-        const word& elementName = specieComposition[i].elementName;
+        const word& elementName = specieComposition[i].name();
 
         if (isotopeAtomicWts_.found(elementName))
         {
@@ -428,24 +428,26 @@ void Foam::chemkinReader::addReaction
     forAll(lhs, i)
     {
         const List<specieElement>& specieComposition =
-            specieComposition_[speciesTable_[lhs[i].index]];
+            speciesComposition_[speciesTable_[lhs[i].index]];
 
         forAll(specieComposition, j)
         {
-            label elementi = elementIndices_[specieComposition[j].elementName];
+            label elementi = elementIndices_[specieComposition[j].name()];
-            nAtoms[elementi] += lhs[i].stoichCoeff*specieComposition[j].nAtoms;
+            nAtoms[elementi] +=
+                lhs[i].stoichCoeff*specieComposition[j].nAtoms();
         }
     }
 
     forAll(rhs, i)
     {
         const List<specieElement>& specieComposition =
-            specieComposition_[speciesTable_[rhs[i].index]];
+            speciesComposition_[speciesTable_[rhs[i].index]];
 
         forAll(specieComposition, j)
         {
-            label elementi = elementIndices_[specieComposition[j].elementName];
+            label elementi = elementIndices_[specieComposition[j].name()];
-            nAtoms[elementi] -= rhs[i].stoichCoeff*specieComposition[j].nAtoms;
+            nAtoms[elementi] -=
+                rhs[i].stoichCoeff*specieComposition[j].nAtoms();
         }
     }
 

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.H

@@ -46,7 +46,6 @@ SourceFiles
 #include "labelList.H"
 #include "speciesTable.H"
 #include "atomicWeights.H"
-
 #include "reactionTypes.H"
 
 #include <FlexLexer.h>
@@ -77,33 +76,6 @@ public:
             gas
         };
 
-        //- species element
-        struct specieElement
-        {
-            word elementName;
-            label nAtoms;
-
-            bool operator==(const specieElement& se) const
-            {
-                return
-                (
-                    nAtoms == se.nAtoms
-                 && elementName == se.elementName
-                );
-            }
-
-            bool operator!=(const specieElement& se) const
-            {
-                return !operator==(se);
-            }
-
-            friend Ostream& operator<<(Ostream& os, const specieElement& se)
-            {
-                os  << se.nAtoms << token::SPACE << se.elementName;
-                return os;
-            }
-        };
-
 
 private:
 
@@ -203,7 +175,7 @@ private:
         HashPtrTable<gasHThermoPhysics> speciesThermo_;
 
         //- Table of species composition
-        HashTable<List<specieElement>> specieComposition_;
+        speciesCompositionTable speciesComposition_;
 
         //- List of the reactions
         ReactionList<gasHThermoPhysics> reactions_;
@@ -369,6 +341,12 @@ public:
             return speciesTable_;
         }
 
+        //- Table of species composition
+        const speciesCompositionTable& specieComposition() const
+        {
+            return speciesComposition_;
+        }
+
         //- Specie phase
         const HashTable<phase>& speciePhase() const
         {
@@ -381,12 +359,6 @@ public:
             return speciesThermo_;
         }
 
-        //- Table of species composition
-        const HashTable<List<specieElement>>& specieComposition() const
-        {
-            return specieComposition_;
-        }
-
         //- List of the reactions
         const ReactionList<gasHThermoPhysics>& reactions() const
         {

src/thermophysicalModels/reactionThermo/chemistryReaders/foamChemistryReader/foamChemistryReader.C

@@ -42,6 +42,85 @@ Foam::speciesTable& Foam::foamChemistryReader<ThermoType>::setSpecies
 }
 
 
+template<class ThermoType>
+void Foam::foamChemistryReader<ThermoType>::readSpeciesComposition()
+{
+    if (!chemDict_.found("elements"))
+    {
+        Info<< "    elements not defined in " << chemDict_.name() << endl;
+        return;
+    }
+
+    wordList e(chemDict_.lookup("elements"));
+    label currentElementIndex(0);
+
+    DynamicList<word> elementNames_;
+    HashTable<label> elementIndices_;
+
+    forAll(e, ei)
+    {
+        if (!elementIndices_.found(e[ei]))
+        {
+            elementIndices_.insert(e[ei], currentElementIndex++);
+            elementNames_.append(e[ei]);
+        }
+        else
+        {
+            IOWarningInFunction(chemDict_)
+                << "element " << e[ei] << " already in table." << endl;
+        }
+    }
+
+    // Loop through all species in thermoDict to retrieve
+    // the species composition
+    forAll(speciesTable_, si)
+    {
+        if (thermoDict_.subDict(speciesTable_[si]).isDict("elements"))
+        {
+            dictionary currentElements
+            (
+                thermoDict_.subDict(speciesTable_[si]).subDict("elements")
+            );
+
+            wordList currentElementsName(currentElements.toc());
+            List<specieElement> currentComposition(currentElementsName.size());
+
+            forAll(currentElementsName, eni)
+            {
+                currentComposition[eni].name() = currentElementsName[eni];
+
+                currentComposition[eni].nAtoms() =
+                    currentElements.lookupOrDefault
+                    (
+                        currentElementsName[eni],
+                        0
+                    );
+            }
+
+            // Add current specie composition to the hash table
+            speciesCompositionTable::iterator specieCompositionIter
+            (
+                speciesComposition_.find(speciesTable_[si])
+            );
+
+            if (specieCompositionIter != speciesComposition_.end())
+            {
+                speciesComposition_.erase(specieCompositionIter);
+            }
+
+            speciesComposition_.insert(speciesTable_[si], currentComposition);
+        }
+        else
+        {
+            FatalIOErrorInFunction(thermoDict_)
+                << "Specie " << speciesTable_[si]
+                << " does not contain element description."
+                << exit(FatalIOError);
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
 
 template<class ThermoType>
@@ -70,7 +149,9 @@ Foam::foamChemistryReader<ThermoType>::foamChemistryReader
     speciesTable_(setSpecies(chemDict_, species)),
     speciesThermo_(thermoDict_),
     reactions_(speciesTable_, speciesThermo_, chemDict_)
-{}
+{
+    readSpeciesComposition();
+}
 
 
 template<class ThermoType>
@@ -95,10 +176,12 @@ Foam::foamChemistryReader<ThermoType>::foamChemistryReader
             fileName(thermoDict.lookup("foamChemistryThermoFile")).expand()
         )()
     ),
-    speciesThermo_(thermoDict_),
     speciesTable_(setSpecies(chemDict_, species)),
+    speciesThermo_(thermoDict_),
     reactions_(speciesTable_, speciesThermo_, chemDict_)
-{}
+{
+    readSpeciesComposition();
+}
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/reactionThermo/chemistryReaders/foamChemistryReader/foamChemistryReader.H

@@ -63,12 +63,21 @@ class foamChemistryReader
         //- Thermo properties dictionary
         dictionary thermoDict_;
 
-        //- Table of the thermodynamic data given in the foamChemistry file
+        //- List of elements
-        HashPtrTable<ThermoType> speciesThermo_;
+        DynamicList<word> elementNames_;
+
+        //- Element indices
+        HashTable<label> elementIndices_;
 
         //- Table of species
         speciesTable& speciesTable_;
 
+        //- Table of species composition
+        speciesCompositionTable speciesComposition_;
+
+        //- Table of the thermodynamic data given in the foamChemistry file
+        HashPtrTable<ThermoType> speciesThermo_;
+
         //- List of the reactions
         ReactionList<ThermoType> reactions_;
 
@@ -78,6 +87,9 @@ class foamChemistryReader
         //- Set the species list
         speciesTable& setSpecies(const dictionary& dict, speciesTable& species);
 
+        //- Read the species composition
+        void readSpeciesComposition();
+
         //- Disallow default bitwise copy construct
         foamChemistryReader(const foamChemistryReader&);
 
@@ -117,12 +129,30 @@ public:
 
     // Member functions
 
+        //- List of elements
+        const wordList& elementNames() const
+        {
+            return elementNames_;
+        }
+
+        //- Element indices
+        const HashTable<label>& elementIndices() const
+        {
+            return elementIndices_;
+        }
+
         //- Table of species
         const speciesTable& species() const
         {
             return speciesTable_;
         }
 
+        //- Table of species composition
+        const speciesCompositionTable& specieComposition() const
+        {
+            return speciesComposition_;
+        }
+
         //- Table of the thermodynamic data given in the foamChemistry file
         const HashPtrTable<ThermoType>& speciesThermo() const
         {

src/thermophysicalModels/specie/specieElement/specieElement.H

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::specieElement
+
+Description
+
+SourceFiles
+    specieElementI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef specieElement_H
+#define specieElement_H
+
+#include "word.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+class specieElement;
+Ostream& operator<<(Ostream&, const specieElement&);
+
+
+/*---------------------------------------------------------------------------*\
+                         Class specieElement Declaration
+\*---------------------------------------------------------------------------*/
+
+class specieElement
+{
+    // Private data
+
+        //- Name of the element
+        word name_;
+
+        //- Number of atoms of this element in the specie
+        label nAtoms_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline specieElement();
+
+        //- Construct from components
+        inline specieElement(const word& name, const label nAtoms);
+
+        //- Construct from Istream
+        inline specieElement(Istream&);
+
+
+    // Member Functions
+
+        //- Return the name of the element
+        inline const word& name() const;
+
+        //- Return non-const access to the name of the element
+        inline word& name();
+
+        //- Return the number of atoms of this element in the specie
+        inline label nAtoms() const;
+
+        //- Return non-const access to the number of atoms of this element
+        //  in the specie
+        inline label& nAtoms();
+
+
+    // Member Operators
+
+        //- Equality comparison
+        inline bool operator==(const specieElement&) const;
+
+        //- Inequality comparison
+        inline bool operator!=(const specieElement&) const;
+
+
+    // IOstream Operators
+
+        inline friend Ostream& operator<<(Ostream&, const specieElement&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "specieElementI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/specie/specieElement/specieElementI.H

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "token.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::specieElement::specieElement()
+{}
+
+
+inline Foam::specieElement::specieElement(const word& name, const label nAtoms)
+:
+    name_(name),
+    nAtoms_(nAtoms)
+{}
+
+
+inline Foam::specieElement::specieElement(Istream& is)
+:
+    name_(is),
+    nAtoms_(readLabel(is))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::word& Foam::specieElement::name() const
+{
+    return name_;
+}
+
+
+inline Foam::word& Foam::specieElement::name()
+{
+    return name_;
+}
+
+
+inline Foam::label Foam::specieElement::nAtoms() const
+{
+    return nAtoms_;
+}
+
+
+inline Foam::label& Foam::specieElement::nAtoms()
+{
+    return nAtoms_;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::specieElement::operator==(const specieElement& se) const
+{
+    return
+    (
+        nAtoms_ == se.nAtoms_
+     && name_ == se.name_
+    );
+}
+
+
+inline bool Foam::specieElement::operator!=(const specieElement& se) const
+{
+    return !operator==(se);
+}
+
+
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+inline Foam::Ostream& Foam::operator<<(Ostream& os, const specieElement& se)
+{
+    os  << se.name() << token::SPACE << se.nAtoms();
+    return os;
+}
+
+
+// ************************************************************************* //

tutorials/combustion/chemFoam/gri/Allclean

@@ -6,7 +6,7 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 cleanCase
 
-rm -rf 0 chemFoam.out validation/OF_vs_CHEMKINII.eps validation/chemkinII
+rm -rf 0 chemFoam.out constant/reactions constant/thermo \
+   validation/OF_vs_CHEMKINII.eps validation/chemkinII
 
 #------------------------------------------------------------------------------
-

tutorials/combustion/chemFoam/gri/Allrun

@@ -7,6 +7,10 @@ cd ${0%/*} || exit 1    # Run from this directory
 # Set application name
 application=`getApplication`
 
+runApplication chemkinToFoam \
+               chemkin/chem.inp chemkin/therm.dat chemkin/transportProperties \
+               constant/reactions constant/thermo
+
 runApplication $application
 
 (cd validation && ./Allrun $*)

tutorials/combustion/chemFoam/gri/constant/thermophysicalProperties

@@ -26,8 +26,8 @@ thermoType
     specie          specie;
 }
 
-CHEMKINFile     "$FOAM_CASE/chemkin/chem.inp";
+chemistryReader foamChemistryReader;
-CHEMKINThermoFile "$FOAM_CASE/chemkin/therm.dat";
+foamChemistryFile "$FOAM_CASE/constant/reactions";
-CHEMKINTransportFile "$FOAM_CASE/chemkin/transportProperties";
+foamChemistryThermoFile "$FOAM_CASE/constant/thermo";
 
 // ************************************************************************* //