/*---------------------------------------------------------------------------*\
  ==  == ====== ====   ====    |
                    \\     ||  | Multiphase Code Repository by HZDR
  ======   //   ||  || ===//   | Website: https://doi.org/10.14278/rodare.767
  ||  ||  //    ||  // || \\   | License: GPL-3.0-or-later
  ==  == ====== ====   ==  ==  |
-------------------------------------------------------------------------------
License
    This file is part of the Multiphase Code Repository by HZDR.

    Copyright (C) 2025 by Helmholtz-Zentrum Dresden-Rossendorf e.V. (HZDR),
    Website: https://hzdr.de

    Multiphase Code Repository by HZDR is based on the free software for
    computational fluid dynamics (CFD) from the OpenFOAM Foundation.
    Copyright (C) 2025 by OpenFOAM Foundation, Website: https://openfoam.org

    If you are interested in which files are original OpenFOAM Foundation
    files, which OpenFOAM Foundation files were modified, and which files were
    newly created, see FILES.md.

    Multiphase Code Repository by HZDR 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.

    Multiphase Code Repository by HZDR 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 the Multiphase Code Repository by HZDR. If not, see
    <http://www.gnu.org/licenses/>.

Class
    Foam::solvers::cipsaMultiphaseEuler

Description
    Solver module for a system of any number of compressible fluid phases with a
    common pressure, but otherwise separate properties. The type of phase model
    is run time selectable and can optionally represent multiple species and
    in-phase reactions. The phase system is also run time selectable and can
    optionally represent different types of momentum, heat and mass transfer.

    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
    pseudo-transient and steady simulations.

    Optional fvModels and fvConstraints are provided to enhance the simulation
    in many ways including adding various sources, Lagrangian
    particles, surface film etc. and constraining or limiting the solution.

    The solver uses the consistent momentum interpolation method of
    Cubero et al. (2014) and a partial-elimination algorithm.

    Reference:
    \verbatim
        Cubero, A., Sánchez-Insa, A., & Fueyo, N. (2014).
        A consistent momentum interpolation method for steady and unsteady
        multiphase flows.
        Computers & Chemical Engineering, 62, 96-107.
    \endverbatim

SourceFiles
    cipsaMultiphaseEuler.C

See also
    Foam::solvers::multiphaseEuler
    Foam::solvers::compressibleVoF
    Foam::solvers::fluidSolver
    Foam::solvers::incompressibleFluid

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

#ifndef cipsaMultiphaseEuler_H
#define cipsaMultiphaseEuler_H

#include "multiphaseEuler.H"

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

namespace Foam
{
namespace solvers
{

/*---------------------------------------------------------------------------*\
                     Class cipsaMultiphaseEuler Declaration
\*---------------------------------------------------------------------------*/

