Update .gitlab-ci.yml

update CI Lint for new GitLab Runner
2018-10-31 15:46:39 +00:00 · 2018-10-31 14:27:04 +00:00 · 2018-10-25 07:31:35 -04:00 · 2018-07-10 15:24:03 +09:00 · 2018-07-10 14:59:50 +09:00 · 2018-04-25 14:56:01 +09:00
192 changed files with 411 additions and 380827 deletions
--- a/.gitignore
+++ b/.gitignore
@ -7,8 +7,6 @@
 *.bak
 *.bak[0-9][0-9]
 \#*\#
 .*.swp
 .*.swo
 # File-browser settings - anywhere
 .directory
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@ -1,9 +1,23 @@
 allmake:
  image: park0d/of4builder
  tags:
    - openfoam4
  before_script:
    - "[[ -d ../ThirdParty-4.x ]] || git clone https://github.com/OpenFOAM/ThirdParty-4.x.git ../ThirdParty-4.x"
    - source etc/bashrc
    - rm -rf platforms
    - wclean all
 allmake:
  script:
    - ./Allwmake -j
  artifacts:
    paths:
      - ./
 cantera:
    image: park0d/of4builder
    tags:
        - openfoam4
    script:
        - whereis libcantera
--- a/applications/solvers/combustion/SLFMFoam/AMC/AMC.C
+++ b/applications/solvers/combustion/SLFMFoam/AMC/AMC.C
@ -1,49 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "AMC.H"
 #include "error.H"
 // * * * * * * * * * * * * * Functions * * * * * * * * * * * * //
 //Amplitude Mapping Closure (from KIVA)
 Foam::tmp<Foam::scalarField> Foam::CMC::AMC(const scalarField& eta)
 {
    tmp<scalarField> tRes(new scalarField(eta.size(), 0.0));
    scalarField& Res = tRes.ref();
    for(label i=0 ; i<eta.size() ; i++)
    {
        Res[i] = AMC(eta[i]);
    }
    return tRes;
 }
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/AMC/AMC.H
+++ b/applications/solvers/combustion/SLFMFoam/AMC/AMC.H
@ -1,74 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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/>.
 Function
    Foam::CMC::AMC
 Description
 SourceFiles
    AMC.C
 \*---------------------------------------------------------------------------*/
 #ifndef AMC_H
 #define AMC_H
 #include "tmp.H"
 #include "scalarField.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 namespace Foam
 {
 namespace CMC
 {
 /*---------------------------------------------------------------------------*\
                         Function AMC Declaration
 \*---------------------------------------------------------------------------*/
    const scalar pi = 3.141592;
    const scalar spi = Foam::sqrt(pi);
    //Amplitude Mapping Closure (from KIVA)
    tmp<scalarField> AMC(const scalarField& eta);
    //Amplitude Mapping Closure (from KIVA)
    //Define exp(-2*(erf^-1(2*eta - 1))^2)
    inline scalar AMC(const scalar eta);
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 } // End namespace CMC
 } // End namespace Foam
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #include "AMCI.H"
 #endif
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/AMC/AMCI.H
+++ b/applications/solvers/combustion/SLFMFoam/AMC/AMCI.H
@ -1,54 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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/>.
 \*---------------------------------------------------------------------------*/
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 //Amplitude Mapping Closure (from KIVA)
 //Define exp(-2*(erf^-1(2*eta - 1))^2)
 inline Foam::scalar Foam::CMC::AMC(const scalar eta)
 {
    const scalar a0 = (2.0*eta - 1.0);
    scalar a1 = 0.5;
    scalar slope;
    scalar da = GREAT;
    while(mag(da) > SMALL)
    {
        slope = (2.0/spi)*Foam::exp(-1.0*Foam::pow(a1,2.0));
        da = (a0 - Foam::erf(a1))/slope;
        a1 = a1 + da;
    }
    return Foam::exp(-2.0*Foam::pow(a1,2.0));
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/BetaFunction/BetaFunction.C
+++ b/applications/solvers/combustion/SLFMFoam/BetaFunction/BetaFunction.C
@ -1,120 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "BetaFunction.H"
 #include "error.H"
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 // const scalar Foam::CMC::BetaFunction::staticData();
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 void Foam::CMC::BetaFunction::_setParameters()
 {
        scalar gamma = mf_*(1.0 - mf_)/(mfVar_ + SMALL) - 1.0;
        if (gamma < SMALL || (mfVar_)/(mf_*(1 - mf_) + SMALL) < 0.001)
        {
            fdelta_ = true;
        }
        else
        {
            a_ = max(mf_*gamma, 0.0);
            b_ = max((1.0 - mf_)*gamma, 0.0);
            if(a_ < 1.0)
            {
                delta0_ = true;
            }
            if(b_ < 1.0)
            {
                delta1_ = true;
            }
            _limitAB();
        }
 }
 void Foam::CMC::BetaFunction::_limitAB()
 {
    scalar fmax = 1.0/(1.0 + (b_ - 1.0)/(a_ - 1.0));
    if(a_ > 500.0)
    {
        a_ = 500.0;
        b_ = (a_ - 1.0 - fmax*(a_ - 2.0))/fmax;
    }
    else if(b_ > 500.0)
    {
        b_ = 500.0;
        a_ = (1.0 + fmax*(b_ - 2.0))/(1.0 - fmax);
    }
 }
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 Foam::CMC::BetaFunction::BetaFunction(const scalar mf, const scalar mfVar)
 :
    mf_(mf),
    mfVar_(mfVar),
    a_(0.0),
    b_(0.0),
    fdelta_(false),
    delta0_(false),
    delta1_(false)
 {
    _setParameters();
 }
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 Foam::CMC::BetaFunction::~BetaFunction()
 {}
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/BetaFunction/BetaFunction.H
+++ b/applications/solvers/combustion/SLFMFoam/BetaFunction/BetaFunction.H
@ -1,199 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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::BetaFunction
 Description
 SourceFiles
    BetaFunctionI.H
    BetaFunction.C
    BetaFunctionIO.C
 \*---------------------------------------------------------------------------*/
 #ifndef BetaFunction_H
 #define BetaFunction_H
 #include "scalar.H"
 #include "Switch.H"
 #include "scalarField.H"
 #include "tmp.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 namespace Foam
 {
 // Forward declaration of classes
 class Istream;
 class Ostream;
 // Forward declaration of friend functions and operators
 namespace CMC
 {
    class BetaFunction;
 }
 Istream& operator>>(Istream&, CMC::BetaFunction&);
 Ostream& operator<<(Ostream&, const CMC::BetaFunction&);
 /*---------------------------------------------------------------------------*\
                         Class BetaFunction Declaration
 \*---------------------------------------------------------------------------*/
 namespace CMC
 {
 class BetaFunction
 {
    // Private data
        //- Mean mixture fraction
        scalar mf_;
        //- Mixture fraction variance
        scalar mfVar_;
        //- Beta pdf parameter alpha
        scalar a_;
        //- Beta pdf parameter beta
        scalar b_;
        //- Flag beta pdf special case: forced-delta function
        Switch fdelta_;
        //- delta function at pure oxidizer stream
        Switch delta0_;
        //- delta function at pure fuel stream
        Switch delta1_;
    // Private Member Functions
        //- Calculate alpha and beta and set special case flags
        void _setParameters();
        //- Limit too large value of beta pdf parameters
        void _limitAB();
        //- Disallow default bitwise copy construct
        // BetaFunction(const BetaFunction&);
        //- Disallow default bitwise assignment
        // void operator=(const BetaFunction&);
 public:
    // Static data members
        //- Static data staticData
        // static const scalar staticData;
    // Constructors
        //- Construct from components
        BetaFunction(const scalar mf, const scalar mfVar);
        //- Construct from Istream
        BetaFunction(Istream&);
        //- Construct as copy
        // BetaFunction(const BetaFunction&);
    //- Destructor
    ~BetaFunction();
    // Member Functions
        tmp<scalarField> etaFunc(const UList<scalar>& eta) const
        {
            return pow(eta, a_-1.0)*pow(1.0-eta, b_-1.0);
        }
        scalar etaFunc(const scalar eta) const
        {
            return Foam::pow(eta, a_-1.0)*Foam::pow(1.0-eta, b_-1.0);
        }
        // Access
            Switch fdelta() const {return fdelta_;}
            Switch delta0() const {return delta0_;}
            Switch delta1() const {return delta1_;}
            scalar mf() const {return mf_;}
            scalar mfVar() const {return mfVar_;}
            scalar alpha() const {return a_;}
            scalar beta() const {return b_;}
        // Check
        // Edit
        // Write
    // Member Operators
        // void operator=(const BetaFunction&);
    // Friend Functions
    // Friend Operators
    // IOstream Operators
        friend Istream& ::Foam::operator>>(Istream&, BetaFunction&);
        friend Ostream& ::Foam::operator<<(Ostream&, const BetaFunction&);
 };
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 } // End namespace CMC
 } // End namespace Foam
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #include "BetaFunctionI.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #endif
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/BetaFunction/BetaFunctionI.H
+++ b/applications/solvers/combustion/SLFMFoam/BetaFunction/BetaFunctionI.H
@ -1,58 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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/>.
 \*---------------------------------------------------------------------------*/
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/BetaFunction/BetaFunctionIO.C
+++ b/applications/solvers/combustion/SLFMFoam/BetaFunction/BetaFunctionIO.C
@ -1,65 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "BetaFunction.H"
 #include "IOstreams.H"
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 Foam::CMC::BetaFunction::BetaFunction(Istream& is)
 {
    // Check state of Istream
    is.check("Foam::CMC::BetaFunction::BetaFunction(Foam::Istream&)");
 }
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 Foam::Istream& Foam::operator>>(Istream& is, CMC::BetaFunction&)
 {
    // Check state of Istream
    is.check
    (
        "Foam::Istream& Foam::operator>>(Foam::Istream&, Foam::CMC::BetaFunction&)"
    );
    return is;
 }
 Foam::Ostream& Foam::operator<<(Ostream& os, const CMC::BetaFunction&)
 {
    // Check state of Ostream
    os.check
    (
        "Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
        "const Foam::CMC::BetaFunction&)"
    );
    return os;
 }
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/BetaGrid/BetaGrid.C
+++ b/applications/solvers/combustion/SLFMFoam/BetaGrid/BetaGrid.C
@ -1,126 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "BetaGrid.H"
 #include "error.H"
 #include "dictionary.H"
 #include "SubList.H"
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 // const scalar Foam::CMC::BetaGrid::staticData();
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 Foam::CMC::BetaGrid::BetaGrid(const dictionary& dict)
 :
    etaCut_(dict.lookup("detailedEta")),
    nSpacings_(dict.lookup("detailedN")),
    etaSpace_(sum(nSpacings_)+3)
 {
    // Check for input list size
    if (etaCut_.size() != nSpacings_.size() + 1)
    {
        FatalErrorInFunction
            << "Number of grid interval points does not match number of grid sizes"
            << abort(FatalError);
    }
    // Check for evenness of number of spacings in each interval
    forAll (nSpacings_, i)
    {
        if ((nSpacings_[i] % 2) == 1)
        {
            FatalErrorInFunction
                << "Number of spacings in intervals should be even"
                << abort(FatalError);
        }
    }
    forAll(nSpacings_, i)
    {
        const label ns = nSpacings_[i];
        const labelList::subList prev(nSpacings_, i);
        const label baseI = sum(prev) + 1;
        const scalar low = etaCut_[i];
        const scalar upp = etaCut_[i+1];
        const scalar delta = (upp - low) / scalar(ns);
        for (label j = 0; j < ns; j++)
        {
            etaSpace_[baseI+j] = low + delta * j;
        }
    }
    etaSpace_[0] = 0.0;
    etaSpace_[etaSpace_.size()-2] = etaCut_.last();
    etaSpace_[etaSpace_.size()-1] = 1.0;
 }
 Foam::CMC::BetaGrid::BetaGrid(const BetaGrid&)
 {}
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 Foam::CMC::BetaGrid::~BetaGrid()
 {}
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
 void Foam::CMC::BetaGrid::operator=(const BetaGrid& rhs)
 {
    // Check for assignment to self
    if (this == &rhs)
    {
        FatalErrorInFunction
            << "Attempted assignment to self"
            << abort(FatalError);
    }
 }
 // * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/BetaGrid/BetaGrid.H
+++ b/applications/solvers/combustion/SLFMFoam/BetaGrid/BetaGrid.H
@ -1,161 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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::BetaGrid
 Description
 SourceFiles
    BetaGridI.H
    BetaGrid.C
    BetaGridIO.C
 \*---------------------------------------------------------------------------*/
 #ifndef BetaGrid_H
 #define BetaGrid_H
 #include "scalarField.H"
 #include "labelList.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 namespace Foam
 {
 // Forward declaration of classes
 class Istream;
 class Ostream;
 // Forward declaration of friend functions and operators
 namespace CMC
 {
    class BetaGrid;
 }
 class dictionary;
 Istream& operator>>(Istream&, CMC::BetaGrid&);
 Ostream& operator<<(Ostream&, const CMC::BetaGrid&);
 /*---------------------------------------------------------------------------*\
                         Class BetaGrid Declaration
 \*---------------------------------------------------------------------------*/
 namespace CMC
 {
 class BetaGrid
 {
    // Private data
        //- Piecewise uniform interval points
        scalarField etaCut_;
        //- Sizes of piecewise uniform intervals
        labelList nSpacings_;
        //- Mixture fraction eta grid points
        scalarField etaSpace_;
    // Private Member Functions
        //- Disallow default bitwise copy construct
        BetaGrid(const BetaGrid&);
        //- Disallow default bitwise assignment
        void operator=(const BetaGrid&);
 public:
    // Static data members
        //- Static data staticData
        // static const scalar staticData;
    // Constructors
        //- Construct from components
        BetaGrid(const dictionary& dict);
        //- Construct from Istream
        BetaGrid(Istream&);
        //- Construct as copy
        // BetaGrid(const BetaGrid&);
    //- Destructor
    ~BetaGrid();
    // Member Functions
        // Access
        const labelList& N() const {return nSpacings_;}
        const scalarField& etaCut() const {return etaCut_;}
        const scalarField& etaSpace() const {return etaSpace_;}
        // Check
        // Edit
        // Write
    // Member Operators
        // void operator=(const BetaGrid&);
    // Friend Functions
    // Friend Operators
    // IOstream Operators
        friend Istream& ::Foam::operator>>(Istream&, BetaGrid&);
        friend Ostream& ::Foam::operator<<(Ostream&, const BetaGrid&);
 };
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 } // End namespace CMC
 } // End namespace Foam
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #include "BetaGridI.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #endif
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/BetaGrid/BetaGridI.H
+++ b/applications/solvers/combustion/SLFMFoam/BetaGrid/BetaGridI.H
@ -1,58 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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/>.
 \*---------------------------------------------------------------------------*/
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/BetaGrid/BetaGridIO.C
+++ b/applications/solvers/combustion/SLFMFoam/BetaGrid/BetaGridIO.C
@ -1,65 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "BetaGrid.H"
 #include "IOstreams.H"
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 Foam::CMC::BetaGrid::BetaGrid(Istream& is)
 {
    // Check state of Istream
    is.check("Foam::CMC::BetaGrid::BetaGrid(Foam::Istream&)");
 }
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 Foam::Istream& Foam::operator>>(Istream& is, CMC::BetaGrid&)
 {
    // Check state of Istream
    is.check
    (
        "Foam::Istream& Foam::operator>>(Foam::Istream&, Foam::CMC::BetaGrid&)"
    );
    return is;
 }
 Foam::Ostream& Foam::operator<<(Ostream& os, const CMC::BetaGrid&)
 {
    // Check state of Ostream
    os.check
    (
        "Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
        "const Foam::CMC::BetaGrid&)"
    );
    return os;
 }
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/BetaIntegrator/BetaIntegrator.C
+++ b/applications/solvers/combustion/SLFMFoam/BetaIntegrator/BetaIntegrator.C
@ -1,180 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "BetaIntegrator.H"
 #include "error.H"
 #include "SubList.H"
 #include "SubField.H"
 #include "interpolateXY.H"
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 // const scalar Foam::CMC::BetaIntegrator::staticData();
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 //extended Simpson's rule (Numerical recipes, 2nd Ed. p.128)
 //for equally spaced and even N intervals (or odd N+1 points)
 Foam::scalar Foam::CMC::BetaIntegrator::simps
 (const scalar xl, const scalar xh, const label n, const UList<scalar>& fx)
    const
 {
    scalar evensum(0.0), oddsum(0.0), sum(0.0);
    scalar h = (xh - xl)/scalar(n);
    for(label i=0 ; i<fx.size() ; i++)
    {
        if(i%2 == 0)
        {
            evensum += fx[i];
        }
        else
        {
            oddsum += fx[i];
        }
    }
    sum = -1.0*fx.first() + 2.0*evensum + 4.0*oddsum - 1.0*fx.last();
    return sum*(h/3.0);
 }
 Foam::scalar Foam::CMC::BetaIntegrator::sumSimps(const UList<scalar>& f) const
 {
    scalar total = 0.0;
    const labelList &N(bg_.N());
    const scalarField &etaCut(bg_.etaCut());
    forAll(N, i)
    {
        const labelList::subList prev(N, i);
        const label baseI = sum(prev) + 1;
        const scalarField::subField subInterval(f, N[i]+1, baseI);
        total += simps(etaCut[i], etaCut[i+1], N[i], subInterval);
    }
    return total;
 }
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 Foam::CMC::BetaIntegrator::BetaIntegrator(const BetaFunction& bf, const BetaGrid &bg)
 :
    etaSpace_(bg.etaSpace()),
    pdfNum_(bg.etaSpace().size(), 0.0),
    pdfDen_(0.0),
    bf_(bf),
    bg_(bg)
 {
    const scalarField::subField domain(etaSpace_, etaSpace_.size()-2, 1);
    scalarField::subField range(pdfNum_, etaSpace_.size()-2, 1);
    if (bf_.fdelta())
    {
        pdfDen_ = 1.0;
    }
    else
    {
        range = bf_.etaFunc(domain);
        if (!bf_.delta0())
        {
            pdfNum_.first() = bf_.etaFunc(etaSpace_.first());
        }
        if (!bf_.delta1())
        {
            pdfNum_.last() = bf_.etaFunc(etaSpace_.last());
        }
        pdfDen_ = sumSimps(pdfNum_)
            + Foam::pow(etaSpace_[1], bf_.alpha())/(bf_.alpha() + SMALL)
            + Foam::pow(etaSpace_[1], bf_.beta())/(bf_.beta() + SMALL);
    }
 }
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 Foam::CMC::BetaIntegrator::~BetaIntegrator()
 {}
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 // beta-pdf weighted integration for given mf, mfVar and f
 Foam::scalar Foam::CMC::BetaIntegrator::betaIntegrate(const scalarField& f) const
 {
    scalar result = 0.0;
    if(bf_.fdelta())
    {
        result = interpolateXY(bf_.mf(), etaSpace_, f);
    }
    else
    {
        scalar num = sumSimps (f * pdfNum_)
            + f.first()*Foam::pow(etaSpace_[1], bf_.alpha())/(bf_.alpha() + SMALL)
            + f.last()*Foam::pow(etaSpace_[1], bf_.beta())/(bf_.beta() + SMALL);
        result = num/pdfDen_;
    }
    return result;
 }
 // * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
 void Foam::CMC::BetaIntegrator::operator=(const BetaIntegrator& rhs)
 {
    // Check for assignment to self
    if (this == &rhs)
    {
        FatalErrorInFunction
            << "Attempted assignment to self"
            << abort(FatalError);
    }
 }
 // * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/BetaIntegrator/BetaIntegrator.H
+++ b/applications/solvers/combustion/SLFMFoam/BetaIntegrator/BetaIntegrator.H
@ -1,168 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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::CMC::BetaIntegrator
 Description
 SourceFiles
    BetaIntegratorI.H
    BetaIntegrator.C
    BetaIntegratorIO.C
 \*---------------------------------------------------------------------------*/
 #ifndef BetaIntegrator_H
 #define BetaIntegrator_H
 #include "scalarField.H"
 #include "BetaFunction.H"
 #include "BetaGrid.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 namespace Foam
 {
 // Forward declaration of classes
 class Istream;
 class Ostream;
 // Forward declaration of friend functions and operators
 namespace CMC
 {
    class BetaIntegrator;
 }
 Istream& operator>>(Istream&, CMC::BetaIntegrator&);
 Ostream& operator<<(Ostream&, const CMC::BetaIntegrator&);
 /*---------------------------------------------------------------------------*\
                         Class BetaIntegrator Declaration
 \*---------------------------------------------------------------------------*/
 namespace CMC
 {
 class BetaIntegrator
 {
    // Private data
        //- Quadrature points
        const scalarField &etaSpace_;
        //- Beta pdf numerators
        scalarField pdfNum_;
        //- Beta pdf denominator
        scalar pdfDen_;
        //- Description of data_
        BetaFunction bf_;
        //- Description of data_
        const BetaGrid &bg_;
    // Private Member Functions
        //- Calculate integral using Simpson's rule
        scalar simps(const scalar xl, const scalar xh, const label n, const UList<scalar>& fx) const;
        //- Calculate integral using Simpson's rule
        scalar sumSimps(const UList<scalar>& f) const;
        //- Disallow default bitwise copy construct
        BetaIntegrator(const BetaIntegrator&);
        //- Disallow default bitwise assignment
        void operator=(const BetaIntegrator&);
 public:
    // Static data members
        //- Static data staticData
        // static const scalar staticData;
    // Constructors
        //- Construct from components
        BetaIntegrator(const BetaFunction& bf, const BetaGrid &bg);
        //- Construct as copy
        // BetaIntegrator(const BetaIntegrator&);
    //- Destructor
    ~BetaIntegrator();
    // Member Functions
        scalar betaIntegrate(const scalarField& f) const;
        // Access
        scalar pdfDen () const {return pdfDen_;}
        // Check
        // Edit
        // Write
    // Member Operators
        // void operator=(const BetaIntegrator&);
    // Friend Functions
    // Friend Operators
    // IOstream Operators
        friend Istream& ::Foam::operator>>(Istream&, BetaIntegrator&);
        friend Ostream& ::Foam::operator<<(Ostream&, const BetaIntegrator&);
 };
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 } // End namespace CMC
 } // End namespace Foam
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #include "BetaIntegratorI.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #endif
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/BetaIntegrator/BetaIntegratorI.H
+++ b/applications/solvers/combustion/SLFMFoam/BetaIntegrator/BetaIntegratorI.H
@ -1,58 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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/>.
 \*---------------------------------------------------------------------------*/
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/BetaIntegrator/BetaIntegratorIO.C
+++ b/applications/solvers/combustion/SLFMFoam/BetaIntegrator/BetaIntegratorIO.C
@ -1,56 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "BetaIntegrator.H"
 #include "IOstreams.H"
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 Foam::Istream& Foam::operator>>(Istream& is, CMC::BetaIntegrator&)
 {
    // Check state of Istream
    is.check
    (
        "Foam::Istream& Foam::operator>>(Foam::Istream&, Foam::CMC::BetaIntegrator&)"
    );
    return is;
 }
 Foam::Ostream& Foam::operator<<(Ostream& os, const CMC::BetaIntegrator&)
 {
    // Check state of Ostream
    os.check
    (
        "Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
        "const Foam::CMC::BetaIntegrator&)"
    );
    return os;
 }
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/FlameStructure/FlameStructure.C
+++ b/applications/solvers/combustion/SLFMFoam/FlameStructure/FlameStructure.C
@ -1,206 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "FlameStructure.H"
 #include "fileName.H"
 #include "IFstream.H"
 #include "hashedWordList.H"
 #include "interpolateXY.H"
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 Foam::SLFM::FlameStructure::FlameStructure(const string& NstDir, scalarList NstList, const scalarField &etaGrid, const hashedWordList &spTable)
 :
    NstDir_(NstDir),
    NstList_(NstList),
    eta_(etaGrid),
    species_(spTable),
    Y_(spTable.size(), scalarFieldArray1d(NstList_.size(), scalarField(eta_.size(), 0.0))),
    W_(spTable.size(), scalarFieldArray1d(NstList_.size(), scalarField(eta_.size(), 0.0))),
    T_(NstList_.size(), scalarField(eta_.size(), 0.0)),
    Q_(NstList_.size(), scalarField(eta_.size(), 0.0)),
    h_(NstList_.size(), scalarField(eta_.size(), 0.0)),
    rho_(NstList_.size(), scalarField(eta_.size(), 0.0)),
    Rgas_(NstList_.size(), scalarField(eta_.size(), 0.0))
 {
    for (label n = 0; n < NstList_.size(); n++)
    {
        scalar N = NstList_[n];
        fileName fname = "sdr"+Foam::name(N)+".inp";
        IFstream fin(NstDir_/fname);
        string gbg;
        label etamax_SLFM, NoSpecies;
        //INPUT FILE...
        fin.getLine(gbg);
        //No.GRID...
        fin.getLine(gbg);
        //Number of eta-point and species
        fin >> etamax_SLFM >> NoSpecies;
        //blank line
        fin.getLine(gbg);
        //PRESSURE...
        fin.getLine(gbg);
        //1.00...
        fin.getLine(gbg);
        //MIXTURE...
        fin.getLine(gbg);
        //read eta space (from SLFM library)
        scalarField etaValue_SLFM(etamax_SLFM, 0.0);
        for(label j = 0 ; j < etamax_SLFM ; j++)
        {
            fin>>etaValue_SLFM[j];
        }
        //blank line
        fin.getLine(gbg);
        //INITIAL...
        fin.getLine(gbg);
        //species loop
        scalarField temp(etamax_SLFM, 0.0);
        for(label i = 0 ; i < NoSpecies ; i++)
        {
            word spName;
            //H2... (species name)
            fin >> spName;
            const label spI (species_[spName]);
            // Info << spI << spName << endl;
            //read species mass fraction (from SLFM library)
            for(label j=0 ; j<etamax_SLFM ; j++)
            {
                scalar Yj = 0.0;
                fin >> Yj;
                temp[j] = max(0.0, Yj);
            }
            Y_[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, temp);
            //read reaction rate (from SLFM library)
            for(label j = 0 ; j < etamax_SLFM ; j++)
            {
                fin>>temp[j];
            }
            W_[spI][n] = interpolateXY(etaGrid, etaValue_SLFM, temp);
            //blank line
            fin.getLine(gbg);
        }
        //read temperature
        {
            //H2... (species name)
            fin.getLine(gbg);
            //read temperature
            for(label j = 0 ; j < etamax_SLFM ; j++)
            {
                scalar Yj = 0.0;
                fin >> Yj;
                temp[j] = max(0.0, Yj);
            }
            T_[n] = interpolateXY(etaGrid, etaValue_SLFM, temp);
            //read heat source
            for(label j = 0 ; j < etamax_SLFM ; j++)
            {
                fin>>temp[j];
            }
            Q_[n] = interpolateXY(etaGrid, etaValue_SLFM, temp);
            //blank line
            fin.getLine(gbg);
        }
        //DENSITY...
        fin.getLine(gbg);
        //read density (from SLFM library)
        for(label j = 0 ; j < etamax_SLFM ; j++)
        {
            fin>>temp[j];
        }
        //Info<<rho_temp<<endl;
        rho_[n] = interpolateXY(etaGrid, etaValue_SLFM, temp);
        //blank line
        fin.getLine(gbg);
        //ENTHALPY...
        fin.getLine(gbg);
        //read enthalpy (from SLFM library)
        for(label j = 0 ; j < etamax_SLFM ; j++)
        {
            fin>>temp[j];
        }
        h_[n] = interpolateXY(etaGrid, etaValue_SLFM, temp);
    }
 }
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 Foam::SLFM::FlameStructure::~FlameStructure()
 {}
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/FlameStructure/FlameStructure.H
+++ b/applications/solvers/combustion/SLFMFoam/FlameStructure/FlameStructure.H
@ -1,215 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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::SLFM::FlameStructure
 Description
 SourceFiles
    FlameStructureI.H
    FlameStructure.C
    FlameStructureIO.C
 \*---------------------------------------------------------------------------*/
 #ifndef FlameStructure_H
 #define FlameStructure_H
 #include "error.H"
 #include "scalarField.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 namespace Foam
 {
 // Forward declaration of classes
 class Istream;
 class Ostream;
 class hashedWordList;
 namespace SLFM
 {
    class FlameStructure;
 }
 // Forward declaration of friend functions and operators
 Istream& operator>>(Istream&, SLFM::FlameStructure&);
 Ostream& operator<<(Ostream&, const SLFM::FlameStructure&);
 namespace SLFM
 {
 /*---------------------------------------------------------------------------*\
                         Class FlameStructure Declaration
 \*---------------------------------------------------------------------------*/
 class FlameStructure
 {
    // Private data
        typedef List<scalarField> scalarFieldArray1d;
        typedef List<scalarFieldArray1d> scalarFieldArray2d;
        typedef List<scalarFieldArray2d> scalarFieldArray3d;
        //- Species mass fraction
        word NstDir_;
        //- Species mass fraction
        const scalarList NstList_;
        //- Species mass fraction
        const scalarField eta_;
        //- Species mass fraction
        const hashedWordList &species_;
        //- Species mass fraction
        scalarFieldArray2d Y_;
        //- Species production rate
        scalarFieldArray2d W_;
        //- Temperature
        scalarFieldArray1d T_;
        //- Heat source
        scalarFieldArray1d Q_;
        //- Enthalpy
        scalarFieldArray1d h_;
        //- Density
        scalarFieldArray1d rho_;
        //- Specific gas constant for mean W
        scalarFieldArray1d Rgas_;
    // Private Member Functions
        //- Disallow Construct null
        FlameStructure();
        //- Disallow default bitwise copy construct
        // FlameStructure(const FlameStructure&);
        //- Disallow default bitwise assignment
        // void operator=(const FlameStructure&);
 public:
    // Static data members
        //- Static data staticData
        // static const dataType staticData;
    // Constructors
        //- Construct from components
        FlameStructure(const string& NstDir, scalarList NstList, const scalarField &etaGrid, const hashedWordList& spTable);
        //- Construct from Istream
        // FlameStructure(Istream&);
        //- Construct as copy
        // FlameStructure(const FlameStructure&);
    // Selectors
        //- Select null constructed
        // static autoPtr<FlameStructure> New();
    //- Destructor
    ~FlameStructure();
    // Member Functions
        // Access
            //- Select null constructed
            scalarFieldArray2d& Y() {return Y_;}
            //- Select null constructed
            scalarFieldArray2d& W() {return W_;}
            //- Select null constructed
            scalarFieldArray1d& T() {return T_;}
            //- Select null constructed
            scalarFieldArray1d& Q() {return Q_;}
            //- Select null constructed
            scalarFieldArray1d& rho() {return rho_;}
            //- Select null constructed
            scalarFieldArray1d& h() {return h_;}
            //- Select null constructed
            scalarFieldArray1d& Rgas() {return Rgas_;}
        // Check
        // Edit
        // Write
    // Member Operators
        // void operator=(const FlameStructure&);
    // Friend Functions
    // Friend Operators
    // IOstream Operators
        friend Istream& operator>>(Istream&, FlameStructure&);
        friend Ostream& operator<<(Ostream&, const FlameStructure&);
 };
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 } // End namespace SLFM
 } // End namespace Foam
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #include "FlameStructureI.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #endif
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/FlameStructure/FlameStructureI.H
+++ b/applications/solvers/combustion/SLFMFoam/FlameStructure/FlameStructureI.H
@ -1,58 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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/>.
 \*---------------------------------------------------------------------------*/
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/FlameStructure/FlameStructureIO.C
+++ b/applications/solvers/combustion/SLFMFoam/FlameStructure/FlameStructureIO.C
@ -1,72 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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 "FlameStructure.H"
 #include "IOstreams.H"
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
    /*
 Foam::SLFM::FlameStructure::FlameStructure(Istream& is)
 :
    base1(is),
    base2(is),
    member1(is),
    member2(is)
 {
    // Check state of Istream
    is.check("Foam::SLFM::FlameStructure::FlameStructure(Foam::Istream&)");
 }
    */
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 Foam::Istream& Foam::operator>>(Istream& is, SLFM::FlameStructure&)
 {
    // Check state of Istream
    is.check
    (
        "Foam::Istream& Foam::operator>>(Foam::Istream&, Foam::SLFM::FlameStructure&)"
    );
    return is;
 }
 Foam::Ostream& Foam::operator<<(Ostream& os, const SLFM::FlameStructure&)
 {
    // Check state of Ostream
    os.check
    (
        "Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
        "const Foam::SLFM::FlameStructure&)"
    );
    return os;
 }
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/Make/files
+++ b/applications/solvers/combustion/SLFMFoam/Make/files
@ -1,17 +0,0 @@
 FlameStructure/FlameStructure.C
 FlameStructure/FlameStructureIO.C
 ./AMC/AMC.C
 ./BetaFunction/BetaFunction.C
 ./BetaFunction/BetaFunctionIO.C
 ./BetaGrid/BetaGrid.C
 ./BetaGrid/BetaGridIO.C
 ./BetaIntegrator/BetaIntegrator.C
 ./BetaIntegrator/BetaIntegratorIO.C
 SLFMFoam.C
 EXE = $(FOAM_APPBIN)/SLFMFoam
--- a/applications/solvers/combustion/SLFMFoam/Make/options
+++ b/applications/solvers/combustion/SLFMFoam/Make/options
@ -1,34 +0,0 @@
 EXE_INC = \
    -I. \
    -I./AMC \
    -I./BetaFunction \
    -I./BetaGrid \
    -I./BetaIntegrator \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
    -I$(LIB_SRC)/ODE/lnInclude \
    -I$(LIB_SRC)/combustionModels/lnInclude \
    -I$(LIB_SRC)/sampling/lnInclude
 EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lreactionThermophysicalModels \
    -lchemistryModel \
    -lODE \
    -lcombustionModels \
    -lfvOptions \
    -lsampling
--- a/applications/solvers/combustion/SLFMFoam/Mixturefraction.H
+++ b/applications/solvers/combustion/SLFMFoam/Mixturefraction.H
@ -1,10 +0,0 @@
 // mixture fraction equation
 fvScalarMatrix mfEqn
 (
    fvm::div(phi, mf)
  - fvm::laplacian((1/Sc)*turbulence->mut(), mf) //let 1/Sc = 1.47
 );
 mfEqn.relax();
 mfEqn.solve();
--- a/applications/solvers/combustion/SLFMFoam/MixturefractionVar.H
+++ b/applications/solvers/combustion/SLFMFoam/MixturefractionVar.H
@ -1,17 +0,0 @@
 //mixture fraction variance equation
 //cf) mfVar should not be ZERO
 SDR = turbulence->epsilon() * mfVar / turbulence->k();
 volVectorField Gradmf = fvc::grad(mf);
 fvScalarMatrix mfVarEqn
 (
    fvm::div(phi, mfVar)
  - fvm::laplacian((1/Sc) * turbulence->mut(), mfVar)
 ==
    2 * ((1/Sc) * turbulence->mut()) * (Gradmf & Gradmf)
  - 2 * rho * SDR
 );
 mfVarEqn.relax();
 mfVarEqn.solve();
--- a/applications/solvers/combustion/SLFMFoam/SLFMFoam.C
+++ b/applications/solvers/combustion/SLFMFoam/SLFMFoam.C
@ -1,118 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013-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/>.
 Application
    SLFMFoam
 Description
    Steady Laminar Flamelet Model(SLFM) solver for turbulent combustion. SLFM
    assumes a turbulent flame composed of thin stretched laminar flamelets.
    Each flamelet library is generated by external program and imported by
    the solver.
 References
    A.Y. Klimenko, R.W. Bilger, Progress in Energy and Combustion Science 25 (1999) 595-687
    N. Peters, Turbulent Combustion, Cambridge University Press (2000)
 Contact
    POSTECH combustion lab.
    huh@postech.ac.kr
 \*---------------------------------------------------------------------------*/
 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
 #include "rhoCombustionModel.H"
 #include "fvOptions.H"
 #include "simpleControl.H"
 #include "FlameStructure/FlameStructure.H"
 #include "interpolateXY.H"
 #include "AMC.H"
 #include "BetaFunction.H"
 #include "BetaGrid.H"
 #include "BetaIntegrator.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 int main(int argc, char *argv[])
 {
    #include "postProcess.H"
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createControl.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
    #include "createFvOptions.H"
    #include "initContinuityErrs.H"
    turbulence->validate();
    //SLFM
    Info<<"Read SLFM library"<<endl;
    #include "readSLFMlib.H" //read SLFM library
    Info<<"Construct Beta-PDF"<<endl;
    #include "preIntegration.H" //AMC C1 coefficient pdf integration
    Info<< "\nStarting time loop\n" << endl;
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    Info<< "\nStarting time loop\n" << endl;
    while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;
        // --- Pressure-velocity SIMPLE corrector loop
        {
            #include "UEqn.H"
            #include "Mixturefraction.H" //mean mixture fraction
            #include "MixturefractionVar.H" //mean mixture fraction variance
            #include "setSDR.H" //calculate scalar dissipation rate
            #include "updateT_RHO.H" //update temperature and density
            #include "pEqn.H"
        }
        turbulence->correct();
        if (runTime.write())
        {
            #include "updateYi.H" //update species
        }
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }
    Info<< "End\n" << endl;
    return 0;
 }
 // ************************************************************************* //
