/*---------------------------------------------------------------------------*\
  ==  == ====== ====   ====    |
                    \\     ||  | 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 "underResolutionModel.H"
#include "phaseSystem.H"
#include "wallDist.H"
#include "fvcSmooth.H"
#include "addonFvcSpread.H"
#include "addToRunTimeSelectionTable.H"
#include "uniformDimensionedFields.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(underResolutionModel, 0);
    defineRunTimeSelectionTable(underResolutionModel, dictionary);
    addNamedToRunTimeSelectionTable
    (
        underResolutionModel,
        underResolutionModel,
        dictionary,
        none
    );
}


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

Foam::tmp<Foam::volScalarField> Foam::underResolutionModel::DeltaI() const
{
    return volScalarField::New
    (
        IOobject::groupName("DeltaI", interface_.name()),
        1.0/max(mag(gradAlphaI()), interface_.fluid().deltaN())
       *pos(interface_.phase1() - interface_.phase1().residualAlpha())
       *pos(interface_.phase2() - interface_.phase2().residualAlpha())
    );
}


Foam::tmp<Foam::volScalarField> Foam::underResolutionModel::Q31
(
    const volScalarField& Lb,
    const volScalarField& gn
) const
{
    return
        sqrt
        (
            0.13*gn*Lb
          + 1.57*interface_.sigma()
           /max
            (
                heavyPhase_.rho()*Lb,
                dimensionedScalar(dimDensity*dimLength, small)
            )
        );
}


Foam::tmp<Foam::volScalarField> Foam::underResolutionModel::Q32
(
    const volScalarField& Lb,
    const volScalarField& gn
) const
{
    return
        sqrt
        (
            7.7e-4*gn*Lb
          + 0.22*interface_.sigma()
           /max
            (
                heavyPhase_.rho()*Lb,
                dimensionedScalar(dimDensity*dimLength, small)
            )
        );
}


Foam::tmp<Foam::volScalarField> Foam::underResolutionModel::blobSize
(
    const volScalarField& q,
    const volScalarField& gn
) const
{
    const volScalarField& epsilon =
        mesh_.lookupObject<volScalarField>
        (
            IOobject::groupName("epsilon", heavyPhase_.name())
        );

    // Turbulence integral length scale (3.4) and (3.6)
    const volScalarField Lb1
    (
        pow3(1.0/sqrt(3.0)*q)
       /max
        (
            epsilon,
            dimensionedScalar(sqr(dimVelocity)/dimTime, small)
        )
    );

    // Mean flow depth via wall distance
    const volScalarField Lb2("Lb2", yWall());

    const uniformDimensionedVectorField& g
    (
        mesh_.lookupObject<uniformDimensionedVectorField>("g")
    );

    // Solve quadratic equation for LbMinus (3.9)
    const dimensionedScalar a
    (
        7.7e-4
       *max(mag(g), dimensionedScalar(dimAcceleration, small))
    );
    const volScalarField b(-sqr(q/(1 + DeltaMu_)));
    const volScalarField c(0.22*interface_.sigma()/heavyPhase_.rho());
    const volScalarField LbMinus
    (
        (
          - b
          - sqrt
            (
                max
                (
                    sqr(b) - 4.0*a*c,
                    dimensionedScalar(pow(dimVelocity, 4), small)
                )
            )
        )
       /(2.0*a)
    );

    // Regime indicator
    const volScalarField FS1(pos(q - Q32(Lb1, gn)));
    const volScalarField FS2(pos(q - Q32(Lb2, gn)));

    // Blob size (3.7) and (3.8)
    return pos0(FS1 - FS2)*Lb1 + neg(FS1 - FS2)*LbMinus;
}


Foam::tmp<Foam::volScalarField> Foam::underResolutionModel::roughness
(
    const volScalarField& yGI,
    const volScalarField& uStarG
) const
{
    // Liquid turbulent velocity
    const volScalarField q(sqrt(2.0*heavyPhase_.k()));

    const uniformDimensionedVectorField& g
    (
        mesh_.lookupObject<uniformDimensionedVectorField>("g")
    );

    // Interface normal gravity contribution (3.3)
    const volScalarField gn
    (
        max(g & interface_.nHat(), dimensionedScalar(dimAcceleration, Zero))
    );

    const volScalarField Lb(blobSize(q, gn));
    const volScalarField q32(Q32(Lb, gn));

    // Turbulence roughness (4.14) and (4.15)
    const volScalarField rt
    (
        IOobject::groupName("turbulenceRoughness", interface_.name()),
        pos(q - q32)*Lb*(q320_ + (q310_ - q320_)*(q - q32)/(Q31(Lb, gn) - q32))
    );

    return
        min
        (
            beta2_*yGI,
            beta1_*sqr(uStarG)
           /max(mag(g), dimensionedScalar(dimAcceleration, small))
          + rt
        );
}


Foam::tmp<Foam::volScalarField> Foam::underResolutionModel::wallLawViscous
(
    const volScalarField& yI,
    const volScalarField& deltaU,
    const volScalarField& nu
) const
{
    return
        sqrt(nu*deltaU/max(yI, dimensionedScalar(dimLength, small)))
       *pos(yI - dimensionedScalar(dimLength, small));
}


Foam::tmp<Foam::volScalarField> Foam::underResolutionModel::wallLawSmooth
(
    const volScalarField& yPlus,
    const volScalarField& yI,
    const volScalarField& deltaU,
    const volScalarField& nu
) const
{
    const volScalarField uPlus(log(max(yPlus, small))/kappa_ + Bs_);

    return
        neg0(yPlus - Cs_)*wallLawViscous(yI, deltaU, nu)
      + pos(yPlus - Cs_)*deltaU/uPlus;
}


Foam::tmp<Foam::volScalarField> Foam::underResolutionModel::wallLawRough
(
    const volScalarField& yPlus,
    const volScalarField& yI,
    const volScalarField& deltaU,
    const volScalarField& nu,
    const volScalarField& Rer
) const
{
    // (4.21) and (4.22)
    const volScalarField xi
    (
        pos(Rer - RerS_)*neg(Rer - RerR_)
       *log(max(Rer/RerS_, small))/log(RerR_/RerS_)
      + pos0(Rer - RerR_)
    );

    // (4.20)
    const volScalarField B
    (
        Bs_
      + (Br_ - Bs_ - log(max(Rer, small))/kappa_)
       *sin(constant::mathematical::piByTwo*xi)
    );

    // (4.19)
    const volScalarField uPlusLM(log(max(yPlus, small)/kappa_) + B);

    const volScalarField condition1
    (
        pos(yPlus - uPlusLM*Rer)*neg(yPlus - Rer)
    );

    // (4.25) and (4.26)
    return deltaU/max(condition1*Br_ + (1 - condition1)*uPlusLM, small);
}


Foam::tmp<Foam::volScalarField> Foam::underResolutionModel::wallLaw
(
    const volScalarField& yPlus,
    const volScalarField& yI,
    const volScalarField& deltaU,
    const volScalarField& nu,
    const volScalarField& Rer
) const
{
    return
        neg(Rer - RerS_)*wallLawSmooth(yPlus, yI, deltaU, nu)
      + pos0(Rer - RerS_)*wallLawRough(yPlus, yI, deltaU, nu, Rer);
}


Foam::tmp<Foam::volScalarField> Foam::underResolutionModel::tau
(
    volScalarField& uStar1,
    volScalarField& uStar2,
    const volScalarField rho1,
    const volScalarField rho2
) const
{
    tmp<volScalarField> tTau
    (
        sqrt(0.5*(sqr(rho1*sqr(uStar1)) + sqr(rho2*sqr(uStar2))))
    );

    uStar1 = sqrt(tTau.ref()/rho1);
    uStar2 = sqrt(tTau.ref()/rho2);

    return tTau;
}


void Foam::underResolutionModel::correctGradAlphaI()
{
    volVectorField gradAlphaI(interface_.gradAlpha());

    volScalarField interfaceIndicator(mag(gradAlphaI));

    smoothen(interfaceIndicator, 4);
    const volScalarField infoLocation
    (
        fvc::spreadInterface(gradAlphaI, interfaceIndicator, nLayers_, small)
    );
    smoothen(gradAlphaI);

    gradAlphaI_ = infoLocation*gradAlphaI;

    if (debug & mesh_.time().writeTime()) gradAlphaI_.write();
}


void Foam::underResolutionModel::correctUIyI()
{
    volScalarField alpha1(interface_.phase1());
    volScalarField alpha2(interface_.phase2());

    smoothen(alpha1);
    smoothen(alpha2);

    // Both individual phases velocities across interface
    U1I_ = interface_.phase1().U();
    U2I_ = interface_.phase2().U();

    const volScalarField infoLocation1
    (
        fvc::spreadInterface(U1I_, alpha1, nLayers_, 0.01)
    );
    const volScalarField infoLocation2
    (
        fvc::spreadInterface(U2I_, alpha2, nLayers_, 0.01)
    );

    U1I_ *= infoLocation1;
    U2I_ *= infoLocation2;

    // Relative velocity across interface
    UrI_ = infoLocation1*infoLocation2*(U2I_ - U1I_);

    // Interface velocity
    const volScalarField mu1(interface_.phase1().fluidThermo().mu());
    const volScalarField mu2(interface_.phase2().fluidThermo().mu());

    UI_ = infoLocation1*infoLocation2*(mu1*U1I_ + mu2*U2I_)/(mu1 + mu2);

    // Distance from interface
    yI_ =
        (alpha1 - dimensionedScalar(dimless, 0.5))
       /max(mag(gradAlphaI()), interface_.fluid().deltaN())
       *pos(alpha1 - interface_.phase1().residualAlpha())
       *pos(alpha2 - interface_.phase2().residualAlpha());

    if (debug & mesh_.time().writeTime())
    {
        U1I_.write();
        U2I_.write();
        UrI_.write();
        UI_.write();
        yI_.write();
    }
}


void Foam::underResolutionModel::correctTauI()
{
    volScalarField y1I(posPart(yI_));
    volScalarField y2I(-negPart(yI_));
    fvc::spread(y1I, interface_.phase1(), nLayers_, 0.001);
    fvc::spread(y2I, interface_.phase2(), nLayers_, 0.001);
    y1I *= pos(y2I - dimensionedScalar(dimLength, small));
    y2I *= pos(y1I - dimensionedScalar(dimLength, small));

    const volScalarField deltaU1
    (
        pos(y1I - dimensionedScalar(dimLength, small))*mag(U1I_ - UI_)
    );
    const volScalarField deltaU2
    (
        pos(y2I - dimensionedScalar(dimLength, small))*mag(U2I_ - UI_)
    );

    const volScalarField nu1(interface_.phase1().fluidThermo().nu());
    const volScalarField nu2(interface_.phase2().fluidThermo().nu());

    volScalarField uStar1(wallLawViscous(y1I, deltaU1, nu1));
    volScalarField uStar2(wallLawViscous(y2I, deltaU2, nu2));

    volScalarField yPlus1(y1I*uStar1/nu1);
    volScalarField yPlus2(y2I*uStar2/nu2);

    for (label i = 0; i < nCorrectors_; i++)
    {
        if (interfaceRoughness_)
        {
            roughness_ =
                interface_.phase1().name() == lightPhase_.name()
              ? roughness(y1I, uStar1)
              : roughness(y2I, uStar2);

            const volScalarField Rer1
            (
                IOobject::groupName("Rer", interface_.name()) + ":"
              + interface_.phase1().name(),
                roughness_*uStar1/nu1
            );
            const volScalarField Rer2
            (
                IOobject::groupName("Rer", interface_.name()) + ":"
              + interface_.phase2().name(),
                roughness_*uStar2/nu2
            );

            uStar1 = wallLaw(yPlus1, y1I, deltaU1, nu1, Rer1);
            uStar2 = wallLaw(yPlus2, y2I, deltaU2, nu2, Rer2);
        }
        else
        {
            uStar1 = wallLawSmooth(yPlus1, y1I, deltaU1, nu1);
            uStar2 = wallLawSmooth(yPlus2, y2I, deltaU2, nu2);
        }

        tauI_ =
            tau
            (
                uStar1,
                uStar2,
                interface_.phase1().rho(),
                interface_.phase2().rho()
            );

        yPlus1 = y1I*uStar1/nu1;
        yPlus2 = y2I*uStar2/nu2;
    }

    if (debug & mesh_.time().writeTime()) tauI_.write();
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::underResolutionModel::underResolutionModel
(
    const dictionary& dict,
    const phaseInterface& interface,
    const bool registerObject
)
:
    regIOobject
    (
        IOobject
        (
            IOobject::groupName(typeName, interface.name()),
            interface.mesh().time().name(),
            interface.mesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            registerObject
        )
    ),

    lightPhase_
    (
        gAverage(interface.phase1().rho()())
      < gAverage(interface.phase2().rho()())
      ? interface.phase1()
      : interface.phase2()
    ),

    heavyPhase_
    (
        gAverage(interface.phase1().rho()())
      < gAverage(interface.phase2().rho()())
      ? interface.phase2()
      : interface.phase1()
    ),

    interfaceRoughness_
    (
        dict.lookupOrDefault<Switch>("interfaceRoughness", false)
    ),

    nCorrectors_(dict.lookupOrDefault<label>("nCorrectors", 3)),

    Cs_(10.67),

    Bs_(4.9),

    kappa_(0.41),

    RerS_(2.25),

    RerR_(90),

    Br_(8.5),

    beta1_(0.05),

    beta2_(5.0),

    q310_(0.15),

    q320_(0.05),

    DeltaMu_(0.01),

    gradAlphaI_
    (
        IOobject
        (
            IOobject::groupName("gradAlphaI", interface.name()),
            interface.mesh().time().name(),
            interface.mesh()
        ),
        interface.mesh(),
        dimensionedVector(inv(dimLength), Zero)
    ),

    yI_
    (
        IOobject
        (
            IOobject::groupName("yI", interface.name()),
            interface.mesh().time().name(),
            interface.mesh()
        ),
        interface.mesh(),
        dimensionedScalar(dimLength, Zero)
    ),

    U1I_
    (
        IOobject
        (
            IOobject::groupName("U1I", interface.name()),
            interface.mesh().time().name(),
            interface.mesh()
        ),
        interface.mesh(),
        dimensionedVector(dimVelocity, Zero)
    ),

    U2I_
    (
        IOobject
        (
            IOobject::groupName("U2I", interface.name()),
            interface.mesh().time().name(),
            interface.mesh()
        ),
        interface.mesh(),
        dimensionedVector(dimVelocity, Zero)
    ),

    UrI_
    (
        IOobject
        (
            IOobject::groupName("UrI", interface.name()),
            interface.mesh().time().name(),
            interface.mesh()
        ),
        interface.mesh(),
        dimensionedVector(dimVelocity, Zero)
    ),

    UI_
    (
        IOobject
        (
            IOobject::groupName("UI", interface.name()),
            interface.mesh().time().name(),
            interface.mesh()
        ),
        interface.mesh(),
        dimensionedVector(dimVelocity, Zero)
    ),

    tauI_
    (
        IOobject
        (
            IOobject::groupName("tauI", interface.name()),
            interface.mesh().time().name(),
            interface.mesh()
        ),
        interface.mesh(),
        dimensionedScalar(dimPressure, Zero)
    ),

    roughness_
    (
        IOobject
        (
            IOobject::groupName("roughness", interface.name()),
            interface.mesh().time().name(),
            interface.mesh()
        ),
        interface.mesh(),
        dimensionedScalar(dimLength, Zero)
    ),

    interface_(interface.phase1(), interface.phase2()),

    mesh_(interface.mesh()),

    underResolution_
    (
        IOobject
        (
            IOobject::groupName("underResolution", interface.name()),
            mesh_.time().name(),
            mesh_
        ),
        mesh_,
        dimensionedScalar(dimless, Zero)
    ),

    nLayers_(dict.lookupOrDefault<label>("nLayers", 5))
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::underResolutionModel::~underResolutionModel()
{}


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

const Foam::volScalarField& Foam::underResolutionModel::yWall() const
{
    return wallDist::New(mesh_).y();
}


void Foam::underResolutionModel::correct()
{
    correctGradAlphaI();
    correctUIyI();
    correctTauI();
}


bool Foam::underResolutionModel::writeData(Ostream& os) const
{
    return os.good();
}


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