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

Class
    Foam::fv::CrankNicolsonDdtScheme

Group
    grpFvDdtSchemes

Description
    Second-oder Crank-Nicolson implicit ddt using the current and
    previous time-step fields as well as the previous time-step ddt.

    The Crank-Nicolson scheme is often unstable for complex flows in complex
    geometries and it is necessary to "off-centre" the scheme to stabilize it
    while retaining greater temporal accuracy than the first-order
    Euler-implicit scheme.  Off-centering is specified via the mandatory
    coefficient \c ocCoeff in the range [0,1] following the scheme name e.g.
    \verbatim
    ddtSchemes
    {
        default         CrankNicolson 0.9;
    }
    \endverbatim
    or with an optional "ramp" function to transition from the Euler scheme to
    Crank-Nicolson over a initial period to avoid start-up problems, e.g.
    \verbatim
    ddtSchemes
    {
        default         CrankNicolson
        ocCoeff
        {
            type scale;
            scale linearRamp;
            duration 0.01;
            value 0.9;
        };
    }
    \endverbatim
    With a coefficient of 1 the scheme is fully centred and second-order,
    with a coefficient of 0 the scheme is equivalent to Euler-implicit.
    A coefficient of 0.9 has been found to be suitable for a range of cases for
    which higher-order accuracy in time is needed and provides similar accuracy
    and stability to the backward scheme.  However, the backward scheme has
    been found to be more robust and provides formal second-order accuracy in
    time.

    The advantage of the Crank-Nicolson scheme over backward is that only the
    new and old-time values are needed, the additional terms relating to the
    fluxes and sources are evaluated at the mid-point of the time-step which
    provides the opportunity to limit the fluxes in such a way as to ensure
    boundedness while maintaining greater accuracy in time compared to the
    Euler-implicit scheme.  This approach is now used with MULES in the
    interFoam family of solvers.  Boundedness cannot be guaranteed with the
    backward scheme.

Note
    The Crank-Nicolson coefficient for the implicit part of the RHS is related
    to the off-centering coefficient by
    \verbatim
        cnCoeff = 1.0/(1.0 + ocCoeff);
    \endverbatim

See also
    Foam::Euler
    Foam::backward

SourceFiles
    CrankNicolsonDdtScheme.C

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

#ifndef gpuCrankNicolsonDdtScheme_H
#define gpuCrankNicolsonDdtScheme_H

#include "gpuddtScheme.H"
#include "Function1.H"

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