--- a/applications/solvers/combustion/SLFMFoam/UEqn.H
+++ b/applications/solvers/combustion/SLFMFoam/UEqn.H
@ -1,20 +0,0 @@
    MRF.correctBoundaryVelocity(U);
    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(U)
     ==
        rho()*g
      + fvOptions(rho, U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();
    UEqn.relax();
    fvOptions.constrain(UEqn);
    solve(UEqn == -fvc::grad(p));
    fvOptions.correct(U);
--- a/applications/solvers/combustion/SLFMFoam/createFieldRefs.H
+++ b/applications/solvers/combustion/SLFMFoam/createFieldRefs.H
@ -1,2 +0,0 @@
 const volScalarField& psi = thermo.psi();
 //const volScalarField& T = thermo.T();
--- a/applications/solvers/combustion/SLFMFoam/createFields.H
+++ b/applications/solvers/combustion/SLFMFoam/createFields.H
@ -1,374 +0,0 @@
 #include "readGravitationalAcceleration.H"
 Info<< "Creating combustion model\n" << endl;
 autoPtr<combustionModels::rhoCombustionModel> combustion
 (
    combustionModels::rhoCombustionModel::New(mesh)
 );
 rhoReactionThermo& thermo = combustion->thermo();
 thermo.validate(args.executable(), "h", "e");
 basicSpecieMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
 const word inertSpecie(thermo.lookup("inertSpecie"));
 if (!composition.contains(inertSpecie))
 {
    FatalErrorInFunction
        << "Specified inert specie '" << inertSpecie << "' not found in "
        << "species list. Available species:" << composition.species()
        << exit(FatalError);
 }
 volScalarField& p = thermo.p();
 volScalarField& T = thermo.T();
 wordList rhoBoundaryTypes
 (
 	T.boundaryField().size(),
 	zeroGradientFvPatchScalarField::typeName
 );
 forAll(composition.Y(inertSpecie).boundaryField(), patchi)
 {
 	rhoBoundaryTypes[patchi]
        = composition.Y(inertSpecie).boundaryField()[patchi].type();
 }
 volScalarField rho
 (
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    thermo.rho(),
    rhoBoundaryTypes
 );
 Info<< "\nReading field U\n" << endl;
 volVectorField U
 (
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
 );
 #include "compressibleCreatePhi.H"
 mesh.setFluxRequired(p.name());
 dimensionedScalar rhoMax
 (
    dimensionedScalar::lookupOrDefault
    (
        "rhoMax",
        simple.dict(),
        dimDensity,
        GREAT
    )
 );
 dimensionedScalar rhoMin
 (
    dimensionedScalar::lookupOrDefault
    (
        "rhoMin",
        simple.dict(),
        dimDensity,
        0
    )
 );
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
    compressible::turbulenceModel::New
    (
        rho,
        U,
        phi,
        thermo
    )
 );
 // Set the turbulence into the combustion model
 combustion->setTurbulence(turbulence());
 Info<< "Creating multi-variate interpolation scheme\n" << endl;
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 forAll(Y, i)
 {
    fields.add(Y[i]);
 }
 fields.add(thermo.he());
 // SLFM-related fields
 Info<<"Creating field mf, mfVar, SDR\n"<<endl;
 volScalarField mf //mixture fraction
 (
    IOobject
    (
        "mf",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
 );
 volScalarField mfVar //mixture fraction variance
 (
    IOobject
    (
        "mfVar",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
 );
 volScalarField SDR //scalar dissipation rate
 (
    IOobject
    (
        "SDR",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("SDR", dimEnergy/dimTime/dimEnergy, 0.0)
 );
 volScalarField jlc //eta-index, lower
 (
    IOobject
    (
        "jlc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("jlc", dimless, 0.0)
 );
 volScalarField jhc //eta-index, higher
 (
    IOobject
    (
        "jhc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("jhc", dimless, 0.0)
 );
 volScalarField jfc //eta-factor
 (
    IOobject
    (
        "jfc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("jfc", dimless, 0.0)
 );
 volScalarField vlc //var-index, lower
 (
    IOobject
    (
        "vlc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("vlc", dimless, 0.0)
 );
 volScalarField vhc //var-index, higher
 (
    IOobject
    (
        "vhc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("vhc", dimless, 0.0)
 );
 volScalarField vfc //var-factor
 (
    IOobject
    (
        "vfc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("vfc", dimless, 0.0)
 );
 volScalarField nlc //sdr-index, lower
 (
    IOobject
    (
        "nlc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("nlc", dimless, 0.0)
 );
 volScalarField nhc //sdr-index, higher
 (
    IOobject
    (
        "nhc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("nhc", dimless, 0.0)
 );
 volScalarField nfc //sdr-factor
 (
    IOobject
    (
        "nfc",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("nfc", dimless, 0.0)
 );
 //SLFM property read
 IOdictionary SLFMdict //SLFM properties dictionary
 (
    IOobject
    (
        "SLFMdict",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ,
        IOobject::NO_WRITE
    )
 );
 const wordList postSpecieNames(SLFMdict.lookup("postSpecieNames"));
 labelList postSpecieIndices(postSpecieNames.size(),-1);
 forAll(postSpecieNames, i)
 {
    const word &spName (postSpecieNames[i]);
    const label spI (composition.species()[spName]);
    postSpecieIndices[i] = spI;
 }
 scalar stoiMF(readScalar(SLFMdict.lookup("stoiMF")));
 scalar stretchFac(readScalar(SLFMdict.lookup("stretchFactor")));
 scalar lowerN(readScalar(SLFMdict.lookup("lowerN")));
 scalar upperN(readScalar(SLFMdict.lookup("upperN")));
 scalar Sc(readScalar(SLFMdict.lookup("Sc")));
 const label etamax(lowerN+upperN+1);
 word outletName(SLFMdict.lookup("outletName"));
 labelList nIntervals (SLFMdict.lookup("detailedN"));
 CMC::BetaGrid bg(SLFMdict);
 const scalarField &etaGrid(bg.etaSpace());
 const scalarField amcGridValue(CMC::AMC(etaGrid));
 Info << etaGrid << endl;
 scalarField etaValue(etamax,0);
 scalarField etaIndex(etamax,0);
 scalarField deta_h(etamax,0);
 scalarField deta_l(etamax,0);
 scalarField deta_m(etamax,0);
 #include "setEtaSpace.H"
 Info<<"\nNumber of eta-point = "<<etamax<<endl;
 scalar NVar(readScalar(SLFMdict.lookup("NVar")));
 scalarField varValue(NVar+1,0);
 scalarField varIndex(NVar+1,0);
 #include "setVarSpace.H"
 PtrList<volScalarField> Neta(etamax);
 for(label j=0 ; j<etamax ; j++)
 {
    word Netai = "Neta_"+Foam::name(j);
    Neta.set
    (
        j,
        new volScalarField
        (
            IOobject
            (
                Netai,
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            dimensionedScalar(Netai, dimless, 0.0)
        )
    );
 }
 #include "createMRF.H"
--- a/applications/solvers/combustion/SLFMFoam/pEqn.H
+++ b/applications/solvers/combustion/SLFMFoam/pEqn.H
@ -1,55 +0,0 @@
 {
    // Thermodynamic density needs to be updated by psi*d(p) after the
    // pressure solution - done in 2 parts. Part 1:
    thermo.rho() -= psi*p;
    volScalarField rAU(1.0/UEqn.A());
    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
    tUEqn.clear();
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        fvc::interpolate(rho)*fvc::flux(HbyA)
    );
    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
    while (simple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvc::div(phiHbyA)
          - fvm::laplacian(rhorAUf, p)
         ==
            fvOptions(psi, p, rho.name())
        );
        pEqn.solve();
        if (simple.finalNonOrthogonalIter())
        {
            phi = phiHbyA + pEqn.flux();
        }
    }
    p.relax();
    // Second part of thermodynamic density update
    thermo.rho() += psi*p;
    #include "compressibleContinuityErrs.H"
    U = HbyA - rAU*fvc::grad(p);
    U.correctBoundaryConditions();
    fvOptions.correct(U);
    rho = max(rho, rhoMin);
    rho = min(rho, rhoMax);
    rho.relax();
    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
 }
--- a/applications/solvers/combustion/SLFMFoam/preIntegration.H
+++ b/applications/solvers/combustion/SLFMFoam/preIntegration.H
@ -1,78 +0,0 @@
 scalarFieldArray1d C1table //AMC C1 coefficient table
 (
    etamax,
    scalarField(NVar+1, 0)
 );
 scalarFieldArray2d Ttable
 (
    etamax,
    scalarFieldArray1d
    (
        NVar+1,
        scalarField(NstList.size(), 0)
    )
 );
 scalarFieldArray2d Rtable
 (
    etamax,
    scalarFieldArray1d
    (
        NVar+1,
        scalarField(NstList.size(), 0)
    )
 );
 scalarFieldArray3d Ytable
 (
    Y.size(),
    scalarFieldArray2d
    (
        etamax,
        scalarFieldArray1d
        (
            NVar+1,
            scalarField(NstList.size(), 0)
        )
    )
 );
 Info<<"Construct C1 coefficient table"<<endl;
 for(label j = 1; j < etamax - 1; j++)
 {
    for(label v = 1; v < varValue.size() - 1; v++)
    {
        const scalar mf = etaValue[j];
        const scalar mfVar = varValue[v]*mf*(1.0-mf);
        CMC::BetaFunction bf(mf, mfVar);
        CMC::BetaIntegrator bi(bf, bg);
        C1table[j][v] = 1.0/bi.betaIntegrate(amcGridValue);
    }
 }
 Info<<"Construct Temp and Rgas table"<<endl;
 for(label j=0 ; j<etamax ; j++)
 {
    for(label v=0 ; v<=NVar ; v++)
    {
        const scalar mf = etaValue[j];
        const scalar mfVar = varValue[v]*mf*(1.0-mf);
        CMC::BetaFunction bf(mf, mfVar);
        CMC::BetaIntegrator bi(bf, bg);
        for(label n=0 ; n<NstList.size() ; n++)
        {
            Ttable[j][v][n] = bi.betaIntegrate(T_SLFM[n]);
            Rtable[j][v][n] = bi.betaIntegrate(Rgas_SLFM[n]);
            forAll(Y, yi)
            {
                Ytable[yi][j][v][n] = bi.betaIntegrate(Y_SLFM[yi][n]);
            }
        }
    }
    Info<<scalar(j+1)/etamax*100<<"%"<<endl;
 }
--- a/applications/solvers/combustion/SLFMFoam/readSLFMlib.H
+++ b/applications/solvers/combustion/SLFMFoam/readSLFMlib.H
@ -1,33 +0,0 @@
 Switch Equilibrium = SLFMdict.lookup("Equilibrium");
 word NstFolder = SLFMdict.lookup("NstFolder");
 scalarField NstList = SLFMdict.lookup("NstList");
 scalarField NstIndex(NstList.size(),0);
 forAll(NstList, n)
 {
    NstIndex[n] = n;
 }
 typedef List<scalarField> scalarFieldArray1d;
 typedef List<scalarFieldArray1d> scalarFieldArray2d;
 typedef List<scalarFieldArray2d> scalarFieldArray3d;
 SLFM::FlameStructure slfmLibrary(runTime.constant()/NstFolder, NstList, etaGrid, composition.species());
 scalarFieldArray2d& Y_SLFM (slfmLibrary.Y());
 scalarFieldArray2d& W_SLFM (slfmLibrary.W());
 scalarFieldArray1d& T_SLFM (slfmLibrary.T());
 scalarFieldArray1d& Q_SLFM (slfmLibrary.Q());
 scalarFieldArray1d& rho_SLFM (slfmLibrary.rho());
 scalarFieldArray1d& h_SLFM (slfmLibrary.h());
 scalarFieldArray1d& Rgas_SLFM (slfmLibrary.Rgas());
 for(label n=0 ; n<NstList.size() ; n++)
 {
    for(label i=0 ; i<Y.size() ; i++)
    {
        Rgas_SLFM[n] += Y_SLFM[i][n]/composition.W(i); // = 1/meanMW
    }
    Rgas_SLFM[n] = 8314.46*Rgas_SLFM[n]; //universal gas constant = 8314.46 [J/K.Kmole]
 }
--- a/applications/solvers/combustion/SLFMFoam/setEtaSpace.H
+++ b/applications/solvers/combustion/SLFMFoam/setEtaSpace.H
@ -1,51 +0,0 @@
 //set eta-space
 scalar upperSum(0), upperDelta(0);
 scalar lowerSum(0), lowerDelta(0);
 //summation of geometric progression (eta > stoi)
 upperSum = Foam::pow(stretchFac,upperN)-1.0;
 upperSum = upperSum/(stretchFac-1.0);
 //summation of geometric progression (eta < stoi)
 lowerSum = Foam::pow(stretchFac,lowerN)-1.0;
 lowerSum = lowerSum/(stretchFac-1.0);
 //set reference delta eta
 upperDelta = (1.0-stoiMF)/upperSum;
 lowerDelta = stoiMF/lowerSum;
 //set boundary and stoichiometric value of eta
 etaValue[0] = 0.0;
 etaValue[lowerN] = stoiMF;
 etaValue.last() = 1.0;
 //set internal value of eta (eta > stoi)
 for(label j=0 ; j<upperN-1 ; j++)
 {
    etaValue[lowerN+j+1] = etaValue[lowerN+j] + upperDelta*pow(stretchFac,j);
 }
 //set internal value of eta (eta < stoi)
 for(label j=0 ; j<lowerN-1 ; j++)
 {
    etaValue[lowerN-j-1] = etaValue[lowerN-j] - lowerDelta*pow(stretchFac,j);
 }
 Info<<etaValue<<endl;
 //calculate delta eta
 for(label j=1 ; j<etamax-1 ; j++)
 {
    deta_h[j] = etaValue[j+1] - etaValue[j];
    deta_l[j] = etaValue[j] - etaValue[j-1];
    deta_m[j] = 0.5*(etaValue[j+1] - etaValue[j-1]);
 }
 //Info<<deta_h<<endl;
 //Info<<deta_l<<endl;
 //Info<<deta_m<<endl;
 for(label j=0 ; j<etamax ; j++)
 {
    etaIndex[j] = j;
 }
--- a/applications/solvers/combustion/SLFMFoam/setSDR.H
+++ b/applications/solvers/combustion/SLFMFoam/setSDR.H
@ -1,79 +0,0 @@
 scalarField AMCcoeff(CMC::AMC(etaValue));
 //cell loop
 forAll(mf, cellI)
 {
    scalar jl(0), jh(0), vl(0), vh(0);
    scalar jfac(0), vfac(0);
    scalar C1coeff(0);
    //find eta-index, factor
    if(mf[cellI] < etaValue[1])
    {
        jl = 0;
        jh = 1;
    }
    else if(mf[cellI] > etaValue[etamax-2])
    {
        jl = etamax-2;
        jh = etamax-1;
    }
    else
    {
        jl = label( interpolateXY(mf[cellI], etaValue, etaIndex) );
        jh = jl+1;
    }
    jfac = max(0.0, (mf[cellI]-etaValue[jl])/(etaValue[jh]-etaValue[jl]));
    //find var-index, factor
    scalar scaledVar = min(0.99999, mfVar[cellI]/(mf[cellI]*(1.0-mf[cellI])+SMALL));
    if(scaledVar < varValue[1])
    {
        vl = 0;
        vh = 1;
    }
    else if(scaledVar > varValue[NVar-1])
    {
        vl = NVar-1;
        vh = NVar;
    }
    else
    {
        vl = label( interpolateXY(scaledVar, varValue, varIndex) );
        vh = vl+1;
    }
    vfac = max(0.0, (scaledVar-varValue[vl])/(varValue[vh]-varValue[vl]));
    //Bi-linear interpolation on j and v
    //Numerical recipes, 2nd Ed. p.117
    C1coeff
    = C1table[jl][vl]*(1-jfac)*(1-vfac)
    + C1table[jh][vl]*(jfac)*(1-vfac)
    + C1table[jh][vh]*(jfac)*(vfac)
    + C1table[jl][vh]*(1-jfac)*(vfac);
    jlc[cellI] = jl;
    jhc[cellI] = jh;
    jfc[cellI] = jfac;
    vlc[cellI] = vl;
    vhc[cellI] = vh;
    vfc[cellI] = vfac;
    scalar Coeff = AMCcoeff[lowerN];
    if (Equilibrium == true)
    {
        Neta[lowerN][cellI] = 0; //Equilibrium
    }
    else
    {
        Neta[lowerN][cellI] = SDR[cellI]*Coeff*C1coeff; //CSDR at stoichiometic m.f.
    }
 }
--- a/applications/solvers/combustion/SLFMFoam/setVarSpace.H
+++ b/applications/solvers/combustion/SLFMFoam/setVarSpace.H
@ -1,23 +0,0 @@
 //set normalized variance space
 {
    scalar sum(0), del(0), fac(1.075);
    sum = Foam::pow(fac,NVar) - 1.0;
    sum = sum/(fac-1.0);
    del = 1.0/sum;
    varValue[0] = 0.0;
    varValue[NVar] = 1.0;
    for(label v=0 ; v<NVar-1 ; v++)
    {
        varValue[v+1] = varValue[v] + del*Foam::pow(fac,v);
    }
    for(label v=0 ; v<=NVar ; v++)
    {
        varIndex[v] = v;
    }
    Info<<varValue<<endl;
 }
--- a/applications/solvers/combustion/SLFMFoam/updateT_RHO.H
+++ b/applications/solvers/combustion/SLFMFoam/updateT_RHO.H
@ -1,60 +0,0 @@
 //force the saving of the old-time values
 //important!!
 rho.oldTime();
 forAll(rho, cellI)
 {
    scalar& Nst = Neta[lowerN][cellI];
    //find Nst index
    if(Nst < NstList.first())
    {
        nlc[cellI] = 0;
        nfc[cellI] = 0;
    }
    else if(Nst > NstList[NstList.size()-2])
    {
        nhc[cellI] = NstList.size()-2; //just before extinction
        nfc[cellI] = 1;
    }
    else
    {
        nlc[cellI] = label( interpolateXY(Nst, NstList, NstIndex) );
        nhc[cellI] = nlc[cellI]+1;
        nfc[cellI] = (Nst-NstList[nlc[cellI]])
                     /(NstList[nhc[cellI]]-NstList[nlc[cellI]]);
    }
    //update Favre mean temp. and rho(=1/RT)
    scalar jl = jlc[cellI], jh = jhc[cellI];
    scalar vl = vlc[cellI], vh = vhc[cellI];
    scalar nl = nlc[cellI], nh = nhc[cellI];
    scalar jfac = jfc[cellI], vfac = vfc[cellI], nfac = nfc[cellI];
    scalar T00 = Ttable[jl][vl][nl]*(1-jfac)+Ttable[jh][vl][nl]*jfac;
    scalar T01 = Ttable[jl][vl][nh]*(1-jfac)+Ttable[jh][vl][nh]*jfac;
    scalar T10 = Ttable[jl][vh][nl]*(1-jfac)+Ttable[jh][vh][nl]*jfac;
    scalar T11 = Ttable[jl][vh][nh]*(1-jfac)+Ttable[jh][vh][nh]*jfac;
    scalar T0 = T00*(1-vfac)+T10*vfac;
    scalar T1 = T01*(1-vfac)+T11*vfac;
    T[cellI] = T0*(1-nfac)+T1*nfac;
    scalar R00 = Rtable[jl][vl][nl]*(1-jfac)+Rtable[jh][vl][nl]*jfac;
    scalar R01 = Rtable[jl][vl][nh]*(1-jfac)+Rtable[jh][vl][nh]*jfac;
    scalar R10 = Rtable[jl][vh][nl]*(1-jfac)+Rtable[jh][vh][nl]*jfac;
    scalar R11 = Rtable[jl][vh][nh]*(1-jfac)+Rtable[jh][vh][nh]*jfac;
    scalar R0 = R00*(1-vfac)+R10*vfac;
    scalar R1 = R01*(1-vfac)+R11*vfac;
    scalar Rgas = R0*(1-nfac)+R1*nfac;
    rho[cellI] = p[cellI]/(Rgas*T[cellI]);
 }
 //correct processor boundary value of T and rho
 T.correctBoundaryConditions();
 rho.correctBoundaryConditions();
--- a/applications/solvers/combustion/SLFMFoam/updateYi.H
+++ b/applications/solvers/combustion/SLFMFoam/updateYi.H
@ -1,56 +0,0 @@
 //force the saving of the old-time values
 //important!!
 forAll(rho, cellI)
 {
    scalar& Nst = Neta[lowerN][cellI];
    //find Nst index
    if(Nst < NstList.first())
    {
        nlc[cellI] = 0;
        nfc[cellI] = 0;
    }
    else if(Nst > NstList[NstList.size()-2])
    {
        nhc[cellI] = NstList.size()-2; //just before extinction
        nfc[cellI] = 1;
    }
    else
    {
        nlc[cellI] = label( interpolateXY(Nst, NstList, NstIndex) );
        nhc[cellI] = nlc[cellI]+1;
        nfc[cellI] = (Nst-NstList[nlc[cellI]])
                     /(NstList[nhc[cellI]]-NstList[nlc[cellI]]);
    }
    //update Favre mean temp. and rho(=1/RT)
    scalar jl = jlc[cellI], jh = jhc[cellI];
    scalar vl = vlc[cellI], vh = vhc[cellI];
    scalar nl = nlc[cellI], nh = nhc[cellI];
    scalar jfac = jfc[cellI], vfac = vfc[cellI], nfac = nfc[cellI];
    forAll(Y, yi)
    {
        scalar y00 = Ytable[yi][jl][vl][nl]*(1-jfac)+Ytable[yi][jh][vl][nl]*jfac;
        scalar y01 = Ytable[yi][jl][vl][nh]*(1-jfac)+Ytable[yi][jh][vl][nh]*jfac;
        scalar y10 = Ytable[yi][jl][vh][nl]*(1-jfac)+Ytable[yi][jh][vh][nl]*jfac;
        scalar y11 = Ytable[yi][jl][vh][nh]*(1-jfac)+Ytable[yi][jh][vh][nh]*jfac;
        scalar y0 = y00*(1-vfac)+y10*vfac;
        scalar y1 = y01*(1-vfac)+y11*vfac;
        scalar ygas = y0*(1-nfac)+y1*nfac;
        Y[yi][cellI] = ygas;
    }
 }
 //correct processor boundary value of T and rho
 forAll(Y, yi)
 {
    Y[yi].correctBoundaryConditions();
 }
--- a/applications/solvers/lagrangian/ThermoMPPICFoam/EEqn.H
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/EEqn.H
@ -1,40 +0,0 @@
 {
    volScalarField& he = thermo.he();
    fvScalarMatrix EEqn
    (
        fvm::ddt(alphac,rho, he) + mvConvection->fvmDiv(alphacPhi, he)
      + fvc::ddt(alphac,rho, K) + fvc::div(alphacPhi, K)
      + (
            he.name() == "e"
          ? fvc::div
            (
                fvc::absolute(phi/fvc::interpolate(rho), U),
                p,
                "div(phiv,p)"
            )
          : -dpdt
        )
      - fvm::laplacian(alphac*turbulence->alphaEff(), he)
     ==
        alphac*rho*(U&g)
      + parcels.Sh(he)
      + radiation->Sh(thermo)
      + combustion->Sh()
      + fvOptions(alphac,rho, he)
    );
    EEqn.relax();
    fvOptions.constrain(EEqn);
    EEqn.solve();
    fvOptions.correct(he);
    thermo.correct();
    radiation->correct();
    Info<< "T gas min/max   = " << min(T).value() << ", "
        << max(T).value() << endl;
 }
--- a/applications/solvers/lagrangian/ThermoMPPICFoam/Make/files
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/Make/files
@ -1,3 +0,0 @@
 ThermoMPPICFoam.C
 EXE = $(FOAM_APPBIN)/ThermoMPPICFoam
--- a/applications/solvers/lagrangian/ThermoMPPICFoam/Make/options
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/Make/options
@ -1,56 +0,0 @@
 EXE_INC = \
    -I. \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
    -I$(LIB_SRC)/ODE/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
    -I$(LIB_SRC)/combustionModels/lnInclude \
    -I$(LIB_SRC)/sampling/lnInclude \
    -I$(FOAM_SOLVERS)/combustion/reactingFoam
 EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lliquidProperties \
    -lliquidMixtureProperties \
    -lsolidProperties \
    -lsolidMixtureProperties \
    -lthermophysicalFunctions \
    -lreactionThermophysicalModels \
    -lSLGThermo \
    -lchemistryModel \
    -lradiationModels \
    -lODE \
    -lregionModels \
    -lsurfaceFilmModels \
    -lcombustionModels \
    -lfvOptions \
    -lsampling
--- a/applications/solvers/lagrangian/ThermoMPPICFoam/ThermoMPPICFoam.C
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/ThermoMPPICFoam.C
@ -1,147 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-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/>.
 Application
    ThermoMPPICFoam
 Description
    Transient solver for compressible, turbulent flow with a reacting,
    multiphase particle cloud, and optional sources/constraints.
 \*---------------------------------------------------------------------------*/
 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
 #include "basicThermoKinematicMPPICCloud.H"
 #include "rhoCombustionModel.H"
 #include "radiationModel.H"
 #include "fvOptions.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 int main(int argc, char *argv[])
 {
    #include "postProcess.H"
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createControl.H"
    #include "createTimeControls.H"
    #include "createRDeltaT.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
    #include "createFvOptions.H"
    #include "initContinuityErrs.H"
    turbulence->validate();
    #include "compressibleCourantNo.H"
    #include "setInitialDeltaT.H"
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    Info<< "\nStarting time loop\n" << endl;
    while (runTime.run())
    {
        #include "readTimeControls.H"
        #include "compressibleCourantNo.H"
        #include "setDeltaT.H"
        runTime++;
        Info<< "Time = " << runTime.timeName() << nl << endl;
        parcels.evolve();
        // Update continuous phase volume fraction field
        alphac = max(1.0 - parcels.theta(), alphacMin);
        alphac.correctBoundaryConditions();
        alphacf = fvc::interpolate(alphac);
        alphacPhi = alphacf*phi;
        fvVectorMatrix cloudSU(parcels.SU(U));
        volVectorField cloudVolSUSu
        (
            IOobject
            (
                "cloudVolSUSu",
                runTime.timeName(),
                mesh
            ),
            mesh,
            dimensionedVector
            (
                "0",
                cloudSU.dimensions()/dimVolume,
                vector::zero
            ),
            zeroGradientFvPatchVectorField::typeName
        );
        cloudVolSUSu.primitiveFieldRef() = -cloudSU.source()/mesh.V();
        cloudVolSUSu.correctBoundaryConditions();
        cloudSU.source() = vector::zero;
        #include "rhoEqn.H"
        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UEqn.H"
            #include "YEqn.H"
            #include "EEqn.H"
            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }
            if (pimple.turbCorr())
            {
                turbulence->correct();
            }
        }
        rho = thermo.rho();
        runTime.write();
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }
    Info<< "End\n" << endl;
    return 0;
 }
 // ************************************************************************* //
--- a/applications/solvers/lagrangian/ThermoMPPICFoam/UEqn.H
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/UEqn.H
@ -1,43 +0,0 @@
    fvVectorMatrix UEqn
    (
        fvm::ddt(alphac, rho, U)
      + fvm::div(alphacPhi, U)
      - fvm::laplacian(alphac*turbulence->muEff(),U)
      - fvc::div(alphac*turbulence->muEff()*dev2(fvc::grad(U)().T()))
      - fvm::Sp(fvc::ddt(alphac,rho) + fvc::div(alphacPhi), U)
     ==
      cloudSU
      + fvOptions(alphac,rho, U)
    );
    UEqn.relax();
    volScalarField rAU(1.0/UEqn.A());
    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
    surfaceScalarField phicForces
    (
        (fvc::interpolate(rAU*cloudVolSUSu) & mesh.Sf())
        + alphacf*rhorAUf*(g & mesh.Sf())
    );
    fvOptions.constrain(UEqn);
    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                phicForces/rAUf
            )
          - fvc::grad(p)
        );
        fvOptions.correct(U);
        K = 0.5*magSqr(U);
    }
--- a/applications/solvers/lagrangian/ThermoMPPICFoam/YEqn.H
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/YEqn.H
@ -1,53 +0,0 @@
 tmp<fv::convectionScheme<scalar>> mvConvection
 (
    fv::convectionScheme<scalar>::New
    (
        mesh,
        fields,
        phi,
        mesh.divScheme("div(phi,Yi_h)")
    )
 );
 {
    combustion->correct();
    dQ = combustion->dQ();
    label inertIndex = -1;
    volScalarField Yt(0.0*Y[0]);
    forAll(Y, i)
    {
        if (Y[i].name() != inertSpecie)
        {
            volScalarField& Yi = Y[i];
            fvScalarMatrix YEqn
            (
                fvm::ddt(alphac,rho, Yi)
              + mvConvection->fvmDiv(alphacPhi, Yi)
              - fvm::laplacian(alphac*turbulence->muEff(), Yi)
             ==
                combustion->R(Yi)
              + fvOptions(alphac,rho, Yi)
            );
            YEqn.relax();
            fvOptions.constrain(YEqn);
            YEqn.solve(mesh.solver("Yi"));
            fvOptions.correct(Yi);
            Yi.max(0.0);
            Yt += Yi;
        }
        else
        {
            inertIndex = i;
        }
    }
    Y[inertIndex] = scalar(1) - Yt;
    Y[inertIndex].max(0.0);
 }
--- a/applications/solvers/lagrangian/ThermoMPPICFoam/createClouds.H
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/createClouds.H
@ -1,26 +0,0 @@
 Info<< "\nConstructing reacting cloud" << endl;
 basicThermoKinematicMPPICCloud parcels
 (
    "reactingCloud",
    rho,
    U,
    slgThermo.thermo().mu(),
    g
 );
 // Continuouse phase volume fraction lower limit
 scalar alphacMin
 (
    1.0
    - readScalar
    (
        parcels.particleProperties().subDict("constantProperties").lookup("alphaMax")
    )
 );
 // Update alphac from the particle locations
 alphac = max(1.0 - parcels.theta(), alphacMin);
 alphac.correctBoundaryConditions();
 alphacf = fvc::interpolate(alphac);
 alphacPhi = alphacf*phi;
--- a/applications/solvers/lagrangian/ThermoMPPICFoam/createFieldRefs.H
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/createFieldRefs.H
@ -1,2 +0,0 @@
 const volScalarField& T = thermo.T();
 const volScalarField& psi = thermo.psi();
--- a/applications/solvers/lagrangian/ThermoMPPICFoam/createFields.H
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/createFields.H
@ -1,158 +0,0 @@
 #include "readGravitationalAcceleration.H"
 Info<< "Creating combustion model\n" << endl;
 autoPtr<combustionModels::rhoCombustionModel> combustion
 (
    combustionModels::rhoCombustionModel::New(mesh)
 );
 rhoReactionThermo& thermo = combustion->thermo();
 thermo.validate(args.executable(), "h", "e");
 SLGThermo slgThermo(mesh, thermo);
 basicSpecieMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
 const word inertSpecie(thermo.lookup("inertSpecie"));
 if (!composition.contains(inertSpecie))
 {
    FatalErrorInFunction
        << "Specified inert specie '" << inertSpecie << "' not found in "
        << "species list. Available species:" << composition.species()
        << exit(FatalError);
 }
 volScalarField& p = thermo.p();
 volScalarField rho
 (
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    thermo.rho()
 );
 Info<< "\nReading field U\n" << endl;
 volVectorField U
 (
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
 );
 #include "compressibleCreatePhi.H"
 mesh.setFluxRequired(p.name());
 dimensionedScalar rhoMax
 (
    dimensionedScalar::lookupOrDefault
    (
        "rhoMax",
        pimple.dict(),
        dimDensity,
        GREAT
    )
 );
 dimensionedScalar rhoMin
 (
    dimensionedScalar::lookupOrDefault
    (
        "rhoMin",
        pimple.dict(),
        dimDensity,
        0
    )
 );
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
    compressible::turbulenceModel::New
    (
        rho,
        U,
        phi,
        thermo
    )
 );
 // Set the turbulence into the combustion model
 combustion->setTurbulence(turbulence());
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt
 (
    IOobject
    (
        "dpdt",
        runTime.timeName(),
        mesh
    ),
    mesh,
    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
 );
 Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
 Info<< "Creating multi-variate interpolation scheme\n" << endl;
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 forAll(Y, i)
 {
    fields.add(Y[i]);
 }
 fields.add(thermo.he());
 volScalarField dQ
 (
    IOobject
    (
        "dQ",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );
 volScalarField alphac
 (
    IOobject
    (
        "alphac",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
        mesh,
        dimensionedScalar("0", dimless, 1)
 );
 surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));
 surfaceScalarField alphacPhi("alphacPhi", alphacf*phi);
 volScalarField alphacrho("alphacrho", alphac*rho);
 #include "createMRF.H"
 #include "createRadiationModel.H"
 #include "createClouds.H"
--- a/applications/solvers/lagrangian/ThermoMPPICFoam/pEqn.H
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/pEqn.H
@ -1,81 +0,0 @@
 {
    rho = thermo.rho();
    // Thermodynamic density needs to be updated by psi*d(p) after the
    // pressure solution - done in 2 parts. Part 1:
    thermo.rho() -= psi*p;
    rhorAUf= fvc::interpolate(rho*rAU);
    surfaceScalarField rhof("rhof", fvc::interpolate(rho));
    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        (
         (fvc::interpolate(rho*HbyA) & mesh.Sf())
         + rhorAUf*fvc::ddtCorr(rho, U, phi)
         + phicForces*rhof
        )
    );
       MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
    fvScalarMatrix pDDtEqn
    (
        fvc::ddt(alphac,rho) + alphac*psi*correction(fvm::ddt(p))
      + fvc::div(alphacf*phiHbyA)
     ==
        fvOptions(alphac*psi, p, rho.name())
    );
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            pDDtEqn
          - fvm::laplacian(alphacf*rhorAUf, p)
        );
        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
        if (pimple.finalNonOrthogonalIter())
        {
           phi = phiHbyA + pEqn.flux()/alphacf;
           alphacPhi = alphacf*phi;
           p.relax();
           thermo.rho() += psi*p;
           #include "rhoEqn.H"
           #include "compressibleContinuityErrs.H"
           U = HbyA + rAU*fvc::reconstruct(phicForces/rAUf) - rAU*fvc::grad(p);
         }
    }
    U.correctBoundaryConditions();
    fvOptions.correct(U);
    K = 0.5*magSqr(U);
    if (thermo.dpdt())
    {
        dpdt = fvc::ddt(p);
    }
    rho = thermo.rho();
    rho = max(rho, rhoMin);
    rho = min(rho, rhoMax);
    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
 }
--- a/applications/solvers/lagrangian/ThermoMPPICFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/rhoEqn.H
@ -1,52 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 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/>.
 Global
    rhoEqn
 Description
    Solve the continuity for density.
 \*---------------------------------------------------------------------------*/
 {
    fvScalarMatrix rhoEqn
    (
        fvm::ddt(alphac,rho)
      + fvc::div(alphacPhi)
      ==
       fvOptions(alphac,rho)
    );
    rhoEqn.solve();
    fvOptions.correct(rho);
    Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value()
        << endl;
    rho = max(rho, rhoMin);
    rho = min(rho, rhoMax);
 }
 // ************************************************************************* //
--- a/applications/solvers/lagrangian/ThermoMPPICFoam/setRDeltaT.H
+++ b/applications/solvers/lagrangian/ThermoMPPICFoam/setRDeltaT.H
@ -1,134 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-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/>.
 \*---------------------------------------------------------------------------*/
 {
    volScalarField& rDeltaT = trDeltaT.ref();
    const dictionary& pimpleDict = pimple.dict();
    // Maximum flow Courant number
    scalar maxCo(readScalar(pimpleDict.lookup("maxCo")));
    // Maximum time scale
    scalar maxDeltaT(pimpleDict.lookupOrDefault<scalar>("maxDeltaT", GREAT));
    // Smoothing parameter (0-1) when smoothing iterations > 0
    scalar rDeltaTSmoothingCoeff
    (
        pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
    );
    // Damping coefficient (1-0)
    scalar rDeltaTDampingCoeff
    (
        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 0.2)
    );
    // Maximum change in cell temperature per iteration
    // (relative to previous value)
    scalar alphaTemp(pimpleDict.lookupOrDefault("alphaTemp", 0.05));
    Info<< "Time scales min/max:" << endl;
    // Cache old reciprocal time scale field
    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
    // Flow time scale
    {
        rDeltaT.ref() =
        (
            fvc::surfaceSum(mag(phi))()()
           /((2*maxCo)*mesh.V()*rho())
        );
        // Limit the largest time scale
        rDeltaT.max(1/maxDeltaT);
        Info<< "    Flow        = "
            << gMin(1/rDeltaT.primitiveField()) << ", "
            << gMax(1/rDeltaT.primitiveField()) << endl;
    }
    // Reaction source time scale
    {
        volScalarField::Internal rDeltaTT
        (
            mag
            (
               parcels.hsTrans()/(mesh.V()*runTime.deltaT())
             + combustion->Sh()()
            )
           /(
               alphaTemp
              *rho()
              *thermo.Cp()()()
              *T()
           )
        );
        Info<< "    Temperature = "
            << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
            << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
        rDeltaT.ref() = max
        (
            rDeltaT(),
            rDeltaTT
        );
    }
    // Update tho boundary values of the reciprocal time-step
    rDeltaT.correctBoundaryConditions();
    // Spatially smooth the time scale field
    if (rDeltaTSmoothingCoeff < 1.0)
    {
        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
    }
    // Limit rate of change of time scale
    // - reduce as much as required
    // - only increase at a fraction of old time scale
    if
    (
        rDeltaTDampingCoeff < 1.0
     && runTime.timeIndex() > runTime.startTimeIndex() + 1
    )
    {
        rDeltaT = max
        (
            rDeltaT,
            (scalar(1.0) - rDeltaTDampingCoeff)*rDeltaT0
        );
    }
    Info<< "    Overall     = "
        << gMin(1/rDeltaT.primitiveField())
        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
 }
 // ************************************************************************* //
--- a/applications/solvers/lagrangian/simpleCoalcombustionFoam/EEqn.H
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/EEqn.H
@ -1,32 +0,0 @@
 {
    volScalarField& he = thermo.he();
    fvScalarMatrix EEqn
    (
        mvConvection->fvmDiv(phi, he)
      + (
            he.name() == "e"
          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
        )
      - fvm::laplacian(turbulence->alphaEff(), he)
     ==
        parcels.Sh(he)
      + radiation->Sh(thermo)
      + combustion->Sh()
      + fvOptions(rho, he)
    );
    EEqn.relax();
    fvOptions.constrain(EEqn);
    EEqn.solve();
    fvOptions.correct(he);
    thermo.correct();
    radiation->correct();
    Info<< "T gas min/max   = " << min(T).value() << ", "
        << max(T).value() << endl;
 }
--- a/applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/files
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/files
@ -1,3 +0,0 @@
 simpleCoalcombustionFoam.C
 EXE = $(FOAM_APPBIN)/simpleCoalcombustionFoam
--- a/applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/options
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/Make/options
@ -1,57 +0,0 @@
 EXE_INC = \
    -I. \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I${LIB_SRC}/meshTools/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
    -I$(LIB_SRC)/ODE/lnInclude \
    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
    -I$(LIB_SRC)/combustionModels/lnInclude \
    -I$(LIB_SRC)/sampling/lnInclude \
    -I$(FOAM_SOLVERS)/combustion/reactingFoam
 EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
    -llagrangian \
    -llagrangianIntermediate \
    -lcoalCombustion \
    -llagrangianTurbulence \
    -lspecie \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lliquidProperties \
    -lliquidMixtureProperties \
    -lsolidProperties \
    -lsolidMixtureProperties \
    -lthermophysicalFunctions \
    -lreactionThermophysicalModels \
    -lSLGThermo \
    -lchemistryModel \
    -lradiationModels \
    -lODE \
    -lregionModels \
    -lsurfaceFilmModels \
    -lcombustionModels \
    -lfvOptions \
    -lsampling
--- a/applications/solvers/lagrangian/simpleCoalcombustionFoam/UEqn.H
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/UEqn.H
@ -1,21 +0,0 @@
    MRF.correctBoundaryVelocity(U);
    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(U)
     ==
        rho()*g
      + parcels.SU(U)
      + fvOptions(rho, U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();
    UEqn.relax();
    fvOptions.constrain(UEqn);
    solve(UEqn == -fvc::grad(p));
    fvOptions.correct(U);
--- a/applications/solvers/lagrangian/simpleCoalcombustionFoam/YEqn.H
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/YEqn.H
@ -1,53 +0,0 @@
 tmp<fv::convectionScheme<scalar>> mvConvection
 (
    fv::convectionScheme<scalar>::New
    (
        mesh,
        fields,
        phi,
        mesh.divScheme("div(phi,Yi_h)")
    )
 );
 {
    combustion->correct();
    dQ = combustion->dQ();
    label inertIndex = -1;
    volScalarField Yt(0.0*Y[0]);
    forAll(Y, i)
    {
        if (Y[i].name() != inertSpecie)
        {
            volScalarField& Yi = Y[i];
            fvScalarMatrix YEqn
            (
                mvConvection->fvmDiv(phi, Yi)
              - fvm::laplacian(turbulence->muEff(), Yi)
             ==
                parcels.Srho(i)
              + combustion->R(Yi)
              + fvOptions(rho, Yi)
            );
            YEqn.relax();
            fvOptions.constrain(YEqn);
            YEqn.solve(mesh.solver("Yi"));
            fvOptions.correct(Yi);
            Yi.max(0.0);
            Yt += Yi;
        }
        else
        {
            inertIndex = i;
        }
    }
    Y[inertIndex] = scalar(1) - Yt;
    Y[inertIndex].max(0.0);
 }
--- a/applications/solvers/lagrangian/simpleCoalcombustionFoam/createClouds.H
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/createClouds.H
@ -1,8 +0,0 @@
 Info<< "\nConstructing reacting cloud" << endl;
 coalCloudList parcels
 (
    rho,
    U,
    g,
    slgThermo
 );
--- a/applications/solvers/lagrangian/simpleCoalcombustionFoam/createFieldRefs.H
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/createFieldRefs.H
@ -1,2 +0,0 @@
 const volScalarField& T = thermo.T();
 const volScalarField& psi = thermo.psi();
--- a/applications/solvers/lagrangian/simpleCoalcombustionFoam/createFields.H
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/createFields.H
@ -1,123 +0,0 @@
 #include "readGravitationalAcceleration.H"
 Info<< "Creating combustion model\n" << endl;
 autoPtr<combustionModels::rhoCombustionModel> combustion
 (
    combustionModels::rhoCombustionModel::New(mesh)
 );
 rhoReactionThermo& thermo = combustion->thermo();
 thermo.validate(args.executable(), "h", "e");
 SLGThermo slgThermo(mesh, thermo);
 basicSpecieMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
 const word inertSpecie(thermo.lookup("inertSpecie"));
 if (!composition.contains(inertSpecie))
 {
    FatalErrorInFunction
        << "Specified inert specie '" << inertSpecie << "' not found in "
        << "species list. Available species:" << composition.species()
        << exit(FatalError);
 }
 volScalarField& p = thermo.p();
 volScalarField rho
 (
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    thermo.rho()
 );
 Info<< "\nReading field U\n" << endl;
 volVectorField U
 (
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
 );
 #include "compressibleCreatePhi.H"
 mesh.setFluxRequired(p.name());
 dimensionedScalar rhoMax
 (
    dimensionedScalar::lookupOrDefault
    (
        "rhoMax",
        simple.dict(),
        dimDensity,
        GREAT
    )
 );
 dimensionedScalar rhoMin
 (
    dimensionedScalar::lookupOrDefault
    (
        "rhoMin",
        simple.dict(),
        dimDensity,
        0
    )
 );
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
    compressible::turbulenceModel::New
    (
        rho,
        U,
        phi,
        thermo
    )
 );
 // Set the turbulence into the combustion model
 combustion->setTurbulence(turbulence());
 Info<< "Creating multi-variate interpolation scheme\n" << endl;
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 forAll(Y, i)
 {
    fields.add(Y[i]);
 }
 fields.add(thermo.he());
 volScalarField dQ
 (
    IOobject
    (
        "dQ",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );
 #include "createMRF.H"
 #include "createRadiationModel.H"
 #include "createClouds.H"
--- a/applications/solvers/lagrangian/simpleCoalcombustionFoam/pEqn.H
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/pEqn.H
@ -1,57 +0,0 @@
 {
    // Thermodynamic density needs to be updated by psi*d(p) after the
    // pressure solution - done in 2 parts. Part 1:
    thermo.rho() -= psi*p;
    volScalarField rAU(1.0/UEqn.A());
    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
    tUEqn.clear();
    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        fvc::interpolate(rho)*fvc::flux(HbyA)
    );
    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
    while (simple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvc::div(phiHbyA)
          - fvm::laplacian(rhorAUf, p)
         ==
            parcels.Srho()
          + fvOptions(psi, p, rho.name())
        );
        pEqn.solve();
        if (simple.finalNonOrthogonalIter())
        {
            phi = phiHbyA + pEqn.flux();
        }
    }
    p.relax();
    // Second part of thermodynamic density update
    thermo.rho() += psi*p;
    #include "compressibleContinuityErrs.H"
    U = HbyA - rAU*fvc::grad(p);
    U.correctBoundaryConditions();
    fvOptions.correct(U);
    rho = thermo.rho();
    rho = max(rho, rhoMin);
    rho = min(rho, rhoMax);
    rho.relax();
    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
 }
