/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
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
    rhoCentralFoam

Group
    grpCompressibleSolvers

Description
    Density-based compressible flow solver based on
    central-upwind schemes of Kurganov and Tadmor with
    support for mesh-motion and topology changes.

    Chemistry reaction solver for multi-gas.

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

#include "gpufvCFD.H"
#include "dynamicFvMesh.H"
#include "psiReactionThermo.H"
#include "reactingMixture.H"
#include "thermoPhysicsTypes.H"
#include "CombustionModel.H"
#include "gpumultivariateScheme.H"
#include "turbulentFluidThermoModel.H"
#include "fixedRhoFvPatchScalargpuField.H"
#include "directionInterpolate.H"
#include "gpulocalEulerDdtScheme.H"
#include "gpufvcSmooth.H"
#include <sys/time.h>

struct my_timer
{
    struct timeval start_time, end_time;
    double time_use;

    void start()
    {
        gettimeofday(&start_time, NULL);
    }

    void stop()
    {
        gettimeofday(&end_time, NULL);
        time_use = (end_time.tv_sec - start_time.tv_sec) + (double)(end_time.tv_usec - start_time.tv_usec)/1000000.0;
    }
};

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
template<class ThermoType>
struct compHsFunctor
{
    ThermoType *speciesData;
    const label id;
    compHsFunctor(ThermoType* _speciesData, const label _id):speciesData(_speciesData),id(_id){}
    __host__ __device__
    scalar operator()(const thrust::tuple<scalar,scalar>& t)
    {
        const scalar p = thrust::get<0>(t);
        const scalar T = thrust::get<1>(t);
        return speciesData[id].Hs(p,T);
    }

};
}