class cipsaMultiphaseEuler
:
    public multiphaseEuler
{
protected:

    // Controls

        //- Partial elimination drag contribution optimisation
        //  Defaults to true
        Switch partialElimination;

        //- Partial elimination heat transfer contribution optimisation for
        //  energy equation
        //  Defaults to false
        Switch energyPartialElimination;

        //- Add diffusive term to stabilise small phase fraction states
        Switch diffusiveStabilisation;


    // Flow properties

        //- Mixture velocity
        volVectorField U_;


    // Thermophysical properties

        //- Mixture density
        volScalarField rho_;


    // Cached temporary equations

        // Momentum sub-equations

            //- Temporary convective-diffusive phase momentum matrices
            //  shared between the momentum predictor and pressure corrector
            PtrList<fvVectorMatrix> HUEqns;

            //- Temporary temporal phase momentum matrices
            //  shared between the momentum predictor and pressure corrector
            PtrList<fvVectorMatrix> ddtUEqns;

            //- Temporary mass change phase momentum matrices
            //  shared between the momentum predictor and pressure corrector
            PtrList<fvVectorMatrix> dmdtUEqns;

            //- Temporary fvModels phase momentum matrices
            //  shared between the momentum predictor and pressure corrector
            PtrList<fvVectorMatrix> fvModelsUEqns;


    // Cached temporary fields

        // Cell-based coefficients

            //- Reciprocal cell-based convective-diffusive implicit coefficient
            PtrList<volScalarField> rHUs;

            //- Reciprocal cell-based implicit coefficients for partial
            //  elimination drag handling
            PtrList<volScalarField> rADs;

            //- Cell-based implicit coefficient of drag
            PtrList<volScalarField> dragDiag;

            //- Cell-based implicit coefficient of temporal virtual mass
            //  contribution
            PtrList<volScalarField> vmDiag;


        // Face-based coefficients

            //- Face-interpolated phase fractions
            PtrList<surfaceScalarField> alphafs;

            //- Reciprocal face-based convective-diffusive implicit coefficient
            PtrList<surfaceScalarField> rHUfs;

            //- Reciprocal face-based dimensionless momentum equation diagonal
            PtrList<surfaceScalarField> rAUfs;

            //- Reciprocal face-based implicit coefficients for partial
            //  elimination drag handling
            PtrList<surfaceScalarField> rADfs;

            //- Temporary temporal phase coefficients
            //  shared between the momentum predictor and pressure corrector
            PtrList<surfaceScalarField> ddtPhiCoeffs;

            //- Face-based implicit coefficient of drag
            PtrList<surfaceScalarField> dragDiagf;

            //- Face-based pressure coefficient of drag
            PtrList<surfaceScalarField> dragPCoefff;

            //- Face-based implicit coefficient of temporal virtual mass
            //  contribution
            PtrList<surfaceScalarField> vmDiagf;


        //- Factor for consistently relaxing momentum equations
        List<scalar> relaxU;


    // Protected Member Functions

        //- Read controls
        virtual bool read();


private:

    // Field construction

        //- Fill up gaps in a phase-indexed list of fields with zeros
        template
        <
            class Type,
            template<class> class PatchField,
            class GeoMesh
        >
        void fillFields
        (
            const word& name,
            const dimensionSet& dims,
            PtrList<GeometricField<Type, GeoMesh, PatchField>>& fieldList
        ) const;


    // Private Member Functions

        //- Energy and temperature predictor
        virtual void energyPredictor();

        //- Convective diffusive momentum sub-equations
        void updateHUEqn();

        //- Stabilise convective diffusive momentum sub-equations
        void stabiliseHUEqn();

        //- Set factors for consistent relaxation of momentum equations
        void relaxUEqn();

        //- Temporal momentum sub-equations
        void updateDdtUEqn();

        //- Momentum sub-equations containing interfacial momentum transfer due
        //  to mass transfer
        void updateDmdtUEqn();

        //- Momentum sub-equations containing fvModels
        void updateFvModelsUEqn();

        //- Implicit drag coefficients containing partial elimination
        void updateDragDiag();

        //- Cell-based explicit drag coefficients containing partial elimination
        PtrList<volVectorField> calcDragSource();

        //- Face-based explicit drag coefficients containing partial elimination
        PtrList<surfaceScalarField> calcDragSourcef();

        //- Implicit virtual mass ddt coefficients
        void updateVmDiag();

        //- Cell-based explicit virtual mass ddt coefficients
        PtrList<volVectorField> calcVmSource();

        //- Face-based explicit virtual mass ddt coefficients
        PtrList<surfaceScalarField> calcVmSourcef();

        //- Convective virtual mass contribution
        void updateVmConvectionEqn();

        //- Phase-specific densities for phase-specific buoyancy contribution
        PtrList<volScalarField> calcRhoAlphaBuoyancy();

        //- Cell-based phase-specific buoyancy contribution
        PtrList<volVectorField> calcBuoyancySource
        (
            PtrList<volScalarField>& rhoAlphaBuoyancy
        );

        //- Face-based phase-specific buoyancy contribution
        PtrList<surfaceScalarField> calcBuoyancySourcef
        (
            PtrList<volScalarField>& rhoAlphaBuoyancy
        );


public:

    // Public Data

        //- Reference to the mixture velocity
        const volVectorField& U;

        //- Reference to the mixture density
        const volScalarField& rho;


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


    // Constructors

        //- Construct from region mesh
        cipsaMultiphaseEuler(fvMesh& mesh);

        //- Disallow default bitwise copy construction
        cipsaMultiphaseEuler(const cipsaMultiphaseEuler&) = delete;


    //- Destructor
    virtual ~cipsaMultiphaseEuler();


    // Member Functions

        //- Called at the start of the time-step, before the PIMPLE loop
        virtual void preSolve();

        //- Called at the start of the PIMPLE loop to move the mesh
        virtual void moveMesh()
        {};

        //- Called at the start of the PIMPLE loop
        virtual void prePredictor();

        //- Construct and optionally solve the momentum equation
        virtual void momentumPredictor();

        //- Construct and solve the pressure equation in the PISO loop
        virtual void pressureCorrector();

        //- Called after the PIMPLE loop at the end of the time-step
        virtual void postSolve()
        {};


    // Member Operators

        //- Disallow default bitwise assignment
        void operator=(const cipsaMultiphaseEuler&) = delete;
};


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

} // End namespace solvers
} // End namespace Foam

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

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

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

#endif

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