#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(), "ha", "ea");

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 "createRadiationModel.H"

#include "createMRF.H"