/*---------------------------------------------------------------------------*\
  ==  == ====== ====   ====    |
                    \\     ||  | 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/>.

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

#include "cipsaMultiphaseEuler.H"
#include "constrainHbyA.H"
#include "constrainPressure.H"
#include "findRefCell.H"
#include "fvcDiv.H"
#include "fvcSnGrad.H"
#include "fvmLaplacian.H"
#include "fvcFlux.H"
#include "fvcReconstruct.H"

// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //

void Foam::solvers::cipsaMultiphaseEuler::pressureCorrector()
{
    volScalarField& p_rgh = p_rgh_;

    rho_ = fluid.rho();

    // Correct p_rgh for consistency with p and the updated densities due to
    // updated temperature
    p_rgh = p - rho*buoyancy.gh - buoyancy.pRef;

    // Construct pressure "diffusivity"
    surfaceScalarField rAUf
    (
        IOobject
        (
            "rAUf",
            runTime.name(),
            mesh
        ),
        mesh,
        dimensionedScalar(dimTime/dimDensity, Zero)
    );

    // rAs, rAUfs, alpharAUfs, rAUf WITH field relaxation of UAlpha
    rAUfs.setSize(phases.size());
    rAs.clear();
    if (fluid.implicitPhasePressure())
    {
        rAs.setSize(phases.size());
    }
    PtrList<surfaceScalarField> alpharAUfs(phases.size());
    forAll(movingPhases, movingPhasei)
    {
        const phaseModel& phase = movingPhases[movingPhasei];

        if (fluid.implicitPhasePressure())
        {
            rAs.set
            (
                phase.index(),
                1.0
               /(
                    HUEqns[phase.index()].A()
                   /max(relaxU[phase.index()], small)
                  + ddtUEqns[phase.index()].A()
                  + dmdtUEqns[phase.index()].A()
                  + fvModelsUEqns[phase.index()].A()
                  + dragDiag[phase.index()]
                  + vmDiag[phase.index()]
                )
            );
        }

        rAUfs.set
        (
            phase.index(),
            1.0
           /(
                1.0/max(relaxU[phase.index()], small)
              + rHUfs[phase.index()]
               *(
                    fvc::interpolate(ddtUEqns[phase.index()].A())
                  + fvc::interpolate(dmdtUEqns[phase.index()].A())
                  + fvc::interpolate(fvModelsUEqns[phase.index()].A())
                  + dragDiagf[phase.index()]
                  + vmDiagf[phase.index()]
                )
            )
        );

        // Phase diagonal coefficients
        alpharAUfs.set
        (
            phase.index(),
            rHUfs[phase.index()]*rAUfs[phase.index()]
           *(alphafs[phase.index()] + dragPCoefff[phase.index()])
        );

        // Pressure "diffusivity"
        rAUf += alphafs[phase.index()]*alpharAUfs[phase.index()];
    }

    rAUf = mag(rAUf);

    // Mass transfer rates
    PtrList<volScalarField::Internal> dmdts(populationBalanceSystem_.dmdts());

    // --- Pressure corrector loop
    while (pimple.correct())
    {
        // Correct fixed-flux BCs to be consistent with the velocity BCs
        fluid_.correctBoundaryFlux();

        PtrList<surfaceScalarField> phigs(phases.size());
        {
            // Mixture density gradient times gravitational acceleration
            const surfaceScalarField ghSnGradRho
            (
                "ghSnGradRho",
                buoyancy.ghf*fvc::snGrad(rho)*mesh.magSf()
            );

            // phigs WITH field relaxation of UAlpha
            forAll(movingPhases, movingPhasei)
            {
                const phaseModel& phase = movingPhases[movingPhasei];

                // Combined buoyancy and force fluxes
                phigs.set
                (
                    phase.index(),
                    alpharAUfs[phase.index()]
                   *(
                      - ghSnGradRho
                      + fluid.resolvedInterface().surfaceTension(phase)
                       *mesh.magSf()
                    )
                  - rHUfs[phase.index()]*rAUfs[phase.index()]
                   *phase.turbulentNormalStress()*mesh.magSf()
                );
            }
        }

        // Update explicit drag terms for PISO loop
        PtrList<volVectorField> dragSource(calcDragSource());
        PtrList<surfaceScalarField> dragSourcef(calcDragSourcef());

        // Explicit part of the time derivative of the virtual mass for
        // PISO loop
        PtrList<volVectorField> vmSource(calcVmSource());
        PtrList<surfaceScalarField> vmSourcef(calcVmSourcef());

        // Update lift, wall lubrication and turbulent dispersion forces for
        // PISO loop
        PtrList<surfaceScalarField> Ffs(momentumTransferSystem_.Ffs());

        // Buoyancy for dispersed phases and parallel to underresolved
        // interfaces
        PtrList<volScalarField> rhoAlphaBuoyancy(calcRhoAlphaBuoyancy());
        PtrList<volVectorField> buoyancySource
        (
            calcBuoyancySource(rhoAlphaBuoyancy)
        );
        PtrList<surfaceScalarField> buoyancySourcef
        (
            calcBuoyancySourcef(rhoAlphaBuoyancy)
        );
        rhoAlphaBuoyancy.clear();

        // Construct Total predicted flux
        surfaceScalarField phiHbyA
        (
            IOobject
            (
                "phiHbyA",
                runTime.name(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh,
            dimensionedScalar(dimVolumetricFlux, Zero)
        );

        PtrList<volVectorField> HbyAs(phases.size());
        PtrList<surfaceScalarField> phiHbyAs(phases.size());

        forAll(movingPhases, movingPhasei)
        {
            const phaseModel& phase = movingPhases[movingPhasei];
            const volVectorField& UAlpha = phase.URef();

            // Pseudo convected velocities
            HbyAs.set
            (
                phase.index(),
                rHUs[phase.index()]*HUEqns[phase.index()].H()
            );
            constrainHbyA(HbyAs[phase.index()], UAlpha, p_rgh);

            // Predicted fluxes for each phase
            phiHbyAs.set
            (
                phase.index(),
                constrainPhiHbyA
                (
                    rAUfs[phase.index()]
                   *(
                        fvc::flux(HbyAs[phase.index()])
                      + (1 - relaxU[phase.index()])
                       /max(relaxU[phase.index()], small)
                       *MRF.absolute
                        (
                            phases[phase.index()].phi()().prevIter()
                        )
                      + rHUfs[phase.index()]
                       *(
                            ddtPhiCoeffs[phase.index()]
                          + fvc::flux
                            (
                                dmdtUEqns[phase.index()].H()
                              + fvModelsUEqns[phase.index()].H()
                            )
                          - Ffs[phase.index()]
                          + buoyancySourcef[phase.index()]
                        )
                      + dragSourcef[phase.index()]
                      + vmSourcef[phase.index()]
                    ),
                    UAlpha,
                    p_rgh
                )
              + phigs[phase.index()]
            );

            phiHbyA += alphafs[phase.index()]*phiHbyAs[phase.index()];
        }

        MRF.makeRelative(phiHbyA);

        // Update the fixedFluxPressure BCs to ensure flux consistency
        {
            surfaceScalarField::Boundary phib
            (
                surfaceScalarField::Internal::null(),
                phi.boundaryField()
            );
            phib = 0;

            forAll(movingPhases, movingPhasei)
            {
                const phaseModel& phase = movingPhases_[movingPhasei];

                phib +=
                    alphafs[phase.index()].boundaryField()
                   *phase.phi()().boundaryField();
            }

            setSnGrad<fixedFluxPressureFvPatchScalarField>
            (
                p_rgh.boundaryFieldRef(),
                (
                    phiHbyA.boundaryField() - phib
                )
               /(mesh.magSf().boundaryField()*rAUf.boundaryField())
            );
        }

        // Compressible pressure equations
        PtrList<fvScalarMatrix> pEqnComps(compressibilityEqns(dmdts));

        // Cache p prior to solve for density update
        volScalarField p_rgh_0(p_rgh);

        // Iterate over the pressure equation to correct for
        // non-orthogonality
        while (pimple.correctNonOrthogonal())
        {
            // Construct the transport part of the pressure equation
            fvScalarMatrix pEqnIncomp
            (
                fvc::div(phiHbyA)
              - fvm::laplacian(rAUf, p_rgh)
            );

            // Solve
            {
                fvScalarMatrix pEqn(pEqnIncomp);

                forAll(phases, phasei)
                {
                    pEqn += pEqnComps[phasei];
                }

                if (fluid.incompressible())
                {
                        pEqn.setReference
                        (
                            pressureReference.refCell(),
                            pressureReference.refValue()
                        );
                }

                pEqn.solve();
            }

            // Correct fluxes and velocities on last non-orthogonal
            // iteration
            if (pimple.finalNonOrthogonalIter())
            {
                phi_ = phiHbyA + pEqnIncomp.flux();

                surfaceScalarField mSfGradp
                (
                    "mSfGradp",
                    pEqnIncomp.flux()/rAUf
                );

                forAll(movingPhases, movingPhasei)
                {
                    phaseModel& phase = movingPhases_[movingPhasei];

                    phase.phiRef() =
                        phiHbyAs[phase.index()]
                      + alpharAUfs[phase.index()]*mSfGradp;

                    // Set the phase dilatation rate
                    phase.divU(-pEqnComps[phase.index()] & p_rgh);
                }

                // Optionally relax pressure for velocity correction
                p_rgh.relax();

                mSfGradp = pEqnIncomp.flux()/rAUf;

                // Momentum corrector
                forAll(movingPhases, movingPhasei)
                {
                    phaseModel& phase = movingPhases_[movingPhasei];
                    volVectorField& UAlpha = phase.URef();

                    if (faceMomentum)
                    {
                        UAlpha = fvc::reconstruct(phase.phi());
                    }
                    else
                    {
                        UAlpha =
                            1.0
                           /(
                                HUEqns[phase.index()].A()
                               /max(relaxU[phase.index()], small)
                              + ddtUEqns[phase.index()].A()
                              + dmdtUEqns[phase.index()].A()
                              + fvModelsUEqns[phase.index()].A()
                              + dragDiag[phase.index()]
                              + vmDiag[phase.index()]
                            )
                           *(
                                HUEqns[phase.index()].A()
                               *HbyAs[phase.index()]
                              + (1 - relaxU[phase.index()])
                               /max(relaxU[phase.index()], small)
                               *HUEqns[phase.index()].A()*UAlpha.prevIter()
                              + ddtUEqns[phase.index()].H()
                              + dmdtUEqns[phase.index()].H()
                              + fvModelsUEqns[phase.index()].H()
                              + dragSource[phase.index()]
                              + vmSource[phase.index()]
                              + buoyancySource[phase.index()]
                            )
                            + fvc::reconstruct
                            (
                                phigs[phase.index()]
                              + alpharAUfs[phase.index()]*mSfGradp
                              - rHUfs[phase.index()]*rAUfs[phase.index()]
                               *Ffs[phase.index()]
                            );
                    }

                    UAlpha.correctBoundaryConditions();
                    fvConstraints().constrain(UAlpha);

                    MRF.makeRelative(phase.phiRef());
                }
            }
        }

        // Update the mixture properties and static pressure
        U_ = fluid.U();
        p_ = p_rgh + rho*buoyancy.gh + buoyancy.pRef;
        fvConstraints().constrain(p_);

        // Account for static pressure reference
        if (p_rgh.needReference() && fluid.incompressible())
        {
            p_ += dimensionedScalar
            (
                "p",
                p.dimensions(),
                pressureReference.refValue()
              - getRefCellValue(p, pressureReference.refCell())
            );
        }

        // Limit p_rgh
        p_rgh = p - rho*buoyancy.gh - buoyancy.pRef;

        // Update densities from change in p_rgh
        forAll(phases, phasei)
        {
            phaseModel& phase = phases_[phasei];
            if (!phase.incompressible())
            {
                phase.rho() +=
                    phase.fluidThermo().psi()*(p_rgh - p_rgh_0);
            }
        }

        // Correct p_rgh for consistency with p and the updated densities
        rho_ = fluid.rho();
        p_rgh = p - rho*buoyancy.gh - buoyancy.pRef;
    }

    alphafs.clear();

    HUEqns.clear();
    relaxU.clear();
    ddtUEqns.clear();
    ddtPhiCoeffs.clear();
    dmdtUEqns.clear();
    fvModelsUEqns.clear();

    rHUs.clear();
    rHUfs.clear();
    rADs.clear();
    rADfs.clear();

    dragDiag.clear();
    dragDiagf.clear();
    dragPCoefff.clear();
    vmDiag.clear();
    vmDiagf.clear();

    rAUfs.clear();

    fluid_.correctKinematics();
}


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