/*---------------------------------------------------------------------------*\
  ==  == ====== ====   ====    |
                    \\     ||  | 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 "addonFilmSurfaceVelocityFvPatchVectorField.H"
#include "compressibleMomentumTransportModel.H"
#include "mappedFvPatchBaseBase.H"
#include "addToRunTimeSelectionTable.H"

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

Foam::addonFilmSurfaceVelocityFvPatchVectorField::
addonFilmSurfaceVelocityFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const dictionary& dict
)
:
    mixedFvPatchField<vector>(p, iF, dict, false),
    phaseNames_(dict.lookup<wordList>("phases")),
    Cs_(phaseNames_.size())
{
    refValue() = Zero;
    refGrad() = Zero;
    valueFraction() = 0;

    forAll(phaseNames_, i)
    {
        Cs_[i] =
            dict.lookupOrDefault<scalar>
            (
                IOobject::groupName("Cs", phaseNames_[i]),
                0
            );
    }

    if (dict.found("value"))
    {
        fvPatchVectorField::operator=
        (
            vectorField("value", iF.dimensions(), dict, p.size())
        );
    }
    else
    {
        // If the value entry is not present initialise to zero-gradient
        fvPatchVectorField::operator=(patchInternalField());
    }
}


Foam::addonFilmSurfaceVelocityFvPatchVectorField::
addonFilmSurfaceVelocityFvPatchVectorField
(
    const addonFilmSurfaceVelocityFvPatchVectorField& ptf,
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const fieldMapper& mapper
)
:
    mixedFvPatchField<vector>(ptf, p, iF, mapper),
    phaseNames_(ptf.phaseNames_),
    Cs_(ptf.Cs_)
{}


Foam::addonFilmSurfaceVelocityFvPatchVectorField::
addonFilmSurfaceVelocityFvPatchVectorField
(
    const addonFilmSurfaceVelocityFvPatchVectorField& ptf,
    const DimensionedField<vector, volMesh>& iF
)
:
    mixedFvPatchField<vector>(ptf, iF),
    phaseNames_(ptf.phaseNames_),
    Cs_(ptf.Cs_)
{}


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

void Foam::addonFilmSurfaceVelocityFvPatchVectorField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    // Since we're inside initEvaluate/evaluate there might be processor
    // comms underway. Change the tag we use.
    const int oldTag = UPstream::msgType();
    UPstream::msgType() = oldTag + 1;

    // Get the coupling information from the mappedPatchBase
    const mappedFvPatchBaseBase& mapper =
        mappedFvPatchBaseBase::getMap(patch());
    const fvPatch& patchNbr = mapper.nbrFvPatch();

    // Retrieve the normal vector of the surface patch
    const vectorField n(patch().nf());

    // Lookup the momentum transport model
    const compressibleMomentumTransportModel& transportModel =
        db().lookupType<compressibleMomentumTransportModel>();

    // Patch laminar dynamic viscosity divided by delta
    const tmp<scalarField> muEffByDelta
    (
        transportModel.rho().boundaryField()[patch().index()]
       *transportModel.nuEff(patch().index())
       *patch().deltaCoeffs()
    );

    // Initialize variables for weighted values
    scalarField alphapNbrTot(patch().size(), 0.0);
    scalarField Ds(patch().size(), 0.0);
    vectorField UMixpNbr(patch().size(), vector::zero);

    forAll(phaseNames_, i)
    {
        const word phase(phaseNames_[i]);

        // Neighbour phase patch face-cell velocity
        const vectorField UphasepNbr
        (
            patchNbr.lookupPatchField<volVectorField, scalar>
            (
                IOobject::groupName("U", phase)
            ).patchInternalField()
        );

        // Remove the normal component of the surface vel
        vectorField UphasepNbrProj = mapper.fromNeighbour(UphasepNbr);
        UphasepNbrProj -= n*(n & UphasepNbrProj);

        // Lookup the neighbour momentum transport model
        const compressibleMomentumTransportModel& transportModelNbr =
            mapper.nbrMesh().lookupType<compressibleMomentumTransportModel>
            (
                phase
            );

         // Get the neighbour patch density
        const tmp<scalarField> rhophasepNbr
        (
            mapper.fromNeighbour
            (
                transportModelNbr.rho().boundaryField()[patchNbr.index()]
            )
        );

         // Get the neighbour patch volume fraction
        const tmp<scalarField> tmpAlphaNbr
        (
            mapper.fromNeighbour
            (
                patchNbr.lookupPatchField<volScalarField, scalar>
                (
                    IOobject::groupName("alpha", phase)
                ).patchInternalField()
            )
        );
        scalarField alphapNbr = tmpAlphaNbr();

        if (Cs_[i] > 0)
        {
            // Calculate the drag coefficient from the drag constant
            // and the magnitude of the velocity difference
            Ds += alphapNbr*Cs_[i]*rhophasepNbr()*mag(UphasepNbrProj - *this);
        }
        else
        {
            // Get neighbour patch laminar dynamic viscosity divided by delta
            const tmp<scalarField> muEffByDeltapNbr
            (
                mapper.fromNeighbour
                (
                    transportModelNbr.rho().boundaryField()[patchNbr.index()]
                   *transportModelNbr.nuEff(patchNbr.index())
                   *patchNbr.deltaCoeffs()
                )
            );

            Ds += alphapNbr*muEffByDeltapNbr();
        }

        alphapNbrTot += alphapNbr;
        UMixpNbr += alphapNbr*UphasepNbrProj;
    }

    UMixpNbr = UMixpNbr/alphapNbrTot;
    Ds = Ds/alphapNbrTot;
    refValue() = UMixpNbr;

    // Calculate the value-fraction from the balance between the
    // external fluid drag and/or stress and internal film stresses
    valueFraction() = Ds/(muEffByDelta + Ds);

    mixedFvPatchField<vector>::updateCoeffs();

    // Restore tag
    UPstream::msgType() = oldTag;
}


void Foam::addonFilmSurfaceVelocityFvPatchVectorField::write(Ostream& os) const
{
    fvPatchField<vector>::write(os);

    writeEntry(os, "phases", phaseNames_);
    forAll(phaseNames_, i)
    {
        writeEntry(os, IOobject::groupName("Cs", phaseNames_[i]), Cs_[i]);
    }
    writeEntry(os, "value", *this);
}


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

namespace Foam
{
   makePatchTypeField
   (
       fvPatchVectorField,
       addonFilmSurfaceVelocityFvPatchVectorField
   );
}


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