--- a/applications/solvers/lagrangian/simpleCoalcombustionFoam/simpleCoalcombustionFoam.C
+++ b/applications/solvers/lagrangian/simpleCoalcombustionFoam/simpleCoalcombustionFoam.C
@ -1,93 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013-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/>.
 Application
    simpleCoalcombustionFoam
 Description
    Steady state solver for compressible, turbulent flow with reacting,
    multiphase particle clouds and optional sources/constraints.
 \*---------------------------------------------------------------------------*/
 #include "fvCFD.H"
 #include "turbulentFluidThermoModel.H"
 #include "coalCloudList.H"
 #include "rhoCombustionModel.H"
 #include "radiationModel.H"
 #include "IOporosityModelList.H"
 #include "fvOptions.H"
 #include "SLGThermo.H"
 #include "simpleControl.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 int main(int argc, char *argv[])
 {
    #include "postProcess.H"
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createControl.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
    #include "createFvOptions.H"
    #include "initContinuityErrs.H"
    turbulence->validate();
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    Info<< "\nStarting time loop\n" << endl;
    while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;
        parcels.evolve();
        // --- Pressure-velocity SIMPLE corrector loop
        {
            #include "UEqn.H"
            #include "YEqn.H"
            #include "EEqn.H"
            #include "pEqn.H"
        }
        turbulence->correct();
        runTime.write();
        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }
    Info<< "End\n" << endl;
    return 0;
 }
 // ************************************************************************* //
