/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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/>.

\*---------------------------------------------------------------------------*/

#include "FSD.H"
#include "addToRunTimeSelectionTable.H"
#include "LESModel.H"
#include "fvcGrad.H"
#include "fvcDiv.H"

namespace Foam
{
namespace combustionModels
{

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class CombThermoType, class ThermoType>
FSD<CombThermoType, ThermoType>::FSD
(
    const word& modelType,
    const fvMesh& mesh,
    const word& phaseName
)
:
    singleStepCombustion<CombThermoType, ThermoType>
    (
        modelType,
        mesh,
        phaseName
    ),
    reactionRateFlameArea_
    (
        reactionRateFlameArea::New
        (
            this->coeffs(),
            this->mesh(),
            *this
        )
    ),
    ft_
    (
        IOobject
        (
            IOobject::groupName("ft", phaseName),
            this->mesh().time().timeName(),
            this->mesh(),
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        this->mesh(),
        dimensionedScalar("zero", dimless, 0.0)
    ),
    YFuelFuelStream_(dimensionedScalar("YFuelStream", dimless, 1.0)),
    YO2OxiStream_(dimensionedScalar("YOxiStream", dimless, 0.23)),
    Cv_(readScalar(this->coeffs().lookup("Cv"))),
    C_(5.0),
    ftMin_(0.0),
    ftMax_(1.0),
    ftDim_(300),
    ftVarMin_(readScalar(this->coeffs().lookup("ftVarMin")))
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

template<class CombThermoType, class ThermoType>
FSD<CombThermoType, ThermoType>::~FSD()
{}


// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //

template<class CombThermoType, class ThermoType>
void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
{
    this->singleMixturePtr_->fresCorrect();

    const label fuelI = this->singleMixturePtr_->fuelIndex();

    const volScalarField& YFuel = this->thermoPtr_->composition().Y()[fuelI];

    const volScalarField& YO2 = this->thermoPtr_->composition().Y("O2");

    const dimensionedScalar s = this->singleMixturePtr_->s();

    ft_ =
        (s*YFuel - (YO2 - YO2OxiStream_))/(s*YFuelFuelStream_ + YO2OxiStream_);


    volVectorField nft(fvc::grad(ft_));

    volScalarField mgft(mag(nft));

    surfaceVectorField SfHat(this->mesh().Sf()/this->mesh().magSf());

    volScalarField cAux(scalar(1) - ft_);

    dimensionedScalar dMgft = 1.0e-3*
        (ft_*cAux*mgft)().weightedAverage(this->mesh().V())
       /((ft_*cAux)().weightedAverage(this->mesh().V()) + SMALL)
      + dimensionedScalar("ddMgft", mgft.dimensions(), SMALL);

    mgft += dMgft;

    nft /= mgft;

    const volVectorField& U = YO2.db().lookupObject<volVectorField>("U");

    const volScalarField sigma
    (
        (nft & nft)*fvc::div(U) - (nft & fvc::grad(U) & nft)
    );

    reactionRateFlameArea_->correct(sigma);

    const volScalarField& omegaFuel = reactionRateFlameArea_->omega();


    const scalar ftStoich =
        YO2OxiStream_.value()
       /(
            s.value()*YFuelFuelStream_.value() + YO2OxiStream_.value()
        );

    tmp<volScalarField> tPc
    (
        new volScalarField
        (
            IOobject
            (
                IOobject::groupName("Pc", this->phaseName_),
                U.time().timeName(),
                U.db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            U.mesh(),
            dimensionedScalar("Pc", dimless, 0)
        )
    );

    volScalarField& pc = tPc.ref();

    tmp<volScalarField> tomegaFuel
    (
        new volScalarField
        (
            IOobject
            (
                IOobject::groupName("omegaFuelBar", this->phaseName_),
                U.time().timeName(),
                U.db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            U.mesh(),
            dimensionedScalar
            (
                "omegaFuelBar",
                omegaFuel.dimensions(),
                0
            )
        )
    );

    volScalarField& omegaFuelBar = tomegaFuel.ref();

    // Calculation of the mixture fraction variance (ftVar)
    const compressible::LESModel& lesModel =
        YO2.db().lookupObject<compressible::LESModel>
        (
            turbulenceModel::propertiesName
        );

    const volScalarField& delta = lesModel.delta();
    const volScalarField ftVar(Cv_*sqr(delta)*sqr(mgft));

    // Thickened flame (average flame thickness for counterflow configuration
    // is 1.5 mm)

    volScalarField  deltaF
    (
        lesModel.delta()/dimensionedScalar("flame", dimLength, 1.5e-3)
    );

    // Linear correlation between delta and flame thickness
    volScalarField omegaF(max(deltaF*(4.0/3.0) + (2.0/3.0), scalar(1)));

    scalar deltaFt = 1.0/ftDim_;

    forAll(ft_, celli)
    {
        if (ft_[celli] > ftMin_ && ft_[celli] < ftMax_)
        {
            scalar ftCell = ft_[celli];

            if (ftVar[celli] > ftVarMin_) //sub-grid beta pdf of ft_
            {
                scalar ftVarc = ftVar[celli];
                scalar a =
                    max(ftCell*(ftCell*(1.0 - ftCell)/ftVarc - 1.0), 0.0);
                scalar b = max(a/ftCell - a, 0.0);

                for (int i=1; i<ftDim_; i++)
                {
                    scalar ft = i*deltaFt;
                    pc[celli] += pow(ft, a-1.0)*pow(1.0 - ft, b - 1.0)*deltaFt;
                }

                for (int i=1; i<ftDim_; i++)
                {
                    scalar ft = i*deltaFt;
                    omegaFuelBar[celli] +=
                        omegaFuel[celli]/omegaF[celli]
                       *exp
                        (
                           -sqr(ft - ftStoich)
                           /(2.0*sqr(0.01*omegaF[celli]))
                        )
                       *pow(ft, a - 1.0)
                       *pow(1.0 - ft, b - 1.0)
                       *deltaFt;
                }
                omegaFuelBar[celli] /= max(pc[celli], 1e-4);
            }
            else
            {
                omegaFuelBar[celli] =
                   omegaFuel[celli]/omegaF[celli]
                  *exp(-sqr(ftCell - ftStoich)/(2.0*sqr(0.01*omegaF[celli])));
            }
        }
        else
        {
            omegaFuelBar[celli] = 0.0;
        }
    }


    // Combustion progress variable, c

    List<label> productsIndex(2, label(-1));
    {
        label i = 0;
        forAll(this->singleMixturePtr_->specieProd(), specieI)
        {
            if (this->singleMixturePtr_->specieProd()[specieI] < 0)
            {
                productsIndex[i] = specieI;
                i++;
            }
        }
    }


    // Flamelet probability of the progress c based on IFC (reuse pc)
    scalar YprodTotal = 0;
    forAll(productsIndex, j)
    {
        YprodTotal += this->singleMixturePtr_->Yprod0()[productsIndex[j]];
    }

    forAll(ft_, celli)
    {
        if (ft_[celli] < ftStoich)
        {
            pc[celli] = ft_[celli]*(YprodTotal/ftStoich);
        }
        else
        {
            pc[celli] = (1.0 - ft_[celli])*(YprodTotal/(1.0 - ftStoich));
        }
    }

    tmp<volScalarField> tproducts
    (
        new volScalarField
        (
            IOobject
            (
                IOobject::groupName("products", this->phaseName_),
                U.time().timeName(),
                U.db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            U.mesh(),
            dimensionedScalar("products", dimless, 0)
        )
    );

    volScalarField& products = tproducts.ref();

    forAll(productsIndex, j)
    {
        label specieI = productsIndex[j];
        const volScalarField& Yp = this->thermoPtr_->composition().Y()[specieI];
        products += Yp;
    }

    volScalarField c
    (
        max(scalar(1) - products/max(pc, scalar(1e-5)), scalar(0))
    );

    pc = min(C_*c, scalar(1));

    const volScalarField fres(this->singleMixturePtr_->fres(fuelI));

    this->wFuel_ == mgft*pc*omegaFuelBar;
}


template<class CombThermoType, class ThermoType>
void FSD<CombThermoType, ThermoType>::correct()
{
    this->wFuel_ ==
        dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0);

    if (this->active())
    {
        calculateSourceNorm();
    }
}


template<class CombThermoType, class ThermoType>
bool FSD<CombThermoType, ThermoType>::read()
{
    if (singleStepCombustion<CombThermoType, ThermoType>::read())
    {
        this->coeffs().lookup("Cv") >> Cv_ ;
        this->coeffs().lookup("ftVarMin") >> ftVarMin_;
        reactionRateFlameArea_->read(this->coeffs());
        return true;
    }
    else
    {
        return false;
    }
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace combustionModels
} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //