#include "gpucreateRDeltaT.H"

Info<< "Reading thermophysical properties\n" << endl;
autoPtr<psiReactionThermo> pThermo(psiReactionThermo::New(devicemesh));
psiReactionThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");

/*autoPtr<psiThermo> pThermo
(
    psiThermo::New(mesh)
);
psiThermo& thermo = pThermo();*/

	// switch for cold flow or reacting flow
	bool SolveChemistry = readBool(runTime.controlDict().lookup("SolveChemistry"));

basicSpecieMixture& composition = thermo.composition();
PtrList<volScalargpuField>& Y = composition.Y();

const word inertSpecie(thermo.get<word>("inertSpecie"));
if (!composition.species().found(inertSpecie))
{
    FatalIOErrorIn(args.executable().c_str(), thermo)
        << "Inert specie " << inertSpecie << " not found in available species "
        << composition.species() << exit(FatalIOError);
}

volScalargpuField& e = thermo.he();

Info<< "Reading field U\n" << endl;
volVectorgpuField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    devicemesh
);

volScalargpuField rho
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    thermo.rho()
);

volVectorgpuField rhoU
(
    IOobject
    (
        "rhoU",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    rho*U
);

volScalargpuField rhoE
(
    IOobject
    (
        "rhoE",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    rho*(e + 0.5*magSqr(U))
);

surfaceScalargpuField pos
(
    IOobject
    (
        "pos",
        runTime.timeName(),
        mesh
    ),
    devicemesh,
    dimensionedScalar("pos", dimless, 1.0)
);

surfaceScalargpuField neg
(
    IOobject
    (
        "neg",
        runTime.timeName(),
        mesh
    ),
    devicemesh,
    dimensionedScalar("neg", dimless, -1.0)
);

surfaceScalargpuField phi("phi", fvc::flux(rhoU));

Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
    compressible::turbulenceModel::New
    (
        rho,
        U,
        phi,
        thermo
    )
);

Info<< "Creating combustion model\n" << endl;
autoPtr<CombustionModel<psiReactionThermo>> combustion
(
     CombustionModel<psiReactionThermo>::New(thermo, turbulence())
);

gpumultivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;

forAll(Y, i)
{
    fields.add(Y[i]);
}
fields.add(thermo.he());

Info<< "Creating field kinetic energy K\n" << endl;
volScalargpuField K("K", 0.5*magSqr(U));

PtrList<volScalargpuField> rhoY(Y.size()); 
forAll(Y, i)
{
    rhoY.set
    (
        i,
		new volScalargpuField 
		(
	    	IOobject
	    	(
			"rho"+Y[i].name(),
			runTime.timeName(),
			mesh,
			IOobject::NO_READ,
			IOobject::NO_WRITE
	    	),
	    	rho*Y[i]
		)
    );    
}

volScalargpuField Qdot
(
    IOobject
    (
        "Qdot",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    devicemesh,
    dimensionedScalar(dimEnergy/dimVolume/dimTime, 0.0)
);

const label inertIndex(composition.species()[inertSpecie]);

volScalargpuField Yt
(
    IOobject
    (
        "Yt",
         mesh.time().timeName(),
         mesh,
         IOobject::NO_READ,
         IOobject::NO_WRITE
     ),
     devicemesh,
     dimensionedScalar( dimless, 0)//,
);

Info<< "Reading transport properties" <<endl;
IOdictionary transportProperties
(
    IOobject
    (
        "transportProperties",
        runTime.constant(),
        runTime,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
     )
);
dictionary Schmidt(transportProperties.subDict("Sc"));
scalarField Sc(Y.size(),0.75);
forAll(Sc,i)
{
    const word& name = Y[i].name();
    if(Schmidt.found(name))
    {
        Sc[i]=readScalar(Schmidt.lookup(name));
    }
} 
const scalar Sct(Schmidt.lookupOrDefault<scalar>("Sct",0.7));