--- a/src/combustionModels/EDM/EDM.C
+++ b/src/combustionModels/EDM/EDM.C
@ -1,227 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013-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 "EDM.H"
 #include "fvmSup.H"
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 template<class Type>
 Foam::combustionModels::EDM<Type>::EDM
 (
    const word& modelType,
    const fvMesh& mesh,
    const word& phaseName
 )
 :
    Type(modelType, mesh, phaseName),
    finiteRate_ (this->coeffs().lookupOrDefault("finiteRate", false)),
    A_ (this->coeffs().lookupOrDefault("A", 4.0)),
    B_ (this->coeffs().lookupOrDefault("B", 0.5))
 {
    if (finiteRate_)
    {
        Info<< "    using Finite-rate/Eddy Dissipation Model" << endl;
    }
    else
    {
        Info<< "    using Eddy Dissipation Model" << endl;
    }
    Info<< "    A = " << A_ << endl;
    Info<< "    B = " << B_ << endl;
 }
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 template<class Type>
 Foam::combustionModels::EDM<Type>::~EDM()
 {}
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 template<class Type>
 Foam::tmp<Foam::volScalarField>
 Foam::combustionModels::EDM<Type>::tc() const
 {
    return this->chemistryPtr_->tc();
 }
 template<class Type>
 void Foam::combustionModels::EDM<Type>::correct()
 {
    if (this->active())
    {
        tmp<volScalarField> tk(this->turbulence().k());
        const volScalarField& k = tk();
        tmp<volScalarField> tepsilon(this->turbulence().epsilon());
        const volScalarField& epsilon = tepsilon();
        this->chemistryPtr_->edm(epsilon/k, A_, B_);
        if (finiteRate_)
        {
            // Copy edm rate stored in chemistryPtr_->RR(i)
            PtrList<DimensionedField<scalar, volMesh>> edmRates
            (
                this->thermo().composition().Y().size()
            );
            forAll(edmRates, i)
            {
                edmRates.set
                (
                    i,
                    new DimensionedField<scalar, volMesh>
                    (
                        this->chemistryPtr_->RR(i)
                    )
                );
            }
            this->chemistryPtr_->calculate();
            forAll(edmRates, i)
            {
                DimensionedField<scalar, volMesh> &finiteRatei
                (
                    this->chemistryPtr_->RR(i)
                );
                forAll(finiteRatei, celli)
                {
                    finiteRatei[celli] = min(finiteRatei[celli], edmRates[i][celli]);
                    finiteRatei[celli] =
                        finiteRatei[celli]*edmRates[i][celli]
                        / (finiteRatei[celli] + edmRates[i][celli] + SMALL);
                }
            }
        }
    }
 }
 template<class Type>
 Foam::tmp<Foam::fvScalarMatrix>
 Foam::combustionModels::EDM<Type>::R(volScalarField& Y) const
 {
    tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));
    fvScalarMatrix& Su = tSu.ref();
    if (this->active())
    {
        const label specieI =
            this->thermo().composition().species()[Y.member()];
        Su += this->chemistryPtr_->RR(specieI);
    }
    return tSu;
 }
 template<class Type>
 Foam::tmp<Foam::volScalarField>
 Foam::combustionModels::EDM<Type>::dQ() const
 {
    tmp<volScalarField> tdQ
    (
        new volScalarField
        (
            IOobject
            (
                IOobject::groupName(typeName + ":dQ", this->phaseName_),
                this->mesh().time().timeName(),
                this->mesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            this->mesh(),
            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
        )
    );
    if (this->active())
    {
        tdQ.ref() = this->chemistryPtr_->dQ();
    }
    return tdQ;
 }
 template<class Type>
 Foam::tmp<Foam::volScalarField>
 Foam::combustionModels::EDM<Type>::Sh() const
 {
    tmp<volScalarField> tSh
    (
        new volScalarField
        (
            IOobject
            (
                IOobject::groupName(typeName + ":Sh", this->phaseName_),
                this->mesh().time().timeName(),
                this->mesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            this->mesh(),
            dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
        )
    );
    if (this->active())
    {
        tSh.ref() = this->chemistryPtr_->Sh();
    }
    return tSh;
 }
 template<class Type>
 bool Foam::combustionModels::EDM<Type>::read()
 {
    if (Type::read())
    {
        this->coeffs().lookup("finiteRate") >> finiteRate_;
        return true;
    }
    else
    {
        return false;
    }
 }
 // ************************************************************************* //
--- a/src/combustionModels/EDM/EDM.H
+++ b/src/combustionModels/EDM/EDM.H
@ -1,147 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013-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::combustionModels::EDM
 Description
    Eddy dissipation model with finite-rate chemistry.
                                                                    28.Feb.2017
                                                        Combustion Lab. POSTECH
                                                                      Karam Han
                                              (Version upgraded by Jinwoo Park)
 SourceFiles
    EDM.C
 \*---------------------------------------------------------------------------*/
 #ifndef EDM_H
 #define EDM_H
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 namespace Foam
 {
 namespace combustionModels
 {
 /*---------------------------------------------------------------------------*\
                            Class EDM Declaration
 \*---------------------------------------------------------------------------*/
 template<class Type>
 class EDM
 :
    public Type
 {
    // Private data
        //- finite-rate/EDM using Arrhenius formula
        bool finiteRate_;
        //- EDM Parameter A
        scalar A_;
        //- EDM Parameter B
        scalar B_;
 protected:
    // Protected Member Functions
        //- Return the chemical time scale
        tmp<volScalarField> tc() const;
 private:
    // Private Member Functions
        //- Disallow copy construct
        EDM(const EDM&);
        //- Disallow default bitwise assignment
        void operator=(const EDM&);
 public:
    //- Runtime type information
    TypeName("EDM");
    // Constructors
        //- Construct from components
        EDM
        (
            const word& modelType,
            const fvMesh& mesh,
            const word& phaseName
        );
    //- Destructor
    virtual ~EDM();
    // Member Functions
        // Evolution
            //- Correct combustion rate
            virtual void correct();
            //- Fuel consumption rate matrix.
            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
            //- Heat release rate calculated from fuel consumption rate matrix
            virtual tmp<volScalarField> dQ() const;
            //-  Return source for enthalpy equation [kg/m/s3]
            virtual tmp<volScalarField> Sh() const;
    // IO
            //- Update properties from given dictionary
            virtual bool read();
 };
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 } // End namespace combustionModels
 } // End namespace Foam
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #ifdef NoRepository
    #include "EDM.C"
 #endif
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #endif
 // ************************************************************************* //
--- a/src/combustionModels/EDM/EDMs.C
+++ b/src/combustionModels/EDM/EDMs.C
@ -1,38 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013-2015 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 "makeCombustionTypes.H"
 #include "psiChemistryCombustion.H"
 #include "rhoChemistryCombustion.H"
 #include "EDM.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 makeCombustionTypes(EDM, psiChemistryCombustion, psiCombustionModel);
 makeCombustionTypes(EDM, rhoChemistryCombustion, rhoCombustionModel);
 // ************************************************************************* //
--- a/src/combustionModels/Make/files
+++ b/src/combustionModels/Make/files
@ -13,8 +13,6 @@ rhoCombustionModel/rhoChemistryCombustion/rhoChemistryCombustion.C
 diffusion/diffusions.C
 infinitelyFastChemistry/infinitelyFastChemistrys.C
 EDM/EDMs.C
 PaSR/PaSRs.C
 laminar/laminars.C
--- a/src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H
+++ b/src/lagrangian/coalCombustion/include/makeCoalParcelSurfaceReactionModels.H
@ -34,7 +34,6 @@ License
 #include "COxidationKineticDiffusionLimitedRate.H"
 #include "COxidationHurtMitchell.H"
 #include "COxidationMurphyShaddix.H"
 #include "COxidationKD_CO.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -48,7 +47,6 @@ License
    );                                                                         \
    makeSurfaceReactionModelType(COxidationIntrinsicRate, CloudType);          \
    makeSurfaceReactionModelType(COxidationHurtMitchell, CloudType);           \
    makeSurfaceReactionModelType(COxidationKD_CO, CloudType);                  \
    makeSurfaceReactionModelType(COxidationMurphyShaddix, CloudType);
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKD_CO/COxidationKD_CO.C
+++ b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKD_CO/COxidationKD_CO.C
@ -1,178 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 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 "COxidationKD_CO.H"
 #include "mathematicalConstants.H"
 #include "thermodynamicConstants.H"
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 template<class CloudType>
 Foam::COxidationKD_CO<CloudType>::
 COxidationKD_CO
 (
    const dictionary& dict,
    CloudType& owner
 )
 :
    SurfaceReactionModel<CloudType>(dict, owner, typeName),
    Sb_(readScalar(this->coeffDict().lookup("Sb"))),
    C1_(readScalar(this->coeffDict().lookup("C1"))),
    C2_(readScalar(this->coeffDict().lookup("C2"))),
    E_(readScalar(this->coeffDict().lookup("E"))),
    CsLocalId_(-1),
    O2GlobalId_(owner.composition().carrierId("O2")),
    COGlobalId_(owner.composition().carrierId("CO")),
    WC_(0.0),
    WO2_(0.0),
    HcCO_(0.0)
 {
    // Determine Cs ids
    label idSolid = owner.composition().idSolid();
    CsLocalId_ = owner.composition().localId(idSolid, "C");
    // Set local copies of thermo properties
    WO2_ = owner.thermo().carrier().W(O2GlobalId_);
    const scalar WCO = owner.thermo().carrier().W(COGlobalId_);
    WC_ = WCO - 0.5*WO2_;
    HcCO_ = owner.thermo().carrier().Hc(COGlobalId_);
    const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
    const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
    Info<< "    C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
 }
 template<class CloudType>
 Foam::COxidationKD_CO<CloudType>::
 COxidationKD_CO
 (
    const COxidationKD_CO<CloudType>& srm
 )
 :
    SurfaceReactionModel<CloudType>(srm),
    Sb_(srm.Sb_),
    C1_(srm.C1_),
    C2_(srm.C2_),
    E_(srm.E_),
    CsLocalId_(srm.CsLocalId_),
    O2GlobalId_(srm.O2GlobalId_),
    COGlobalId_(srm.COGlobalId_),
    WC_(srm.WC_),
    WO2_(srm.WO2_),
    HcCO_(srm.HcCO_)
 {}
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 template<class CloudType>
 Foam::COxidationKD_CO<CloudType>::
 ~COxidationKD_CO()
 {}
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 template<class CloudType>
 Foam::scalar Foam::COxidationKD_CO<CloudType>::calculate
 (
    const scalar dt,
    const label cellI,
    const scalar d,
    const scalar T,
    const scalar Tc,
    const scalar pc,
    const scalar rhoc,
    const scalar mass,
    const scalarField& YGas,
    const scalarField& YLiquid,
    const scalarField& YSolid,
    const scalarField& YMixture,
    const scalar N,
    scalarField& dMassGas,
    scalarField& dMassLiquid,
    scalarField& dMassSolid,
    scalarField& dMassSRCarrier
 ) const
 {
    // Fraction of remaining combustible material
    const label idSolid = CloudType::parcelType::SLD;
    const scalar fComb = YMixture[idSolid]*YSolid[CsLocalId_];
    // Surface combustion active combustible fraction is consumed
    if (fComb < SMALL)
    {
        return 0.0;
    }
    const SLGThermo& thermo = this->owner().thermo();
    // Local mass fraction of O2 in the carrier phase
    const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[cellI];
    // Diffusion rate coefficient
    const scalar D0 = C1_/d*pow(0.5*(T + Tc), 0.75);
    // Kinetic rate
    const scalar Rk = C2_*exp(-E_/(constant::thermodynamic::RR*T)); //
    // Particle surface area
    const scalar Ap = constant::mathematical::pi*sqr(d);
    // Change in C mass [kg]
    scalar dmC = Ap*rhoc*constant::thermodynamic::RR*Tc*YO2/WO2_*D0*Rk/(D0 + Rk)*dt;
    // Limit mass transfer by availability of C
    dmC = min(mass*fComb, dmC);
    // Molar consumption
    const scalar dOmega = dmC/WC_;
    // Change in O2 mass [kg]
    const scalar dmO2 = dOmega*Sb_*WO2_;
    // Mass of newly created CO [kg]
    const scalar dmCO = dOmega*(WC_ + Sb_*WO2_);
    // Update local particle C mass
    dMassSolid[CsLocalId_] += dOmega*WC_;
    // Update carrier O2 and CO mass
    dMassSRCarrier[O2GlobalId_] -= dmO2;
    dMassSRCarrier[COGlobalId_] += dmCO;
    const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
    // carrier sensible enthalpy exchange handled via change in mass
    // Heat of reaction [J]
 //Info<<"Heat of reaction(Char)"<<tab<<(dmC*HsC - dmCO*HcCO_)/dmC<<endl;
    return dmC*HsC - dmCO*HcCO_;
 }
 // ************************************************************************* //
--- a/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKD_CO/COxidationKD_CO.H
+++ b/src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKD_CO/COxidationKD_CO.H
@ -1,176 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 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
    COxidationKD_CO
 Description
    Kinetic/diffusion limited rate surface reaction model for coal parcels.
    Limited to:
        C(s) + Sb*O2 -> CO
    where Sb is the stoichiometry of the reaction
 \*---------------------------------------------------------------------------*/
 #ifndef COxidationKD_CO_H
 #define COxidationKD_CO_H
 #include "SurfaceReactionModel.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 namespace Foam
 {
 // Forward class declarations
 template<class CloudType>
 class COxidationKD_CO;
 /*---------------------------------------------------------------------------*\
            Class COxidationKD_CO Declaration
 \*---------------------------------------------------------------------------*/
 template<class CloudType>
 class COxidationKD_CO
 :
    public SurfaceReactionModel<CloudType>
 {
    // Private data
        // Model constants
            //- Stoichiometry of reaction
            const scalar Sb_;
            //- Mass diffusion limited rate constant, C1
            const scalar C1_;
            //- Kinetics limited rate pre-exponential constant, C2
            const scalar C2_;
            //- Kinetics limited rate activation energy
            const scalar E_;
        // Addressing
            //- Cs positions in global/local lists
            label CsLocalId_;
            //- O2 position in global list
            label O2GlobalId_;
            //- CO2 positions in global list
            label COGlobalId_;
        // Local copies of thermo properties
            //- Molecular weight of C [kg/kmol]
            scalar WC_;
            //- Molecular weight of O2 [kg/kmol]
            scalar WO2_;
            //- Formation enthalpy for CO2 [J/kg]
            scalar HcCO_;
 public:
    //- Runtime type information
    TypeName("COxidationKD_CO");
    // Constructors
        //- Construct from dictionary
        COxidationKD_CO
        (
            const dictionary& dict,
            CloudType& owner
        );
        //- Construct copy
        COxidationKD_CO
        (
            const COxidationKD_CO<CloudType>& srm
        );
        //- Construct and return a clone
        virtual autoPtr<SurfaceReactionModel<CloudType> > clone() const
        {
            return autoPtr<SurfaceReactionModel<CloudType> >
            (
                new COxidationKD_CO<CloudType>(*this)
            );
        }
    //- Destructor
    virtual ~COxidationKD_CO();
    // Member Functions
        //- Update surface reactions
        virtual scalar calculate
        (
            const scalar dt,
            const label cellI,
            const scalar d,
            const scalar T,
            const scalar Tc,
            const scalar pc,
            const scalar rhoc,
            const scalar mass,
            const scalarField& YGas,
            const scalarField& YLiquid,
            const scalarField& YSolid,
            const scalarField& YMixture,
            const scalar N,
            scalarField& dMassGas,
            scalarField& dMassLiquid,
            scalarField& dMassSolid,
            scalarField& dMassSRCarrier
        ) const;
 };
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 } // End namespace Foam
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #ifdef NoRepository
    #include "COxidationKD_CO.C"
 #endif
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #endif
 // ************************************************************************* //
--- a/src/lagrangian/intermediate/Make/files
+++ b/src/lagrangian/intermediate/Make/files
@ -50,12 +50,6 @@ $(KINEMATICMPPICPARCEL)/defineBasicKinematicMPPICParcel.C
 $(KINEMATICMPPICPARCEL)/makeBasicKinematicMPPICParcelSubmodels.C
 /* thermo kinematic MPPIC parcel sub-models */
 THERMOMPPICPARCEL=$(DERIVEDPARCELS)/basicThermoKinematicMPPICParcel
 $(THERMOMPPICPARCEL)/defineBasicThermoKinematicMPPICParcel.C
 $(THERMOMPPICPARCEL)/makeBasicThermoKinematicMPPICParcelSubmodels.C
 /* bolt-on models */
 RADIATION=submodels/addOns/radiation
 $(RADIATION)/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C
@ -213,88 +213,6 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
 }
 template<class CloudType>
 Foam::ThermoCloud<CloudType>::ThermoCloud
 (
    const word& cloudName,
    const volScalarField& rho,
    const volVectorField& U,
    const volScalarField& mu,
    const dimensionedVector& g,
    bool readFields
 )
 :
    CloudType
    (
        cloudName,
        rho,
        U,
        mu,
        g,
        false
    ),
    thermoCloud(),
    cloudCopyPtr_(NULL),
    constProps_(this->particleProperties()),
    thermo_(this->mesh().template lookupObject<SLGThermo>("SLGThermo")),
    T_(thermo_.thermo().T()),
    p_(thermo_.thermo().p()),
    heatTransferModel_(NULL),
    TIntegrator_(NULL),
    radiation_(false),
    radAreaP_(NULL),
    radT4_(NULL),
    radAreaPT4_(NULL),
    hsTrans_
    (
        new DimensionedField<scalar, volMesh>
        (
            IOobject
            (
                this->name() + ":hsTrans",
                this->db().time().timeName(),
                this->db(),
                IOobject::READ_IF_PRESENT,
                IOobject::AUTO_WRITE
            ),
            this->mesh(),
            dimensionedScalar("zero", dimEnergy, 0.0)
        )
    ),
    hsCoeff_
    (
        new DimensionedField<scalar, volMesh>
        (
            IOobject
            (
                this->name() + ":hsCoeff",
                this->db().time().timeName(),
                this->db(),
                IOobject::READ_IF_PRESENT,
                IOobject::AUTO_WRITE
            ),
            this->mesh(),
            dimensionedScalar("zero", dimEnergy/dimTemperature, 0.0)
        )
    )
 {
    if (this->solution().active())
    {
        setModels();
        if (readFields)
        {
            parcelType::readFields(*this);
        }
    }
    if (this->solution().resetSourcesOnStartup())
    {
        resetSourceTerms();
    }
 }
 template<class CloudType>
 Foam::ThermoCloud<CloudType>::ThermoCloud
 (
--- a/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H
@ -180,17 +180,6 @@ public:
            bool readFields = true
        );
        //- Construct given carrier gas fields
        ThermoCloud
        (
            const word& cloudName,
            const volScalarField& rho,
            const volVectorField& U,
            const volScalarField& mu,
            const dimensionedVector& g,
            bool readFields = true
        );
        //- Copy constructor with new name
        ThermoCloud(ThermoCloud<CloudType>& c, const word& name);
--- a/src/lagrangian/intermediate/clouds/derived/basicThermoKinematicMPPICCloud/basicThermoKinematicMPPICCloud.H
+++ b/src/lagrangian/intermediate/clouds/derived/basicThermoKinematicMPPICCloud/basicThermoKinematicMPPICCloud.H
@ -1,66 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013 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::basicThermoKinematicMPPICCloud
 Description
    Cloud class to introduce Thermo kinematic MPPIC parcels
 \*---------------------------------------------------------------------------*/
 #ifndef basicThermoKinematicMPPICCloud_H
 #define basicThermoKinematicMPPICCloud_H
 #include "Cloud.H"
 #include "KinematicCloud.H"
 #include "MPPICCloud.H"
 #include "basicThermoKinematicMPPICParcel.H"
 #include "ThermoCloud.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 namespace Foam
 {
    typedef MPPICCloud
    <
       ThermoCloud
       <
           KinematicCloud
           <
                Cloud
                <
                    basicThermoKinematicMPPICParcel
                >
           >
       >
    >
    basicThermoKinematicMPPICCloud;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #endif
 // ************************************************************************* //
--- a/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/basicThermoKinematicMPPICParcel.H
+++ b/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/basicThermoKinematicMPPICParcel.H
@ -1,63 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013-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::basicThermoKinematicMPPICParcel
 Description
    Definition of basic Thermo kinematic MPPIC parcel
 SourceFiles
    basicThermoKinematicMPPICParcel.H
 \*---------------------------------------------------------------------------*/
 #ifndef basicThermoKinematicMPPICParcel_H
 #define basicThermoKinematicMPPICParcel_H
 #include "contiguous.H"
 #include "particle.H"
 #include "KinematicParcel.H"
 #include "MPPICParcel.H"
 #include "ThermoParcel.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 namespace Foam
 {
    typedef MPPICParcel < ThermoParcel < KinematicParcel < particle > > > basicThermoKinematicMPPICParcel;
    template<>
    inline bool contiguous<basicThermoKinematicMPPICParcel>()
    {
        return true;
    }
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #endif
 // ************************************************************************* //
--- a/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/defineBasicThermoKinematicMPPICParcel.C
+++ b/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/defineBasicThermoKinematicMPPICParcel.C
@ -1,38 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013 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 "basicThermoKinematicMPPICParcel.H"
 #include "Cloud.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 namespace Foam
 {
    defineTemplateTypeNameAndDebug(basicThermoKinematicMPPICParcel, 0);
    defineTemplateTypeNameAndDebug(Cloud<basicThermoKinematicMPPICParcel>, 0);
 }
 // ************************************************************************* //
--- a/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/makeBasicThermoKinematicMPPICParcelSubmodels.C
+++ b/src/lagrangian/intermediate/parcels/derived/basicThermoKinematicMPPICParcel/makeBasicThermoKinematicMPPICParcelSubmodels.C
@ -1,67 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2013-2015 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 "basicThermoKinematicMPPICCloud.H"
 #include "makeParcelCloudFunctionObjects.H"
 // Kinematic sub-models
 #include "makeParcelForces.H"
 #include "makeParcelDispersionModels.H"
 #include "makeParcelInjectionModels.H"
 #include "makeParcelPatchInteractionModels.H"
 #include "makeParcelStochasticCollisionModels.H"
 #include "makeParcelSurfaceFilmModels.H"
 // MPPIC sub-models
 #include "makeMPPICParcelDampingModels.H"
 #include "makeMPPICParcelIsotropyModels.H"
 #include "makeMPPICParcelPackingModels.H"
 // Thermo sub-models
 #include "makeParcelHeatTransferModels.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 makeParcelCloudFunctionObjects(basicThermoKinematicMPPICCloud);
 // Kinematic sub-models
 makeParcelForces(basicThermoKinematicMPPICCloud);
 makeParcelDispersionModels(basicThermoKinematicMPPICCloud);
 makeParcelInjectionModels(basicThermoKinematicMPPICCloud);
 makeParcelPatchInteractionModels(basicThermoKinematicMPPICCloud);
 makeParcelStochasticCollisionModels(basicThermoKinematicMPPICCloud);
 makeParcelSurfaceFilmModels(basicThermoKinematicMPPICCloud);
 // MPPIC sub-models
 makeMPPICParcelDampingModels(basicThermoKinematicMPPICCloud);
 makeMPPICParcelIsotropyModels(basicThermoKinematicMPPICCloud);
 makeMPPICParcelPackingModels(basicThermoKinematicMPPICCloud);
 // Thermo sub-models
 makeParcelHeatTransferModels(basicThermoKinematicMPPICCloud);
 // ************************************************************************* //
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C
@ -86,17 +86,4 @@ Foam::basicChemistryModel::~basicChemistryModel()
 {}
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 void Foam::basicChemistryModel::edm
 (
    const scalarField& mixing,
    const scalar A,
    const scalar B
 )
 {
    NotImplemented;
 }
 // ************************************************************************* //
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
@ -159,14 +159,6 @@ public:
                //- Calculates the reaction rates
                virtual void calculate() = 0;
                //- Calculates the eddy dissipation reaction rates
                virtual void edm
                (
                    const scalarField& mixing,
                    const scalar A,
                    const scalar B
                );
                //- Solve the reaction system for the given time step
                //  and return the characteristic time
                virtual scalar solve(const scalar deltaT) = 0;
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
@ -756,104 +756,6 @@ void Foam::chemistryModel<CompType, ThermoType>::calculate()
 }
 template<class CompType, class ThermoType>
 void Foam::chemistryModel<CompType, ThermoType>::edm
 (
    const scalarField& mixing,
    const scalar A,
    const scalar B
 )
 {
    if (!this->chemistry_)
    {
        return;
    }
    const volScalarField rho
    (
        IOobject
        (
            "rho",
            this->time().timeName(),
            this->mesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        this->thermo().rho()
    );
    forAll(rho, celli)
    {
        const scalar rhoi = rho[celli];
        const scalar mixingi = mixing[celli];
        scalarField Yi(nSpecie_, 0.0);
        for (label k=0; k<nSpecie_; k++)
        {
            Yi[k] = Y_[k][celli];
        }
        scalarField om(nSpecie_, 0.0);
        scalar minYR = 1.0;
        scalar YP = 0.0;
        forAll(reactions_, m)
        {
            const Reaction<ThermoType>& R = reactions_[m];
            scalar num = 0.0;
            scalar den = 0.0;
            forAll(R.lhs(), s)
            {
                const label k = R.lhs()[s].index;
                const scalar vk = R.lhs()[s].stoichCoeff + SMALL;
                const scalar Wk = specieThermo_[k].W();
                minYR = min(minYR, Yi[k]/Wk/vk);
            }
            forAll(R.rhs(), s)
            {
                const label k = R.rhs()[s].index;
                const scalar vk = R.rhs()[s].stoichCoeff;
                const scalar Wk = specieThermo_[k].W();
                num += Yi[k];
                den += vk*Wk;
            }
            YP = num/den;
            scalar omegam = min (A*rhoi*mixingi*minYR, A*B*rhoi*mixingi*YP);
            forAll(R.lhs(), s)
            {
                const label k = R.lhs()[s].index;
                const scalar vk = R.lhs()[s].stoichCoeff;
                const scalar Wk = specieThermo_[k].W();
                om[k] -= vk*Wk*omegam;
            }
            forAll(R.rhs(), s)
            {
                const label k = R.rhs()[s].index;
                const scalar vk = R.rhs()[s].stoichCoeff;
                const scalar Wk = specieThermo_[k].W();
                om[k] += vk*Wk*omegam;
            }
        }
        for (label k=0; k<nSpecie_; k++)
        {
            RR_[k][celli] = om[k];
        }
    }
 }
 template<class CompType, class ThermoType>
 template<class DeltaTType>
 Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
@ -191,14 +191,6 @@ public:
        //- Calculates the reaction rates
        virtual void calculate();
        //- Calculates the eddy dissipation reaction rates
        virtual void edm
        (
            const scalarField& mixing,
            const scalar A,
            const scalar B
        );
        // Chemistry model functions (overriding abstract functions in
        // basicChemistryModel.H)
--- a/src/thermophysicalModels/radiation/Make/files
+++ b/src/thermophysicalModels/radiation/Make/files
@ -25,7 +25,6 @@ submodels/absorptionEmissionModel/binaryAbsorptionEmission/binaryAbsorptionEmiss
 submodels/absorptionEmissionModel/greyMeanAbsorptionEmission/greyMeanAbsorptionEmission.C
 submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C
 submodels/absorptionEmissionModel/greyMeanSolidAbsorptionEmission/greyMeanSolidAbsorptionEmission.C
 submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.C
 /* Soot model */
 submodels/sootModel/sootModel/sootModel.C
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.C
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.C
@ -1,206 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 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 "wsggmAbsorptionEmission.H"
 #include "addToRunTimeSelectionTable.H"
 #include "unitConversion.H"
 #include "zeroGradientFvPatchFields.H"
 #include "basicSpecieMixture.H"
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 namespace Foam
 {
    namespace radiation
    {
        defineTypeNameAndDebug(wsggmAbsorptionEmission, 0);
        addToRunTimeSelectionTable
        (
            absorptionEmissionModel,
            wsggmAbsorptionEmission,
            dictionary
        );
    }
 }
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 Foam::radiation::wsggmAbsorptionEmission::wsggmAbsorptionEmission
 (
    const dictionary& dict,
    const fvMesh& mesh
 )
 :
    absorptionEmissionModel(dict, mesh),
    coeffsDict_((dict.subDict(typeName + "Coeffs"))),
    thermo_(mesh.lookupObject<fluidThermo>("thermophysicalProperties")),
    emissivityCoeffs_(coeffsDict_.lookup("emissivityCoeffs")),
    fittingFactors_(coeffsDict_.lookup("fittingFactors")),
    pathLength_(coeffsDict_.lookup("pathLength"))
 {
    if (!isA<basicMultiComponentMixture>(thermo_))
    {
        FatalErrorIn
        (
            "radiation::wsggmAbsorptionEmission::wsggmAbsorptionEmission"
            "("
                "const dictionary&, "
                "const fvMesh&"
            ")"
        )   << "Model requires a multi-component thermo package"
            << abort(FatalError);
    }
 }
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 Foam::radiation::wsggmAbsorptionEmission::~wsggmAbsorptionEmission()
 {}
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 Foam::tmp<Foam::volScalarField>
 Foam::radiation::wsggmAbsorptionEmission::aCont(const label bandI) const
 {
    const basicSpecieMixture& mixture =
        dynamic_cast<const basicSpecieMixture&>(thermo_);
    const volScalarField& T = thermo_.T();
    const volScalarField& p = thermo_.p();
    tmp<volScalarField> ta
    (
        new volScalarField
        (
            IOobject
            (
                "aCont" + name(bandI),
                mesh().time().timeName(),
                mesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh(),
            dimensionedScalar("a", dimless/dimLength, 0.0),
            zeroGradientFvPatchVectorField::typeName
        )
    );
    scalarField& a = ta.ref().primitiveFieldRef();
    label indexCO2, indexH2O;
    scalar wCO2, wH2O;
    forAll(mixture.Y(), specieI)
    {
        if(mixture.Y(specieI).name() == "CO2")
        {
            indexCO2 = specieI;
            wCO2 = mixture.W(indexCO2);
        }
        else if(mixture.Y(specieI).name() == "H2O")
        {
            indexH2O = specieI;
            wH2O = mixture.W(indexH2O);
        }
    }
    forAll(a,celli)
    {
        scalar invWt = 0.0;
        forAll(mixture.Y(), specieI)
        {
            invWt += mixture.Y(specieI)[celli]/mixture.W(specieI);
        }
        scalar meanWt = 1.0/invWt; //Mean molecular weight [kg/kmole]
        // Sum of partial pressures of all absorbing gases [atm]
        scalar Pcell_atm =  paToAtm(meanWt*p[celli]*(mixture.Y(indexCO2)[celli]/wCO2 + mixture.Y(indexH2O)[celli]/wH2O));
        // P*S at Eq.(1) of Ref. [atm*m]
        scalar presPathLength = Pcell_atm*pathLength_.value();
        // Limit cell temperature [K]
        scalar Tcell = min(T[celli],2400.0);
        scalar emissivity = 0.0;
        forAll(emissivityCoeffs_, i) //k_i at Eq.(1) of Ref.
        {
            scalar weightingFactor = 0.0;
            for(label j=0 ; j<4 ; j++)
            {
                //a_epsilon_i at Eq.(3) of Ref.
                weightingFactor += fittingFactors_[i][j]*pow(Tcell, j);
            }
            //epsilon at Eq.(1) of Ref.
            emissivity += weightingFactor * (1.0 - exp(-1.0*emissivityCoeffs_[i] * presPathLength));
        }
        emissivity = min(emissivity, 0.99999);
        //Fluent theory guide(ver. 13) Eq.(5-106) [1/m]
        a[celli] = max(1e-2, (-1.0)*log(1.0-emissivity)/pathLength_.value());
    }
    ta.ref().correctBoundaryConditions();
    return ta;
 }
 Foam::tmp<Foam::volScalarField>
 Foam::radiation::wsggmAbsorptionEmission::eCont(const label bandI) const
 {
   return aCont(bandI);
 }
 Foam::tmp<Foam::volScalarField>
 Foam::radiation::wsggmAbsorptionEmission::ECont(const label bandI) const
 {
    tmp<volScalarField> E
    (
        new volScalarField
        (
            IOobject
            (
                "ECont" + name(bandI),
                mesh_.time().timeName(),
                mesh_,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh_,
            dimensionedScalar("E", dimMass/dimLength/pow3(dimTime), 0.0)
        )
    );
    return E;
 }
 // ************************************************************************* //
--- a/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.H
+++ b/src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wsggmAbsorptionEmission/wsggmAbsorptionEmission.H
@ -1,163 +0,0 @@
 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2012 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::radiation::wsggmAbsorptionEmission
 Description
    wsggmAbsorptionEmission radiation absorption and emission coefficients
    for continuous phase
    Ref. T.F Smith, Z.F. Shen and J.N. Friedman
         Evaluation of Coefficients for the Weighted Sum of Gray Gases Model,
         Transactions of the ASME, Vol. 104, pp. 602-608
    Valid when Ptotal = 1atm, 0.001 < P*S < 10.0atm.m, 600 < T < 2400K
               (where P*S = sum of partial pressures of absorbing gases
                            * mean beam length)
    i.e. dictionary
    \verbatim
        wsggmAbsorptionEmissionCoeffs
        {
            pathLength pathLength [ 0 1  0  0  0  0  0] 0.251; //3.6V/S
            // Pw/Pc = 2.0
            emissivityCoeffs 3(0.4201 6.516 131.9);
            fittingFactors
            3
            (
                4(6.508e-1 -5.551e-4 3.029e-7  -5.353e-11)
                4(-0.2504e-1  6.112e-4  -3.882e-7  6.528e-11)
                4(2.718e-1 -3.118e-4  1.221e-7  -1.612e-11)
            );
        }
    \endverbatim
 SourceFiles
    wsggmAbsorptionEmission.C
 \*---------------------------------------------------------------------------*/
 #ifndef wsggmAbsorptionEmission_H
 #define wsggmAbsorptionEmission_H
 #include "interpolationLookUpTable.H"
 #include "absorptionEmissionModel.H"
 #include "HashTable.H"
 #include "absorptionCoeffs.H"
 #include "fluidThermo.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 namespace Foam
 {
 namespace radiation
 {
 /*---------------------------------------------------------------------------*\
                 Class wsggmAbsorptionEmission Declaration
 \*---------------------------------------------------------------------------*/
 class wsggmAbsorptionEmission
 :
    public absorptionEmissionModel
 {
 public:
 private:
    // Private data
    //- Absorption model dictionary
    dictionary coeffsDict_;
    //- SLG thermo package
    const fluidThermo& thermo_;
    //- absorption coefficitent for individual grey gases
    scalarList emissivityCoeffs_;
    //- fitting Factors
    scalarListList fittingFactors_;
    //characteristic path length [m] from dict
    dimensionedScalar pathLength_;
 public:
    //- Runtime type information
    TypeName("wsggmAbsorptionEmission");
    // Constructors
        //- Construct from components
        wsggmAbsorptionEmission(const dictionary& dict, const fvMesh& mesh);
    //- Destructor
    virtual ~wsggmAbsorptionEmission();
    // Member Functions
        // Access
            // Absorption coefficient
                //- Absorption coefficient for continuous phase
                tmp<volScalarField> aCont(const label bandI = 0) const;
            // Emission coefficient
                //- Emission coefficient for continuous phase
                tmp<volScalarField> eCont(const label bandI = 0) const;
            // Emission contribution
                //- Emission contribution for continuous phase
                tmp<volScalarField> ECont(const label bandI = 0) const;
    // Member Functions
        inline bool isGrey() const
        {
            return true;
        }
 };
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 } // End namespace radiation
 } // End namespace Foam
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #endif
 // ************************************************************************* //
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L
@ -148,7 +148,7 @@ irreversibleReactionDelimiter {space}"=>"{space}
 startPDependentSpecie {space}"("{space}"+"{space}
 pDependentSpecie      {specieName}")"{space}
 reactionCoeffs        {space}{floatNum}{some_space}{floatNum}{some_space}{floatNum}{space}
-reactionKeyword       {space}[A-Za-z](([A-Za-z0-9)*-])|("("[^+]))*{space}
+reactionKeyword       {space}[A-Za-z](([A-Za-z0-9)*-])|("("[^+]))*((\+({some_space}))?){space}
 reactionKeywordSlash  {reactionKeyword}"/"{space}
 thirdBodyEfficiency   {space}{floatNum}{space}"/"{space}
 startReactionCoeffs   {space}"/"{space}
@ -637,7 +637,8 @@ bool finishReaction = false;
                    (
                        currentSpecieName,
                        1.0,
-                        molecularWeight(currentSpecieComposition)
+                        molecularWeight(currentSpecieComposition),
                        chargeNumber(currentSpecieComposition)
                    ),
                    currentLowT,
                    currentHighT,
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C
@ -144,6 +144,28 @@ Foam::scalar Foam::chemkinReader::molecularWeight
 }
 Foam::scalar Foam::chemkinReader::chargeNumber
 (
    const List<specieElement>& specieComposition
 ) const
 {
    scalar nElemCharges = 0.0;
    forAll(specieComposition, i)
    {
        label nAtoms = specieComposition[i].nAtoms;
        const word& elementName = specieComposition[i].elementName;
        if (elementName[0] == 'e')
        {
            nElemCharges -= nAtoms;
        }
    }
    return nElemCharges;
 }
 void Foam::chemkinReader::checkCoeffs
 (
    const scalarList& reactionCoeffs,
--- a/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.H
+++ b/src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.H
@ -254,6 +254,11 @@ private:
            const List<specieElement>& specieComposition
        ) const;
        scalar chargeNumber
        (
            const List<specieElement>& specieComposition
        ) const;
        void finishElements(labelList& currentAtoms);
        void checkCoeffs
--- a/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C
+++ b/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C
@ -68,6 +68,70 @@ Foam::scalar Foam::SpecieMixture<MixtureType>::W
 }
 template<class MixtureType>
 Foam::scalar Foam::SpecieMixture<MixtureType>::z
 (
    const label speciei
 ) const
 {
    return this->getLocalThermo(speciei).z();
 }
 template<class MixtureType>
 Foam::scalar Foam::SpecieMixture<MixtureType>::Qc
 (
    const label speciei
 ) const
 {
    return this->getLocalThermo(speciei).Qc();
 }
 template<class MixtureType>
 Foam::scalar Foam::SpecieMixture<MixtureType>::Qc2
 (
    const label speciei
 ) const
 {
    return this->getLocalThermo(speciei).Qc2();
 }
 template<class MixtureType>
 void Foam::SpecieMixture<MixtureType>::Cp
 (
    scalarField& Cps,
    const scalar p,
    const scalar T
 ) const
 {
    forAll(this->species(), speciei)
    {
        Cps[speciei] = this->getLocalThermo(speciei).Cp(p, T);
    }
    return;
 }
 template<class MixtureType>
 void Foam::SpecieMixture<MixtureType>::Cv
 (
    scalarField& Cvs,
    const scalar p,
    const scalar T
 ) const
 {
    forAll(this->species(), speciei)
    {
        Cvs[speciei] = this->getLocalThermo(speciei).Cv(p, T);
    }
    return;
 }
 template<class MixtureType>
 Foam::scalar Foam::SpecieMixture<MixtureType>::Cp
 (
--- a/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H
+++ b/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H
@ -81,6 +81,34 @@ public:
            //- Molecular weight of the given specie [kg/kmol]
            virtual scalar W(const label speciei) const;
            //- Number of charges of the given specie []
            virtual scalar z(const label speciei) const;
            //- Specific charge of the given specie []
            virtual scalar Qc(const label specieI) const;
            //- Absolute specific charge of the given specie []
            virtual scalar Qc2(const label specieI) const;
        // All species thermo properties
            //- Heat capacity at constant pressure [J/(kg K)]
            virtual void Cp
            (
                scalarField& Cps,
                const scalar p,
                const scalar T
            ) const;
            //- Heat capacity at constant volume [J/(kg K)]
            virtual void Cv
            (
                scalarField& Cvs,
                const scalar p,
                const scalar T
            ) const;
        // Per specie thermo properties
--- a/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.C
+++ b/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.C
@ -81,4 +81,72 @@ Foam::tmp<Foam::volScalarField> Foam::basicSpecieMixture::W() const
 }
 Foam::tmp<Foam::volScalarField> Foam::basicSpecieMixture::Qc() const
 {
    const PtrList<volScalarField>& Y(basicMultiComponentMixture::Y());
    tmp<volScalarField> tQc
    (
        new volScalarField
        (
            IOobject
            (
                IOobject::groupName("Qc", Y[0].group()),
                Y[0].time().timeName(),
                Y[0].mesh()
            ),
            Y[0].mesh(),
            dimensionedScalar("zero", dimCurrent*dimTime/dimMass, 0)
        )
    );
    volScalarField& Qc = tQc.ref();
    dimensionedScalar Qci("Qci", dimCurrent*dimTime/dimMass, 0);
    forAll(Y, i)
    {
        Qci.value() = this->Qc(i);
        Qc += Y[i]*Qci;
        // Qc += Y[i]*this->Qc(i);
    }
    return tQc;
 }
 Foam::tmp<Foam::volScalarField> Foam::basicSpecieMixture::Qc2() const
 {
    const PtrList<volScalarField>& Y(basicMultiComponentMixture::Y());
    tmp<volScalarField> tQc
    (
        new volScalarField
        (
            IOobject
            (
                IOobject::groupName("Qc2", Y[0].group()),
                Y[0].time().timeName(),
                Y[0].mesh()
            ),
            Y[0].mesh(),
            dimensionedScalar("zero", sqr(dimCurrent*dimTime)/dimMass, 0)
        )
    );
    volScalarField& Qc = tQc.ref();
    dimensionedScalar Qc2i("Qc2i", sqr(dimCurrent*dimTime)/dimMass, 0);
    forAll(Y, i)
    {
        Qc2i.value() = this->Qc2(i);
        Qc += Y[i]*Qc2i;
        //Qc += Y[i]*this->Qc2(i);
    }
    return tQc;
 }
 // ************************************************************************* //
--- a/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.H
+++ b/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.H
@ -88,9 +88,43 @@ public:
            //- Molecular weight of the given specie [kg/kmol]
            virtual scalar W(const label speciei) const = 0;
            //- Number of charges of the given specie []
            virtual scalar z(const label speciei) const = 0;
            //- Specific charge of the given specie [C/kg]
            virtual scalar Qc(const label speciei) const = 0;
            //- Absolute specific charge of the given specie [C/kg]
            virtual scalar Qc2(const label speciei) const = 0;
            //- Molecular weight of the mixture [kg/kmol]
            tmp<volScalarField> W() const;
            //- Specific charge of the mixture [C/kg]
            tmp<volScalarField> Qc() const;
            //- Absolute specific charge of the mixture [C^2/kg]
            tmp<volScalarField> Qc2() const;
        // All species thermo properties
            //- Heat capacity at constant pressure [J/(kg K)]
            virtual void Cp
            (
                scalarField& Cps,
                const scalar p,
                const scalar T
            ) const = 0;
            //- Heat capacity at constant volume [J/(kg K)]
            virtual void Cv
            (
                scalarField& Cvs,
                const scalar p,
                const scalar T
            ) const = 0;
        // Per specie thermo properties
--- a/src/thermophysicalModels/specie/specie/specie.C
+++ b/src/thermophysicalModels/specie/specie/specie.C
@ -29,6 +29,24 @@ License
 /* * * * * * * * * * * * * * * public constants  * * * * * * * * * * * * * * */
 //- Universal gas constant (default in [J/(kmol K)])
 const Foam::scalar Foam::specie::RR = constant::physicoChemical::R.value()*1000;
 //- Standard pressure (default in [Pa])
 const Foam::scalar Foam::specie::Pstd = constant::standard::Pstd.value();
 //- Standard temperature (default in [K])
 const Foam::scalar Foam::specie::Tstd = constant::standard::Tstd.value();
 //- Elementary charge (default in [C])
 const Foam::scalar Foam::specie::e = constant::electromagnetic::e.value();
 //- Avogadro number (default in [1/mol])
 const Foam::scalar Foam::specie::NA = constant::physicoChemical::NA.value()*1000;
 //- Boltzmann constant (default in [J/K])
 const Foam::scalar Foam::specie::k = constant::physicoChemical::k.value();
 namespace Foam
 {
    defineTypeNameAndDebug(specie, 0);
@ -41,7 +59,8 @@ Foam::specie::specie(Istream& is)
 :
    name_(is),
    nMoles_(readScalar(is)),
-    molWeight_(readScalar(is))
+    molWeight_(readScalar(is)),
    nCharges_(readScalar(is))
 {
    is.check("specie::specie(Istream& is)");
 }
@ -51,7 +70,8 @@ Foam::specie::specie(const dictionary& dict)
 :
    name_(dict.dictName()),
    nMoles_(readScalar(dict.subDict("specie").lookup("nMoles"))),
-    molWeight_(readScalar(dict.subDict("specie").lookup("molWeight")))
+    molWeight_(readScalar(dict.subDict("specie").lookup("molWeight"))),
    nCharges_(dict.subDict("specie").lookupOrDefault("nCharges", 0.0))
 {}
@ -62,6 +82,7 @@ void Foam::specie::write(Ostream& os) const
    dictionary dict("specie");
    dict.add("nMoles", nMoles_);
    dict.add("molWeight", molWeight_);
    dict.add("nCharges", nCharges_);
    os  << indent << dict.dictName() << dict;
 }
