ENH: added weightedFlux interpolation scheme (#1388)

- This scheme is useful to calculate the face interpolation values for
  the Gauss gradient when the diffussion coefficient is discontinuous
  across a face. This sheme is used for Gauss grad.

src/finiteVolume/Make/files

@@ -317,6 +317,7 @@ $(schemes)/pointLinear/pointLinear.C
 $(schemes)/midPoint/midPoint.C
 $(schemes)/downwind/downwind.C
 $(schemes)/weighted/weighted.C
+$(schemes)/weightedFlux/weightedFlux.C
 $(schemes)/cubic/cubic.C
 $(schemes)/skewCorrected/skewCorrectionVectors.C
 $(schemes)/skewCorrected/skewCorrected.C

src/finiteVolume/interpolation/surfaceInterpolation/schemes/weightedFlux/weightedFlux.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 Norbert Weber, HZDR
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "weightedFlux.H"
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+template<class Type>
+void Foam::weightedFlux<Type>::clearOut()
+{
+    deleteDemandDrivenData(oDelta_);
+    deleteDemandDrivenData(nDelta_);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::weightedFlux<Type>::~weightedFlux()
+{
+    clearOut();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::weightedFlux<Type>::makeDeltas() const
+{
+    const fvMesh& mesh = this->mesh();
+
+    oDelta_ = new surfaceScalarField
+    (
+        IOobject
+        (
+            "oDelta",
+            mesh.pointsInstance(),
+            mesh
+        ),
+        mesh,
+        dimLength
+    );
+    auto& oDelta = *oDelta_;
+
+    nDelta_ = new surfaceScalarField
+    (
+        IOobject
+        (
+            "nDelta",
+            mesh.pointsInstance(),
+            mesh
+        ),
+        mesh,
+        dimLength
+    );
+    auto& nDelta = *nDelta_;
+
+    const labelUList& owner = mesh.owner();
+    const labelUList& neighbour = mesh.neighbour();
+
+    const surfaceVectorField n = mesh.Sf()/mesh.magSf();
+
+    const vectorField& C = mesh.cellCentres();
+    const vectorField& Cf = mesh.faceCentres();
+
+    // all distances are NORMAL to the face,
+    // as in the weighting factors in surfaceInterpolation.C
+    forAll(owner, facei)
+    {
+        oDelta[facei] =
+            mag(n[facei] & (C[owner[facei]] - Cf[facei]));
+        nDelta[facei] =
+            mag(n[facei] & (C[neighbour[facei]] - Cf[facei]));
+    }
+
+    const fvPatchList& patches = mesh.boundary();
+
+    forAll(patches, patchi)
+    {
+        const fvPatch& currPatch = mesh.boundary()[patchi];
+
+        // Patch normal vector
+        const vectorField nPatch = currPatch.Sf()/currPatch.magSf();
+
+        // Processor patch
+        if (currPatch.coupled())
+        {
+            const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
+            const vectorField& pCf = mesh.Cf().boundaryField()[patchi];
+
+            forAll(pOwner, facei)
+            {
+                const label own = pOwner[facei];
+
+                // All distances are NORMAL to the face
+                oDelta.boundaryFieldRef()[patchi][facei] =
+                    mag(nPatch[facei] & (pCf[facei] - C[own]));
+            }
+
+            // Weight = delta_neighbour / delta in ORTHOGONAL direction,
+            nDelta.boundaryFieldRef()[patchi] =
+                currPatch.weights()*oDelta.boundaryFieldRef()[patchi]
+               /(scalar(1) - currPatch.weights());
+        }
+        else
+        {
+            const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
+            const vectorField& pCf = mesh.Cf().boundaryField()[patchi];
+
+            forAll(pOwner, facei)
+            {
+                const label own = pOwner[facei];
+
+                // All distances are NORMAL to the face!
+                oDelta.boundaryFieldRef()[patchi][facei] =
+                    mag(nPatch[facei] & (pCf[facei] - C[own]));
+
+                nDelta.boundaryFieldRef()[patchi][facei] =
+                    mag(nPatch[facei] & (pCf[facei] - C[own]));
+            }
+        }
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>>
+Foam::weightedFlux<Type>::interpolate
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    const fvMesh& mesh = vf.mesh();
+    const surfaceScalarField& oDelta = weightedFlux<Type>::oDelta();
+    const surfaceScalarField& nDelta = weightedFlux<Type>::nDelta();
+
+    auto tresult = tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>::New
+    (
+        IOobject
+        (
+            "weightedFlux::interpolate(" + vf.name() + ')',
+            mesh.time().timeName(),
+            mesh
+        ),
+        mesh,
+        vf.dimensions()
+    );
+    auto& result = tresult.ref();
+
+    const labelUList& owner = mesh.owner();
+    const labelUList& neighbour = mesh.neighbour();
+
+    forAll(result, facei)
+    {
+        const scalar sigmaDeltaO = sigma_[owner[facei]]/oDelta[facei];
+        const scalar sigmaDeltaN = sigma_[neighbour[facei]]/nDelta[facei];
+
+        result[facei] =
+            (vf[owner[facei]]*sigmaDeltaO + vf[neighbour[facei]]*sigmaDeltaN)
+           /(sigmaDeltaO + sigmaDeltaN);
+    }
+
+    // Interpolate patches
+    auto& bfld = result.boundaryFieldRef();
+
+    forAll(bfld, patchi)
+    {
+        fvsPatchField<Type>& pfld = bfld[patchi];
+
+        // If not coupled - simply copy the boundary values of the field
+        if (!pfld.coupled())
+        {
+            pfld = vf.boundaryField()[patchi];
+        }
+        else
+        {
+            // e.g. processor patches have to calculated separately
+
+            const labelUList& pOwner = mesh.boundary()[patchi].faceCells();
+            scalarField sigmaN =
+                sigma_.boundaryField()[patchi].patchNeighbourField();
+
+            Field<Type> vfO = vf.boundaryField()[patchi].patchInternalField();
+            Field<Type> vfN = vf.boundaryField()[patchi].patchNeighbourField();
+
+            forAll(pOwner, facei)
+            {
+                const label own = pOwner[facei];
+
+                const scalar sigmaDeltaO =
+                    sigma_[own]/oDelta.boundaryField()[patchi][facei];
+
+                const scalar sigmaDeltaN =
+                    sigmaN[facei]/nDelta.boundaryField()[patchi][facei];
+
+                pfld[facei] =
+                    (vfO[facei]*sigmaDeltaO + vfN[facei]*sigmaDeltaN)
+                   /(sigmaDeltaO + sigmaDeltaN);
+            }
+        }
+    }
+
+    return tresult;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makeSurfaceInterpolationScheme(weightedFlux)
+}
+
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/schemes/weightedFlux/weightedFlux.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 Norbert Weber, HZDR
+-------------------------------------------------------------------------------
+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::weightedFlux
+
+Description
+    Weighted flux interpolation scheme class.
+
+    This scheme is used to compute fluxes with variable diffusivity or
+    conductivity, as e.g.
+    - a thermal flux: lambda*grad(T)
+    - a mass flux: D*grad(u)
+    - an electric current: -sigma*grad(potential)
+
+    When using the Gauss theorem to compute a gradient, cell centred values
+    need to be interpolated to the faces. Using this scheme, temperature (T)
+    is weighted by thermal conductivity when being interpolated. Similarly,
+    velocity is weighted by diffusivity (D) and the electric potential by
+    the electric conductivity (sigma). Lambda, D or sigma are read from the
+    object registry - the names need to be specified in fvSchemes as e.g.
+
+    \verbatim
+        gradSchemes
+        {
+            grad(T)             Gauss weightedFlux "lambda";
+            grad(u)             Gauss weightedFlux "D";
+            grad(potential)     Gauss weightedFlux "sigma";
+        }
+    \endverbatim
+
+    For more details, see equation 16 and 17 in
+    \verbatim
+        Weber, N., Beckstein, P., Galindo, V., Starace, M. & Weier, T. (2018).
+        Electro-vortex flow simulation using coupled meshes.
+        Computers and Fluids 168, 101-109.
+        doi:10.1016/j.compfluid.2018.03.047
+        https://arxiv.org/pdf/1707.06546.pdf
+    \endverbatim
+
+Note
+    For support, contact Norbert.Weber@hzdr.de
+
+SourceFiles
+    weightedFlux.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef weightedFlux_H
+#define weightedFlux_H
+
+#include "surfaceInterpolationScheme.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class weightedFlux Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class weightedFlux
+:
+    public surfaceInterpolationScheme<Type>
+{
+    // Private Data
+
+        //- Const reference to step-wise pre-gradient factor field
+        const volScalarField& sigma_;
+
+    // Demand-driven data
+
+        //- Face to owner cell distance
+        mutable surfaceScalarField* oDelta_;
+
+        //- Face to neighbour cell distance
+        mutable surfaceScalarField* nDelta_;
+
+
+    // Private Member Functions
+
+        //- Compute face-owner and face-neighbour distance
+        void makeDeltas() const;
+
+        //- No copy assignment
+        void operator=(const weightedFlux&) = delete;
+
+
+protected:
+
+    // Protected Member Functions
+
+        // Storage management
+
+            //- Clear all fields
+            void clearOut();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("weightedFlux");
+
+
+    // Constructors
+
+        //- Construct from Istream.
+        //  The name of the flux field is read from the Istream and looked-up
+        //  from the mesh objectRegistry
+        weightedFlux
+        (
+            const fvMesh& mesh,
+            Istream& is
+        )
+        :
+            surfaceInterpolationScheme<Type>(mesh),
+            sigma_(this->mesh().objectRegistry::template
+                lookupObject<volScalarField>(word(is))),
+            oDelta_(nullptr),
+            nDelta_(nullptr)
+        {}
+
+        //- Construct from faceFlux and Istream
+        weightedFlux
+        (
+            const fvMesh& mesh,
+            const surfaceScalarField& faceFlux,
+            Istream& is
+        )
+        :
+            surfaceInterpolationScheme<Type>(mesh),
+            sigma_(this->mesh().objectRegistry::template
+                lookupObject<volScalarField>(word(is))),
+            oDelta_(nullptr),
+            nDelta_(nullptr)
+        {}
+
+
+    //- Destructor
+    ~weightedFlux();
+
+
+    // Member Functions
+
+        //- Return the interpolation weighting factors
+        tmp<surfaceScalarField> weights
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const
+        {
+            return this->mesh().surfaceInterpolation::weights();
+        }
+
+        //- Return the distance between face and owner cell
+        const surfaceScalarField& oDelta() const
+        {
+            if (!oDelta_)
+            {
+                makeDeltas();
+            }
+
+            return *oDelta_;
+        }
+
+        //- Return the distance between face and neighbour cell
+        const surfaceScalarField& nDelta() const
+        {
+            if (!nDelta_)
+            {
+                makeDeltas();
+            }
+
+            return *nDelta_;
+        }
+
+        //- Interpolate the cell values to faces
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+        interpolate
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //