Clean-up of Gauss div and least squares grad

src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivScheme.C

@@ -79,7 +79,7 @@ tmp
 {
     FatalErrorIn
     (
-       "tmp<BlockLduSystem> fvmUDiv\n"
+        "tmp<BlockLduSystem> gaussDivScheme<Type>::fvmUDiv\n"
         "(\n"
         "    GeometricField<Type, fvPatchField, volMesh>&"
         ")\n"
@@ -109,7 +109,7 @@ tmp
 {
     FatalErrorIn
     (
-       "tmp<BlockLduSystem> fvmUDiv\n"
+        "tmp<BlockLduSystem> gaussDivScheme<Type>::fvmUDiv\n"
         "(\n"
         "    const surfaceScalarField& flux"
         "    const GeometricField<Type, fvPatchField, volMesh>&"

src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/gaussDivSchemes.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "gaussDivScheme.H"
+#include "vectorGaussDivScheme.H"
 #include "fvMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/finiteVolume/finiteVolume/divSchemes/gaussDivScheme/vectorGaussDivScheme.C

@@ -0,0 +1,222 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Specialisation of gaussDivScheme for vectors. Needed for implicit fvmDiv
+    operator for block coupled systems.
+
+\*---------------------------------------------------------------------------*/
+
+#include "vectorGaussDivScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<>
+tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmUDiv
+(
+    const GeometricField<vector, fvPatchField, volMesh>& vf
+) const
+{
+    // Get weights
+    surfaceScalarField weights = this->tinterpScheme_().weights(vf);
+
+    const fvMesh& mesh = vf.mesh();
+
+    tmp<BlockLduSystem<vector, scalar> > tbs
+    (
+        new BlockLduSystem<vector, scalar>(mesh)
+    );
+    BlockLduSystem<vector, scalar>& bs = tbs();
+    scalarField& source = bs.source();
+
+    // Grab ldu parts of block matrix as linear always
+    CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
+    CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
+    CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
+
+    // Internal field
+    const vectorField& SfIn = mesh.Sf().internalField();
+    const scalarField& wIn = weights.internalField();
+
+    l = -wIn*SfIn;
+    u = l + SfIn;
+    bs.negSumDiag();
+
+    // Boundary contributions
+    forAll(vf.boundaryField(), patchI)
+    {
+        const fvPatchVectorField& pf = vf.boundaryField()[patchI];
+        const fvPatch& patch = pf.patch();
+        const vectorField& Sf = patch.Sf();
+        const fvsPatchScalarField& pw = weights.boundaryField()[patchI];
+        const unallocLabelList& fc = patch.faceCells();
+
+        const vectorField internalCoeffs(pf.valueInternalCoeffs(pw));
+
+        // Diag contribution
+        forAll(pf, faceI)
+        {
+            d[fc[faceI]] += cmptMultiply(internalCoeffs[faceI], Sf[faceI]);
+        }
+
+        if (patch.coupled())
+        {
+            CoeffField<vector>::linearTypeField& pcoupleUpper =
+                bs.coupleUpper()[patchI].asLinear();
+            CoeffField<vector>::linearTypeField& pcoupleLower =
+                bs.coupleLower()[patchI].asLinear();
+
+            const vectorField pcl = -pw*Sf;
+            const vectorField pcu = pcl + Sf;
+
+            // Coupling  contributions
+            pcoupleLower -= pcl;
+            pcoupleUpper -= pcu;
+        }
+        else
+        {
+            const vectorField boundaryCoeffs(pf.valueBoundaryCoeffs(pw));
+
+            // Boundary contribution
+            forAll(pf, faceI)
+            {
+                source[fc[faceI]] -= boundaryCoeffs[faceI] & Sf[faceI];
+            }
+        }
+    }
+
+    // Interpolation schemes with corrections not accounted for
+
+    return tbs;
+}
+
+
+template<>
+tmp<BlockLduSystem<vector, scalar> > gaussDivScheme<vector>::fvmUDiv
+(
+    const surfaceScalarField& flux,
+    const GeometricField<vector, fvPatchField, volMesh>& vf
+) const
+{
+    // Get weights
+    surfaceScalarField weights = this->tinterpScheme_().weights(vf);
+
+    const fvMesh& mesh = vf.mesh();
+
+    tmp<BlockLduSystem<vector, scalar> > tbs
+    (
+        new BlockLduSystem<vector, scalar>(mesh)
+    );
+    BlockLduSystem<vector, scalar>& bs = tbs();
+    scalarField& source = bs.source();
+
+    // Grab ldu parts of block matrix as linear always
+    CoeffField<vector>::linearTypeField& d = bs.diag().asLinear();
+    CoeffField<vector>::linearTypeField& u = bs.upper().asLinear();
+    CoeffField<vector>::linearTypeField& l = bs.lower().asLinear();
+
+    // Internal field
+    const vectorField& SfIn = mesh.Sf().internalField();
+    const scalarField& wIn = weights.internalField();
+    const scalarField& fluxIn = flux.internalField();
+
+    l = -wIn*fluxIn*SfIn;
+    u = l + fluxIn*SfIn;
+    bs.negSumDiag();
+
+    // Boundary contributions
+    forAll(vf.boundaryField(), patchI)
+    {
+        const fvPatchVectorField& pf = vf.boundaryField()[patchI];
+        const fvPatch& patch = pf.patch();
+        const vectorField& Sf = patch.Sf();
+        const fvsPatchScalarField& pw = weights.boundaryField()[patchI];
+        const unallocLabelList& fc = patch.faceCells();
+
+        const scalarField& pFlux = flux.boundaryField()[patchI];
+
+        const vectorField internalCoeffs(pf.valueInternalCoeffs(pw));
+
+        // Diag contribution
+        forAll(pf, faceI)
+        {
+            d[fc[faceI]] += cmptMultiply
+            (
+                internalCoeffs[faceI],
+                pFlux[faceI]*Sf[faceI]
+            );
+        }
+
+        if (patch.coupled())
+        {
+            CoeffField<vector>::linearTypeField& pcoupleUpper =
+                bs.coupleUpper()[patchI].asLinear();
+            CoeffField<vector>::linearTypeField& pcoupleLower =
+                bs.coupleLower()[patchI].asLinear();
+
+            const vectorField pcl = -pw*pFlux*Sf;
+            const vectorField pcu = pcl + pFlux*Sf;
+
+            // Coupling  contributions
+            pcoupleLower -= pcl;
+            pcoupleUpper -= pcu;
+        }
+        else
+        {
+            const vectorField boundaryCoeffs(pf.valueBoundaryCoeffs(pw));
+
+            // Boundary contribution
+            forAll(pf, faceI)
+            {
+                source[fc[faceI]] -=
+                (
+                    boundaryCoeffs[faceI] & (pFlux[faceI]*Sf[faceI])
+                );
+            }
+        }
+    }
+
+    // Interpolation schemes with corrections not accounted for
+
+    return tbs;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrads.C

@@ -27,6 +27,7 @@ Description
 
 #include "fvMesh.H"
 #include "gaussGrad.H"
+#include "scalarGaussGrad.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresGrad.C

@@ -94,15 +94,26 @@ leastSquaresGrad<Type>::grad
     const unallocLabelList& own = mesh.owner();
     const unallocLabelList& nei = mesh.neighbour();
 
+    // Get access to internal field
+
+    const Field<Type>& vsfIn = vsf.internalField();
+
+    Field<GradType>& lsGradIn = lsGrad.internalField();
+
+    const vectorField& ownLsIn = ownLs.internalField();
+    const vectorField& neiLsIn = neiLs.internalField();
+
+    register label ownFaceI, neiFaceI;
+
     forAll (own, facei)
     {
-        register label ownFaceI = own[facei];
+        ownFaceI = own[facei];
-        register label neiFaceI = nei[facei];
+        neiFaceI = nei[facei];
 
-        Type deltaVsf = vsf[neiFaceI] - vsf[ownFaceI];
+        Type deltaVsf = vsfIn[neiFaceI] - vsfIn[ownFaceI];
 
-        lsGrad[ownFaceI] += ownLs[facei]*deltaVsf;
+        lsGradIn[ownFaceI] += ownLsIn[facei]*deltaVsf;
-        lsGrad[neiFaceI] -= neiLs[facei]*deltaVsf;
+        lsGradIn[neiFaceI] -= neiLsIn[facei]*deltaVsf;
     }
 
     // Boundary faces
@@ -138,7 +149,6 @@ leastSquaresGrad<Type>::grad
         }
     }
 
-
     lsGrad.correctBoundaryConditions();
     gaussGrad<Type>::correctBoundaryConditions(vsf, lsGrad);
 

src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C

@@ -134,12 +134,23 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
         // Better version of d-vectors: Zeljko Tukovic, 25/Apr/2010
         const vectorField pd = p.delta();
 
+        if (p.coupled())
+        {
             forAll (pd, pFaceI)
             {
                 const vector& d = pd[pFaceI];
                 dd[fc[pFaceI]] += (1.0/magSqr(d))*sqr(d);
             }
         }
+        else
+        {
+            forAll (pd, pFaceI)
+            {
+                const vector& d = pd[pFaceI];
+                dd[fc[pFaceI]] += (1.0/magSqr(d))*sqr(d);
+            }
+        }
+    }
 
     // For easy access of neighbour coupled patch field needed for
     // lsN vectors on implicitly coupled boundaries. VV, 18/June/2014
@@ -221,8 +232,6 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
         }
     }
 