@ -72,7 +93,8 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const specie& st)
 {
    os  << st.name_ << tab
        << st.nMoles_ << tab
-        << st.molWeight_;
+        << st.molWeight_ << tab
        << st.nCharges_;
    os.check("Ostream& operator<<(Ostream& os, const specie& st)");
    return os;
--- a/src/thermophysicalModels/specie/specie/specie.H
+++ b/src/thermophysicalModels/specie/specie/specie.H
@ -77,6 +77,9 @@ class specie
        //- Molecular weight of specie [kg/kmol]
        scalar molWeight_;
        //- Number of elementary charges of specie
        scalar nCharges_;
 public:
@ -84,6 +87,29 @@ public:
    ClassName("specie");
    // Public constants
        // Thermodynamic constants
            //- Universal gas constant [J/(kmol K)]
            static const scalar RR;
            //- Standard pressure [Pa]
            static const scalar Pstd;
            //- Standard temperature [K]
            static const scalar Tstd;
            //- Elementary charge [C]
            static const scalar e;
            //- Avogadro number [1/kmol]
            static const scalar NA;
            //- Boltzmann constant [J/K]
            static const scalar k;
    // Constructors
@ -98,6 +124,23 @@ public:
            const scalar molWeight
        );
        //- Construct from components without name
        inline specie
        (
            const scalar nMoles,
            const scalar molWeight,
            const scalar nCharges
        );
        //- Construct from components with name
        inline specie
        (
            const word& name,
            const scalar nMoles,
            const scalar molWeight,
            const scalar nCharges
        );
        //- Construct as copy
        inline specie(const specie&);