namespace Foam
{

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

namespace fv
{

/*---------------------------------------------------------------------------*\
                       Class CrankNicolsonDdtScheme Declaration
\*---------------------------------------------------------------------------*/

template<class Type>
class gpuCrankNicolsonDdtScheme
:
    public fv::gpuddtScheme<Type>
{
    // Private Data

        //- Class to store the ddt0 fields on the objectRegistry for use in the
        //  next time-step.  The start-time index of the CN scheme is also
        //  stored to help handle the transition from Euler to CN
        template<class GeoField>
        class gpuDDt0Field
        :
            public GeoField
        {
            label startTimeIndex_;

        public:

            //- Constructor from file for restart.
            gpuDDt0Field
            (
                const IOobject& io,
                const gpufvMesh& mesh
            );

            //- Constructor from components, initisalised to zero with given
            //  dimensions.
            gpuDDt0Field
            (
                const IOobject& io,
                const gpufvMesh& mesh,
                const dimensioned<typename GeoField::value_type>& dimType
            );

            //- Return the start-time index
            label startTimeIndex() const;

            //- Cast to the underlying GeoField
            GeoField& operator()();

            //- Assignment to a GeoField
            void operator=(const GeoField& gf);
        };


        //- Off-centering coefficient function
        //  1 -> CN, less than one blends with EI
        autoPtr<Function1<scalar>> ocCoeff_;


    // Private Member Functions

        //- No copy construct
        gpuCrankNicolsonDdtScheme(const gpuCrankNicolsonDdtScheme&) = delete;

        //- No copy assignment
        void operator=(const gpuCrankNicolsonDdtScheme&) = delete;

        template<class GeoField>
        gpuDDt0Field<GeoField>& ddt0_
        (
            const word& name,
            const dimensionSet& dims
        );

        //- Check if the ddt0 needs to be evaluated for this time-step
        template<class GeoField>
        bool evaluate(gpuDDt0Field<GeoField>& ddt0);

        //- Return the coefficient for Euler scheme for the first time-step
        //  for and CN thereafter
        template<class GeoField>
        scalar coef_(const gpuDDt0Field<GeoField>&) const;

        //- Return the old time-step coefficient for Euler scheme for the
        //  second time-step and for CN thereafter
        template<class GeoField>
        scalar coef0_(const gpuDDt0Field<GeoField>&) const;

        //- Return the reciprocal time-step coefficient for Euler for the
        //  first time-step and CN thereafter
        template<class GeoField>
        dimensionedScalar rDtCoef_(const gpuDDt0Field<GeoField>&) const;

        //- Return the reciprocal old time-step coefficient for Euler for the
        //  second time-step and CN thereafter
        template<class GeoField>
        dimensionedScalar rDtCoef0_(const gpuDDt0Field<GeoField>&) const;

        //- Return ddt0 multiplied by the off-centreing coefficient
        template<class GeoField>
        tmp<GeoField> offCentre_(const GeoField& ddt0) const;


public:

    //- Runtime type information
    TypeName("CrankNicolson");


    // Constructors

        //- Construct from mesh
        gpuCrankNicolsonDdtScheme(const gpufvMesh& mesh);

        //- Construct from mesh and Istream
        gpuCrankNicolsonDdtScheme(const gpufvMesh& mesh, Istream& is);


    // Member Functions

        //- Return mesh reference
        const gpufvMesh& mesh() const
        {
            return fv::gpuddtScheme<Type>::mesh();
        }

        //- Return the current off-centreing coefficient
        scalar ocCoeff() const
        {
            return ocCoeff_->value(mesh().time().value());
        }

        tmp<GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>> fvcDdt
        (
            const dimensioned<Type>&
        );

        tmp<GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>> fvcDdt
        (
            const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>&
        );

        tmp<GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>> fvcDdt
        (
            const dimensionedScalar&,
            const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>&
        );

        tmp<GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>> fvcDdt
        (
            const volScalargpuField&,
            const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>&
        );

        tmp<GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>> fvcDdt
        (
            const volScalargpuField& alpha,
            const volScalargpuField& rho,
            const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>& psi
        );

        tmp<gpufvMatrix<Type>> fvmDdt
        (
            const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>&
        );

        tmp<gpufvMatrix<Type>> fvmDdt
        (
            const dimensionedScalar&,
            const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>&
        );

        tmp<gpufvMatrix<Type>> fvmDdt
        (
            const volScalargpuField&,
            const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>&
        );

        tmp<gpufvMatrix<Type>> fvmDdt
        (
            const volScalargpuField& alpha,
            const volScalargpuField& rho,
            const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>& psi
        );

        typedef typename gpuddtScheme<Type>::fluxFieldType fluxFieldType;

        tmp<fluxFieldType> fvcDdtUfCorr
        (
            const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>& U,
            const GeometricgpuField<Type, fvsPatchgpuField, gpusurfaceMesh>& Uf
        );

        tmp<fluxFieldType> fvcDdtPhiCorr
        (
            const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>& U,
            const fluxFieldType& phi
        );

        tmp<fluxFieldType> fvcDdtUfCorr
        (
            const volScalargpuField& rho,
            const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>& U,
            const GeometricgpuField<Type, fvsPatchgpuField, gpusurfaceMesh>& Uf
        );

        tmp<fluxFieldType> fvcDdtPhiCorr
        (
            const volScalargpuField& rho,
            const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>& U,
            const fluxFieldType& phi
        );


        tmp<surfaceScalargpuField> meshPhi
        (
            const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>&
        );
};


template<>
tmp<surfaceScalargpuField> gpuCrankNicolsonDdtScheme<scalar>::fvcDdtUfCorr
(
    const GeometricgpuField<scalar, fvPatchgpuField, gpuvolMesh>& U,
    const GeometricgpuField<scalar, fvsPatchgpuField, gpusurfaceMesh>& Uf
);

template<>
tmp<surfaceScalargpuField> gpuCrankNicolsonDdtScheme<scalar>::fvcDdtPhiCorr
(
    const volScalargpuField& U,
    const surfaceScalargpuField& phi
);

template<>
tmp<surfaceScalargpuField> gpuCrankNicolsonDdtScheme<scalar>::fvcDdtUfCorr
(
    const volScalargpuField& rho,
    const volScalargpuField& U,
    const surfaceScalargpuField& Uf
);

template<>
tmp<surfaceScalargpuField> gpuCrankNicolsonDdtScheme<scalar>::fvcDdtPhiCorr
(
    const volScalargpuField& rho,
    const volScalargpuField& U,
    const surfaceScalargpuField& phi
);


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

} // End namespace fv

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

} // End namespace Foam

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

#ifdef NoRepository
    #include "gpuCrankNicolsonDdtScheme.C"
#endif

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

#endif

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