-# if 0
-
     // Coefficient sign correction on least squares vectors
     // The sign of the coefficient must be positive to ensure correct
     // behaviour in coupled systems.  This is checked by evaluating
@@ -408,7 +417,42 @@ void Foam::leastSquaresVectors::makeLeastSquaresVectors() const
         }
     }
 
-#endif
+    // Manual check of least squares vectors
+
+    vectorField sumLsq(mesh().nCells(), vector::zero);
+    vectorField sumSf(mesh().nCells(), vector::zero);
+
+    const vectorField& sfIn = mesh().Sf().internalField();
+
+    // Least squares vectors on internal faces
+    forAll (owner, faceI)
+    {
+        sumLsq[owner[faceI]] += lsPIn[faceI];
+        sumLsq[neighbour[faceI]] += lsNIn[faceI];
+
+        sumSf[owner[faceI]] += sfIn[faceI];
+        sumSf[neighbour[faceI]] -= sfIn[faceI];
+    }
+
+    // Least squares vectors on boundary faces
+    forAll (lsP.boundaryField(), patchI)
+    {
+        const vectorField& patchLsP = lsP.boundaryField()[patchI];
+        const vectorField& patchSf = mesh().Sf().boundaryField()[patchI];
+
+        const unallocLabelList& fc = mesh().boundary()[patchI].faceCells();
+
+        forAll (fc, pFaceI)
+        {
+            //sumLsq[fc[pFaceI]] += 0.5*patchLsP[pFaceI]; // works!!!
+            sumLsq[fc[pFaceI]] += patchLsP[pFaceI];
+
+            sumSf[fc[pFaceI]] += patchSf[pFaceI];
+        }
+    }
+
+    Info<< "sumLsq min = " << min(sumLsq) << " max = " << max(sumLsq) << endl;
+    Info<< "sumSf min = " << min(sumSf) << " max = " << max(sumSf) << endl;
 
     if (debug)
     {

src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/scalarLeastSquaresGrad.C

@@ -90,13 +90,19 @@ tmp<BlockLduSystem<vector, vector> > leastSquaresGrad<scalar>::fvmGrad
     const surfaceVectorField& ownLs = lsv.pVectors();
     const surfaceVectorField& neiLs = lsv.nVectors();
 
+    // Internal field
+    const vectorField& ownLsIn = ownLs.internalField();
+    const vectorField& neiLsIn = neiLs.internalField();
+
+    register label owner, neighbour;
+
     forAll (nei, faceI)
     {
-        register label owner = own[faceI];
+        owner = own[faceI];
-        register label neighbour = nei[faceI];
+        neighbour = nei[faceI];
 
-        u[faceI] = cellVIn[owner]*ownLs[faceI];
+        u[faceI] = cellVIn[owner]*ownLsIn[faceI];
-        l[faceI] = cellVIn[neighbour]*neiLs[faceI];
+        l[faceI] = cellVIn[neighbour]*neiLsIn[faceI];
 
         // Caution - this is NOT negSumDiag(). VV, 17/July/2014.
         d[owner] -= u[faceI];
@@ -133,11 +139,8 @@ tmp<BlockLduSystem<vector, vector> > leastSquaresGrad<scalar>::fvmGrad
             // Coupling  and diagonal contributions
             forAll (pf, faceI)
             {
-                const vector upper = cellVIn[fc[faceI]]*pownLs[faceI];
+                pcoupleUpper[faceI] -= cellVIn[fc[faceI]]*pownLs[faceI];
-                const vector lower = cellVInNei[faceI]*pneiLs[faceI];
+                pcoupleLower[faceI] -= cellVInNei[faceI]*pneiLs[faceI];
-
-                pcoupleUpper[faceI] -= upper;
-                pcoupleLower[faceI] -= lower;
             }
         }
         else
@@ -149,7 +152,9 @@ tmp<BlockLduSystem<vector, vector> > leastSquaresGrad<scalar>::fvmGrad
             forAll (pf, faceI)
             {
                 const label cellI = fc[faceI];
+
                 d[cellI] += cellVIn[cellI]*pownLs[faceI]*internalCoeffs[faceI];
+
                 source[cellI] -= cellVIn[cellI]*pownLs[faceI]*
                     boundaryCoeffs[faceI];
             }