@ -127,6 +170,27 @@ public:
            //- Gas constant [J/(kg K)]
            inline scalar R() const;
            //- Charge number
            inline scalar z() const;
        // Electric charge function
            // Mole specific properties
                //- Electric charge [C/kmol]
                inline scalar qc() const;
            // Mass specific properties
                //- Electric charge [C/kg]
                inline scalar Qc() const;
                //- Electric charge squared [C^2/kg]
                inline scalar Qc2() const;
        // I-O
--- a/src/thermophysicalModels/specie/specie/specieI.H
+++ b/src/thermophysicalModels/specie/specie/specieI.H
@ -41,7 +41,8 @@ inline specie::specie
 :
    name_(name),
    nMoles_(nMoles),
-    molWeight_(molWeight)
+    molWeight_(molWeight),
    nCharges_(0.0)
 {}
@ -52,7 +53,36 @@ inline specie::specie
 )
 :
    nMoles_(nMoles),
-    molWeight_(molWeight)
+    molWeight_(molWeight),
    nCharges_(0.0)
 {}
 inline specie::specie
 (
    const word& name,
    const scalar nMoles,
    const scalar molWeight,
    const scalar nCharges
 )
 :
    name_(name),
    nMoles_(nMoles),
    molWeight_(molWeight),
    nCharges_(nCharges)
 {}
 inline specie::specie
 (
    const scalar nMoles,
    const scalar molWeight,
    const scalar nCharges
 )
 :
    nMoles_(nMoles),
    molWeight_(molWeight),
    nCharges_(nCharges)
 {}
@ -62,7 +92,8 @@ inline specie::specie(const specie& st)
 :
    name_(st.name_),
    nMoles_(st.nMoles_),
-    molWeight_(st.molWeight_)
+    molWeight_(st.molWeight_),
    nCharges_(st.nCharges_)
 {}
@ -70,7 +101,8 @@ inline specie::specie(const word& name, const specie& st)
 :
    name_(name),
    nMoles_(st.nMoles_),
-    molWeight_(st.molWeight_)
+    molWeight_(st.molWeight_),
    nCharges_(st.nCharges_)
 {}
@ -100,6 +132,30 @@ inline scalar specie::R() const
 }
 inline scalar specie::z() const
 {
    return nCharges_;
 }
 inline scalar specie::qc() const
 {
    return nCharges_ * NA * e;
 }
 inline scalar specie::Qc() const
 {
    return nCharges_ * NA * e / molWeight_;
 }
 inline scalar specie::Qc2() const
 {
    return nCharges_ * nCharges_ * NA * e * e / molWeight_;
 }
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 inline void specie::operator=(const specie& st)
@ -107,6 +163,7 @@ inline void specie::operator=(const specie& st)
    //name_ = st.name_;
    nMoles_ = st.nMoles_;
    molWeight_ = st.molWeight_;
    nCharges_ = st.nCharges_;
 }
@ -114,6 +171,10 @@ inline void specie::operator+=(const specie& st)
 {
    scalar sumNmoles = max(nMoles_ + st.nMoles_, SMALL);
    nCharges_ =
        nMoles_/sumNmoles*nCharges_
      + st.nMoles_/sumNmoles*st.nCharges_;
    molWeight_ =
        nMoles_/sumNmoles*molWeight_
      + st.nMoles_/sumNmoles*st.molWeight_;
@ -130,6 +191,10 @@ inline void specie::operator-=(const specie& st)
        diffnMoles = SMALL;
    }
    nCharges_ =
        nMoles_/diffnMoles*nCharges_
      - st.nMoles_/diffnMoles*st.nCharges_;
    molWeight_ =
        nMoles_/diffnMoles*molWeight_
      - st.nMoles_/diffnMoles*st.molWeight_;
@ -154,7 +219,9 @@ inline specie operator+(const specie& st1, const specie& st2)
    (
        sumNmoles,
        st1.nMoles_/sumNmoles*st1.molWeight_
-      + st2.nMoles_/sumNmoles*st2.molWeight_
+      + st2.nMoles_/sumNmoles*st2.molWeight_,
        st1.nMoles_/sumNmoles*st1.nCharges_
      + st2.nMoles_/sumNmoles*st2.nCharges_
    );
 }
@ -171,7 +238,9 @@ inline specie operator-(const specie& st1, const specie& st2)
    (
        diffNmoles,
        st1.nMoles_/diffNmoles*st1.molWeight_
-      - st2.nMoles_/diffNmoles*st2.molWeight_
+      - st2.nMoles_/diffNmoles*st2.molWeight_,
        st1.nMoles_/diffNmoles*st1.nCharges_
      - st2.nMoles_/diffNmoles*st2.nCharges_
    );
 }
@ -181,7 +250,8 @@ inline specie operator*(const scalar s, const specie& st)
    return specie
    (
        s*st.nMoles_,
-        st.molWeight_
+        st.molWeight_,
        st.nCharges_
    );
 }
--- a/src/thermophysicalModels/specie/thermo/thermo/thermo.H
+++ b/src/thermophysicalModels/specie/thermo/thermo/thermo.H
@ -336,6 +336,7 @@ public:
            ) const;
        // I-O
            //- Write to Ostream
--- a/tutorials/combustion/SLFMFoam/FlameD/0/CH4
+++ b/tutorials/combustion/SLFMFoam/FlameD/0/CH4
@ -1,65 +0,0 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
 |  \\    /   O peration     | Version:  4.x                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
    version         2;
    format          ascii;
    class           volScalarField;
    location        "0";
    object          CH4;
 }
 dimensions      [ 0 0 0 0 0 0 0 ];
 internalField   uniform 0;
 boundaryField
 {
    mainjet
    {
        type            fixedValue;
        value           uniform 0.155956;
    }
    pilot
    {
        type            fixedValue;
        value           uniform 0;
    }
    coflow
    {
        type            fixedValue;
        value           uniform 0;
    }
    outlet
    {
        type            zeroGradient;
    }
    burnerwall
    {
        type            zeroGradient;
    }
    slipwall
    {
        type            slip;
    }
    axis
    {
        type            symmetryPlane;
    }
    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
 }
 // ************************************************************************* //
--- a/tutorials/combustion/SLFMFoam/FlameD/0/CO
+++ b/tutorials/combustion/SLFMFoam/FlameD/0/CO
@ -1,65 +0,0 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
 |  \\    /   O peration     | Version:  4.x                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
    version         2;
    format          ascii;
    class           volScalarField;
    location        "0/Fluid";
    object          CO;
 }
 dimensions      [ 0 0 0 0 0 0 0 ];
 internalField   uniform 0;
 boundaryField
 {
    axis
    {
        type            symmetryPlane;
    }
    mainjet
    {
        type            fixedValue;
        value           uniform 0;
    }
    pilot
    {
        type            fixedValue;
        value           uniform 0.00407;
    }
    coflow
    {
        type            fixedValue;
        value           uniform 0;
    }
    outlet
    {
        type            zeroGradient;
    }
    burnerwall
    {
        type            zeroGradient;
    }
    slipwall
    {
        type            slip;
    }
    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
 }
 // ************************************************************************* //
--- a/tutorials/combustion/SLFMFoam/FlameD/0/CO2
+++ b/tutorials/combustion/SLFMFoam/FlameD/0/CO2
@ -1,65 +0,0 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
 |  \\    /   O peration     | Version:  4.x                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
    version         2;
    format          ascii;
    class           volScalarField;
    location        "0";
    object          CO2;
 }
 dimensions      [ 0 0 0 0 0 0 0 ];
 internalField   uniform 0;
 boundaryField
 {
    mainjet
    {
        type            fixedValue;
        value           uniform 0;
    }
    pilot
    {
        type            fixedValue;
        value           uniform 0.1098;
    }
    coflow
    {
        type            fixedValue;
        value           uniform 0;
    }
    outlet
    {
        type            zeroGradient;
    }
    burnerwall
    {
        type            zeroGradient;
    }
    slipwall
    {
        type            slip;
    }
    axis
    {
        type            symmetryPlane;
    }
    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
 }
 // ************************************************************************* //
--- a/tutorials/combustion/SLFMFoam/FlameD/0/H
+++ b/tutorials/combustion/SLFMFoam/FlameD/0/H
@ -1,65 +0,0 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
 |  \\    /   O peration     | Version:  4.x                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
    version         2;
    format          ascii;
    class           volScalarField;
    location        "0";
    object          H;
 }
 dimensions      [ 0 0 0 0 0 0 0 ];
 internalField   uniform 0;
 boundaryField
 {
    mainjet
    {
        type            fixedValue;
        value           uniform 0;
    }
    pilot
    {
        type            fixedValue;
        value           uniform 2.48e-05;
    }
    coflow
    {
        type            fixedValue;
        value           uniform 0;
    }
    outlet
    {
        type            zeroGradient;
    }
    burnerwall
    {
        type            zeroGradient;
    }
    slipwall
    {
        type            slip;
    }
    axis
    {
        type            symmetryPlane;
    }
    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
 }
 // ************************************************************************* //
--- a/tutorials/combustion/SLFMFoam/FlameD/0/H2
+++ b/tutorials/combustion/SLFMFoam/FlameD/0/H2
@ -1,65 +0,0 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
 |  \\    /   O peration     | Version:  4.x                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
    version         2;
    format          ascii;
    class           volScalarField;
    location        "0";
    object          H2;
 }
 dimensions      [ 0 0 0 0 0 0 0 ];
 internalField   uniform 0;
 boundaryField
 {
    mainjet
    {
        type            fixedValue;
        value           uniform 0;
    }
    pilot
    {
        type            fixedValue;
        value           uniform 0.000129;
    }
    coflow
    {
        type            fixedValue;
        value           uniform 0;
    }
    outlet
    {
        type            zeroGradient;
    }
    burnerwall
    {
        type            zeroGradient;
    }
    slipwall
    {
        type            slip;
    }
    axis
    {
        type            symmetryPlane;
    }
    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
 }
 // ************************************************************************* //
--- a/tutorials/combustion/SLFMFoam/FlameD/0/H2O
+++ b/tutorials/combustion/SLFMFoam/FlameD/0/H2O
@ -1,65 +0,0 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
 |  \\    /   O peration     | Version:  4.x                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
    version         2;
    format          ascii;
    class           volScalarField;
    location        "0";
    object          H2O;
 }
 dimensions      [ 0 0 0 0 0 0 0 ];
 internalField   uniform 0;
 boundaryField
 {
    mainjet
    {
        type            fixedValue;
        value           uniform 0;
    }
    pilot
    {
        type            fixedValue;
        value           uniform 0.0942;
    }
    coflow
    {
        type            fixedValue;
        value           uniform 0;
    }
    outlet
    {
        type            zeroGradient;
    }
    burnerwall
    {
        type            zeroGradient;
    }
    slipwall
    {
        type            slip;
    }
    axis
    {
        type            symmetryPlane;
    }
    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
 }
 // ************************************************************************* //
--- a/tutorials/combustion/SLFMFoam/FlameD/0/N2
+++ b/tutorials/combustion/SLFMFoam/FlameD/0/N2
@ -1,65 +0,0 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
 |  \\    /   O peration     | Version:  4.x                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
    version         2;
    format          ascii;
    class           volScalarField;
    location        "0";
    object          N2;
 }
 dimensions      [ 0 0 0 0 0 0 0 ];
 internalField   uniform 0.767;
 boundaryField
 {
    mainjet
    {
        type            fixedValue;
        value           uniform 0.647382;
    }
    pilot
    {
        type            fixedValue;
        value           uniform 0.7342;
    }
    coflow
    {
        type            fixedValue;
        value           uniform 0.767;
    }
    outlet
    {
        type            zeroGradient;
    }
    burnerwall
    {
        type            zeroGradient;
    }
    slipwall
    {
        type            slip;
    }
    axis
    {
        type            symmetryPlane;
    }
    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
 }
 // ************************************************************************* //
--- a/tutorials/combustion/SLFMFoam/FlameD/0/NO
+++ b/tutorials/combustion/SLFMFoam/FlameD/0/NO
@ -1,65 +0,0 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
 |  \\    /   O peration     | Version:  4.x                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
    version         2;
    format          ascii;
    class           volScalarField;
    location        "0";
    object          NO;
 }
 dimensions      [ 0 0 0 0 0 0 0 ];
 internalField   uniform 0;
 boundaryField
 {
    mainjet
    {
        type            fixedValue;
        value           uniform 0;
    }
    pilot
    {
        type            fixedValue;
        value           uniform 4.8e-06;
    }
    coflow
    {
        type            fixedValue;
        value           uniform 0;
    }
    outlet
    {
        type            zeroGradient;
    }
    burnerwall
    {
        type            zeroGradient;
    }
    slipwall
    {
        type            slip;
    }
    axis
    {
        type            symmetryPlane;
    }
    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
 }
 // ************************************************************************* //
--- a/tutorials/combustion/SLFMFoam/FlameD/0/O
+++ b/tutorials/combustion/SLFMFoam/FlameD/0/O
@ -1,65 +0,0 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
 |  \\    /   O peration     | Version:  4.x                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
    version         2;
    format          ascii;
    class           volScalarField;
    location        "0";
    object          O;
 }
 dimensions      [ 0 0 0 0 0 0 0 ];
 internalField   uniform 0;
 boundaryField
 {
    mainjet
    {
        type            fixedValue;
        value           uniform 0;
    }
    pilot
    {
        type            fixedValue;
        value           uniform 0.000747;
    }
    coflow
    {
        type            fixedValue;
        value           uniform 0;
    }
    outlet
    {
        type            zeroGradient;
    }
    burnerwall
    {
        type            zeroGradient;
    }
    slipwall
    {
        type            slip;
    }
    axis
    {
        type            symmetryPlane;
    }
    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
 }
 // ************************************************************************* //
--- a/Show more
+++ b/Show more
Author	SHA1	Message	Date
Yeongdo Park	01b0b1396b	Update .gitlab-ci.yml	2018-10-31 15:46:39 +00:00
Yeongdo Park	9f1a9c60d1	Update .gitlab-ci.yml	2018-10-31 14:27:04 +00:00
Yeongdo Park	75c8f48c51	update CI Lint for new GitLab Runner	2018-10-25 07:31:35 -04:00
ignis	b5bebc52ec	New methods for all species Cp and Cv	2018-07-10 15:24:03 +09:00
ignis	6b3e8c972b	Qc and Qc2 field calculation dimension correctly set	2018-07-10 14:59:50 +09:00
ignis	8670f39b9f	specie and specieMixtures feature for electrical property calculation	2018-04-25 14:56:01 +09:00
ignis	c06d324891	mixture got specific charge field calculation	2018-04-19 13:39:45 +09:00
ignis	3efddbac8a	new attr nCharges_ and related methods of specie and chemkinReader mod for cation name	2018-04-18 11:09:25 +09:00
		`@ -1,2 +0,0 @@`
			`const volScalarField& psi = thermo.psi();`
			`//const volScalarField& T = thermo.T();`
		`@ -1,3 +0,0 @@`
			`ThermoMPPICFoam.C`

			`EXE = $(FOAM_APPBIN)/ThermoMPPICFoam`
		`@ -1,2 +0,0 @@`
			`const volScalarField& T = thermo.T();`
			`const volScalarField& psi = thermo.psi();`
		`@ -1,3 +0,0 @@`
			`simpleCoalcombustionFoam.C`

			`EXE = $(FOAM_APPBIN)/simpleCoalcombustionFoam`