int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Density-based compressible flow solver based on"
        " central-upwind schemes of Kurganov and Tadmor with"
        " support for mesh-motion and topology changes."
        "Chemistry reaction solver for multi-gas."
    );

    #define NO_CONTROL
    #include "gpupostProcess.H"

    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createDynamicFvMesh.H"
    #include "gpucreateMesh.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
    #include "createTimeControls.H"
    #include "thermodata.H"
    double totaldensitytime=0.0,totalmomentumtime=0.0,totalspeciestime=0.0,totalenergytime=0.0,totalODEtime=0.0,totalturbulencetime=0.0;    
    turbulence->validate();

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

    #include "readFluxScheme.H"

    const dimensionedScalar v_zero(dimVolume/dimTime, Zero);

    // Courant numbers used to adjust the time-step
    scalar CoNum = 0.0;
    scalar meanCoNum = 0.0;

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.run())
    {
        #include "readTimeControls.H"

        if (!LTS)
        {
            #include "setDeltaT.H"

            ++runTime;

            // Do any mesh changes
            mesh.update();
        }

        // --- Directed interpolation of primitive fields onto faces

        surfaceScalargpuField rho_pos(interpolate(rho, pos));
        surfaceScalargpuField rho_neg(interpolate(rho, neg));

        surfaceVectorgpuField rhoU_pos(interpolate(rhoU, pos, U.name()));
        surfaceVectorgpuField rhoU_neg(interpolate(rhoU, neg, U.name()));

        volScalargpuField rPsi("rPsi", 1.0/psi);
        surfaceScalargpuField rPsi_pos(interpolate(rPsi, pos, T.name()));
        surfaceScalargpuField rPsi_neg(interpolate(rPsi, neg, T.name()));

        surfaceScalargpuField e_pos(interpolate(e, pos, T.name()));
        surfaceScalargpuField e_neg(interpolate(e, neg, T.name()));

        surfaceVectorgpuField U_pos("U_pos", rhoU_pos/rho_pos);
        surfaceVectorgpuField U_neg("U_neg", rhoU_neg/rho_neg);

        surfaceScalargpuField p_pos("p_pos", rho_pos*rPsi_pos);
        surfaceScalargpuField p_neg("p_neg", rho_neg*rPsi_neg);

	    PtrList<surfaceScalargpuField> Y_pos(Y.size()), Y_neg(Y.size()); 
        forAll(Y,i)
	    {
	        Y_pos.set(i,new surfaceScalargpuField (fvc::interpolate(Y[i], pos, "reconstruct(Y)")));
	        Y_neg.set(i,new surfaceScalargpuField (fvc::interpolate(Y[i], neg, "reconstruct(Y)")));
	    }

        surfaceScalargpuField phiv_pos("phiv_pos", U_pos & devicemesh.Sf());
        // Note: extracted out the orientation so becomes unoriented
        phiv_pos.setOriented(false);
        surfaceScalargpuField phiv_neg("phiv_neg", U_neg & devicemesh.Sf());
        phiv_neg.setOriented(false);

        // Make fluxes relative to mesh-motion
        if (mesh.moving())
        {
            surfaceScalargpuField meshPhi(devicemesh.phi());
            meshPhi.setOriented(false);
            phiv_pos -= meshPhi;
            phiv_neg -= meshPhi;
        }

        volScalargpuField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
        surfaceScalargpuField cSf_pos
        (
            "cSf_pos",
            interpolate(c, pos, T.name())*devicemesh.magSf()
        );

        surfaceScalargpuField cSf_neg
        (
            "cSf_neg",
            interpolate(c, neg, T.name())*devicemesh.magSf()
        );

        surfaceScalargpuField ap
        (
            "ap",
            max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
        );

        surfaceScalargpuField am
        (
            "am",
            min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
        );

        surfaceScalargpuField a_pos("a_pos", ap/(ap - am));

        surfaceScalargpuField amaxSf("amaxSf", max(mag(am), mag(ap)));

        surfaceScalargpuField aSf("aSf", am*a_pos);

        if (fluxScheme == "Tadmor")
        {
            aSf = -0.5*amaxSf;
            a_pos = 0.5;
        }

        surfaceScalargpuField a_neg("a_neg", 1.0 - a_pos);

        phiv_pos *= a_pos;
        phiv_neg *= a_neg;

        surfaceScalargpuField aphiv_pos("aphiv_pos", phiv_pos - aSf);
        surfaceScalargpuField aphiv_neg("aphiv_neg", phiv_neg + aSf);

        // Reuse amaxSf for the maximum positive and negative fluxes
        // estimated by the central scheme
        amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));

        #include "centralCourantNo.H"

        if (LTS)
        {
            #include "setRDeltaT.H"

            ++runTime;
        }

        Info<< "Time = " << runTime.timeName() << nl << endl;

        phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;

        surfaceVectorgpuField phiU(aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg);
        // Note: reassembled orientation from the pos and neg parts so becomes
        // oriented
        phiU.setOriented(true);

        surfaceVectorgpuField phiUp(phiU + (a_pos*p_pos + a_neg*p_neg)*devicemesh.Sf());

        surfaceScalargpuField phiEp
        (
            "phiEp",
            aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
          + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
          + aSf*p_pos - aSf*p_neg
        );

        // Make flux for pressure-work absolute
        if (mesh.moving())
        {
            surfaceScalargpuField meshPhi(devicemesh.phi());
            meshPhi.setOriented(false);
            phiEp += meshPhi*(a_pos*p_pos + a_neg*p_neg);
        }

        volScalargpuField muEff("muEff", turbulence->muEff());
        volTensorgpuField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));
        my_timer tm_density;tm_density.start();
        // --- Solve density
        solve(fvm::ddt(rho) + fvc::div(phi));Info<<"------Timecost density = "<<tm_density.time_use<<"s-------------------"<< nl << endl;totaldensitytime+=tm_density.time_use;
        rho.max(0.0001);	
		rho.min(50.0); 
        my_timer tm_momentum;tm_momentum.start();
        // --- Solve momentum
        solve(fvm::ddt(rhoU) + fvc::div(phiUp));Info<<"------Timecost momentum = "<<tm_momentum.time_use<<"s-------------------"<< nl << endl;totalmomentumtime+=tm_momentum.time_use;

        U.ref() =
            rhoU()
           /rho();
        U.correctBoundaryConditions();
        rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField();

        if (!inviscid)
        {
            solve
            (
                fvm::ddt(rho, U) - fvc::ddt(rho, U)
              - fvm::laplacian(muEff, U)
              - fvc::div(tauMC)
            );
            rhoU = rho*U;
        }

		//for(label ic=0; ic<3; ic++)
        //{	
	    //   volScalarField Ui=U.component(ic);
	    //   Ui.max(dimensionedScalar("U_limit", U.dimensions(), -2000.0));
	    //   Ui.min(dimensionedScalar("U_limit", U.dimensions(), 2000.0));
	    //   U.replace(ic,Ui);
        //   Info<< "min/max(U[" << ic<< "]) = " << min(Ui).value() << ", " << max(Ui).value() << endl;
        //}

		Info<< "min/max(U) = " << min(mag(U)).value() << ", " << max(mag(U)).value() << endl;
		my_timer tm_species;tm_species.start();
        // --- Solve species equation        
		PtrList<surfaceScalargpuField> phiY(Y.size()); 
	    forAll(Y, i) 
	    { 
	        phiY.set
	       (
		      i,
		      new surfaceScalargpuField
		      ( 
		         (aphiv_pos*(rho_pos*Y_pos[i]) + aphiv_neg*(rho_neg*Y_neg[i]))	 
		      )	
	        ); 
	    }
		
	    Yt=0.0;  
        Y[inertIndex]=1.0;
	    forAll(Y, i)
	    {
			if (Y[i].name() != inertSpecie)
		    {   
				solve(fvm::ddt(rhoY[i]) + fvc::div(phiY[i]));
				Y[i] = rhoY[i]/rho;  
				Y[i].correctBoundaryConditions();
				Yt += Y[i];
				Y[i].min(1.0);
		    }
		}
		Y[inertIndex] -=Yt;
		Y[inertIndex].correctBoundaryConditions();
	    Yt.correctBoundaryConditions();
        Info<< "min/max(Yt after convection) = " << min(Yt).value() << ", " << max(Yt).value() << endl;
        
		forAll(Y, i)
	    {
            Y[i].correctBoundaryConditions();
	        rhoY[i].boundaryFieldRef() == rho.boundaryField()*Y[i].boundaryField();
	    }
	    if(!inviscid)
	    {
	        forAll(Y,i)
			{
			DiffsY[i]=turbulence->mut()/Sct+turbulence->mu()/Sc[i];
			}
			Y[inertIndex]=1.0;
	        forAll(Y, i)
	        {
		       if (Y[i].name() != inertSpecie)
		       {   
		          volScalargpuField& Yi = Y[i];
		          gpufvScalarMatrix YiEqn
		          (
			         fvm::ddt(rho, Yi)- fvc::ddt(rho, Yi) 
		             == fvm::laplacian(DiffsY[i], Yi)
		          );
		          YiEqn.relax();   
				  YiEqn.solve(mesh.solver("Yi"));
				  //Yi.max(0.0); 
                  Yi.min(1.0);    
				  Y[inertIndex] -=Y[i];
		       }
	        }
            Yt=0.0;    
			Yt.correctBoundaryConditions();
	        forAll(Y,i)	    
            {   
                Yt +=Y[i];
                Y[i].correctBoundaryConditions();  
                rhoY[i] = rho*Y[i];
            }
            Yt.correctBoundaryConditions();
	    }
        Info<< "min/max(Yt after diffusion) = " << min(Yt).value() << ", " << max(Yt).value() << endl;tm_species.stop();Info<<"------Timecost species = "<<tm_species.time_use<<"s-------------------"<< nl << endl;totalspeciestime+=tm_species.time_use;
        my_timer tm_energy;tm_energy.start();
        // --- Solve energy
        surfaceScalargpuField sigmaDotU
        (
            "sigmaDotU",
            (
                fvc::interpolate(muEff)*devicemesh.magSf()*fvc::snGrad(U)
              + fvc::dotInterpolate(devicemesh.Sf(), tauMC)
            )
          & (a_pos*U_pos + a_neg*U_neg)
        );

        solve
        (
            fvm::ddt(rhoE)
          + fvc::div(phiEp)
          - fvc::div(sigmaDotU)
        );

        e = rhoE/rho - 0.5*magSqr(U);
        e.correctBoundaryConditions();
        thermo.correct();
        rhoE.boundaryFieldRef() ==
            rho.boundaryField()*
            (
                e.boundaryField() + 0.5*magSqr(U.boundaryField())
            );

        if (!inviscid)
        {
			//sensible enthalpy diffusion term induced by species diffusion
			forAll(HsY,i)
			{
                thrust::transform(
                        thrust::make_zip_iterator(thrust::make_tuple(p.begin(),T.begin())),
                        thrust::make_zip_iterator(thrust::make_tuple(p.end(),T.end())),
                        HsY[i].begin(),
                        compHsFunctor<gasEThermoPhysics>(raw_pointer_cast(gSpeciesData),i)
                        ); 
                        
				HsY[i].correctBoundaryConditions();
			}
			volVectorgpuField HsDiff("HsDiff",DiffsY[0]*HsY[0]*fvc::grad(Y[0]));
			
			for(label i=1;i<Y.size();i++)
			{
				HsDiff +=DiffsY[i]*HsY[i]*fvc::grad(Y[i]);	 
			}
			surfaceScalargpuField HsDiffTerm=(devicemesh.Sf()&(fvc::interpolate(HsDiff)))(); 
			
            solve
            (
                fvm::ddt(rho, e) - fvc::ddt(rho, e)
              - fvc::div(HsDiffTerm)
              - fvm::laplacian(turbulence->alphaEff(), e)
            );
            thermo.correct();
            rhoE = rho*(e + 0.5*magSqr(U));
        }
        Info<< "min/max(T) after convection & diffusion = " << min(T).value() << ", " << max(T).value() << endl;
        Info<< "min/max(p) after convection & diffusion = " << min(p).value() << ", " << max(p).value() << endl;tm_energy.stop();Info<<"------Timecost energy = "<<tm_energy.time_use<<"s-------------------"<< nl << endl;totalenergytime+=tm_energy.time_use;

        if(SolveChemistry)
        {my_timer tm_ODE;tm_ODE.start();
            combustion->correct();
			Qdot = combustion->Qdot();	tm_ODE.stop();Info<<"------Timecost ODE = "<<tm_ODE.time_use<<"s-------------------"<< nl << endl;totalODEtime+=tm_ODE.time_use;
			Qdot.correctBoundaryConditions();	
            Info << "Qdot = " << gSum(Qdot.internalField())<<endl;

            Yt = 0.0;   
            forAll(Y,i)
            {
				if(i!=inertIndex && composition.active(i))
				{
					solve(fvm::ddt(rhoY[i])-fvc::ddt(rhoY[i]) == combustion->R2(Y[i]));
                    Y[i]=rhoY[i]/rho;   
					Y[i].max(0.0);
				    Yt += Y[i];
				}
            }
            Yt.correctBoundaryConditions();
			Y[inertIndex] = scalar(1) - Yt;
			Y[inertIndex].max(0.0);
			
			Yt = 0.0;
			forAll(Y, i)
			{
				Yt += Y[i];
			}
            Info<<"min/max(Yt after reaction) = "<<min(Yt).value()<<" , "<<max(Yt).value()<<endl;
            forAll(Y, i)
	        {
            	Y[i].correctBoundaryConditions();
                rhoY[i].boundaryFieldRef() == rho.boundaryField()*Y[i].boundaryField();
            }
	
			solve(fvm::ddt(rhoE)-fvc::ddt(rhoE) == combustion->Qdot());

            e = rhoE/rho - 0.5*magSqr(U); 
            e.correctBoundaryConditions(); 
            thermo.correct();
            rhoE.boundaryFieldRef() == 
				rho.boundaryField()*
				(
					e.boundaryField() + 0.5*magSqr(U.boundaryField())
				);
        }
	
        p.ref() =
            rho()
           /psi();
        p.correctBoundaryConditions();
        rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();

		if(SolveChemistry)
		{
           Info<< "min/max(T) after combustion = " << min(T).value() << ", " << max(T).value() << endl;
           Info<< "min/max(p) after combustion = " << min(p).value() << ", " << max(p).value() << endl; 
		}
        my_timer tm_turbulence;tm_turbulence.start();
        turbulence->correct();
        tm_turbulence.stop();Info<<"------Timecost turbulence = "<<tm_turbulence.time_use<<"s-------------------"<< nl << endl;totalturbulencetime+=tm_turbulence.time_use;
        runTime.write();

        runTime.printExecutionTime(Info);
    }

    Info<< "End\n" << endl;
Info<<"------TotalTimecost density = "<<totaldensitytime<<"s-------------------"<< nl << endl;
Info<<"------TotalTimecost momentum = "<<totalmomentumtime<<"s-------------------"<< nl << endl;
Info<<"------TotalTimecost species = "<<totalspeciestime<<"s-------------------"<< nl << endl;
Info<<"------TotalTimecost energy = "<<totalenergytime<<"s-------------------"<< nl << endl;
Info<<"------TotalTimecost ODE = "<<totalODEtime<<"s-------------------"<< nl << endl;
Info<<"------TotalTimecost turbulence = "<<totalturbulencetime<<"s-------------------"<< nl << endl;
    return 0;
}

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