Conservative GGI bridging updates and fixes

Scaling the fluxes and killing non-orthogonal correction vectors in case
bridging is used, both in order to ensure global conservation across bridged GGI
with potentially partially overlapping faces.

src/finiteVolume/Make/files

@@ -119,6 +119,7 @@ $(constraintFvPatchFields)/symmetry/symmetryFvPatchFields.C
 $(constraintFvPatchFields)/wedge/wedgeFvPatchFields.C
 $(constraintFvPatchFields)/wedge/wedgeFvPatchScalarField.C
 $(constraintFvPatchFields)/wedge/wedgeFvPatchVectorNFields.C
+$(constraintFvPatchFields)/ggi/ggiFvPatchScalarField.C
 $(constraintFvPatchFields)/ggi/ggiFvPatchFields.C
 $(constraintFvPatchFields)/ggi/ggiFvPatchVectorNFields.C
 $(constraintFvPatchFields)/jumpGgi/jumpGgiFvPatchFields.C

src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C

@@ -38,6 +38,7 @@ Note on parallelisation
 #include "ggiFvPatchField.H"
 #include "symmTransformField.H"
 #include "coeffFields.H"
+#include "fvMatrices.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -520,12 +521,6 @@ void Foam::ggiFvPatchField<Type>::manipulateValueCoeffs
 
         // Scale partially overlapping boundary coeffs to ensure conservation
         ggiPatch_.shadow().scalePartialFaces(slaveBC);
-
-//        Info<< "MANIPULATING VALUE COEFFS" << endl;
-//        Info<< "Slave internal coeffs: " << slaveIC << nl
-//            << "Slave boundary coeffs: " << slaveBC << nl
-//            << "Master internal coeffs: " << masterIC << nl
-//            << "Master boundary coeffs: " << masterBC << endl;
     }
 }
 
@@ -552,12 +547,6 @@ void Foam::ggiFvPatchField<Type>::manipulateGradientCoeffs
         Field<Type>& slaveIC = matrix.internalCoeffs()[sPatchI];
         Field<Type>& slaveBC = matrix.boundaryCoeffs()[sPatchI];
 
-//        Info<< "ORIGINAL COEFFS" << endl;
-//        Info<< "Slave internal coeffs: " << slaveIC << nl
-//            << "Slave boundary coeffs: " << slaveBC << nl
-//            << "Master internal coeffs: " << masterIC << nl
-//            << "Master boundary coeffs: " << masterBC << endl;
-
         // Get surface area magnitudes
         const scalarField& magSfMaster = ggiPatch_.magSf();
         const scalarField& magSfSlave = ggiPatch_.shadow().magSf();
@@ -579,12 +568,6 @@ void Foam::ggiFvPatchField<Type>::manipulateGradientCoeffs
         const Field<Type> slaveInterpolatedBC =
             magSfMaster*ggiPatch_.interpolate(slaveBC/magSfSlave);
 
-//        Info<< "INTERPOLATED COEFFS" << endl;
-//        Info<< "Slave internal coeffs: " << slaveInterpolatedIC << nl
-//            << "Slave boundary coeffs: " << slaveInterpolatedBC << nl
-//            << "Master internal coeffs: " << masterInterpolatedIC << nl
-//            << "Master boundary coeffs: " << masterInterpolatedBC << endl;
-
 
         // Set partially covered master coeffs using slave data
         ggiPatch_.setPartialFaces(slaveInterpolatedBC, masterIC);
@@ -610,12 +593,6 @@ void Foam::ggiFvPatchField<Type>::manipulateGradientCoeffs
         // Add to partially overlapping slave internal coeffs to ensure
         // conservation. Note: boundary coeffs must be negated
         ggiPatch_.shadow().addToPartialFaces((-slaveBC)(), slaveIC);
-
-//        Info<< "MANIPULATING GRADIENT COEFFS" << endl;
-//        Info<< "Slave internal coeffs: " << slaveIC << nl
-//            << "Slave boundary coeffs: " << slaveBC << nl
-//            << "Master internal coeffs: " << masterIC << nl
-//            << "Master boundary coeffs: " << masterBC << endl;
     }
 }
 

src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.H

@@ -290,6 +290,8 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "ggiFvPatchScalarField.H"
+
 #ifdef NoRepository
 #   include "ggiFvPatchField.C"
 #endif

src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchScalarField.C

@@ -25,13 +25,14 @@ Author
     Vuko Vukcevic, Wikki Ltd.  All rights reserved
 
 Description
-    Specilalisation of manipulateMatrix member function for scalars needed to
+    Specilalisation of patchFlux member function for scalars needed to
     ensure conservation across GGI patches that are partially covered if the
     bridge overlap switch is on.
 
 \*---------------------------------------------------------------------------*/
 
 #include "ggiFvPatchField.H"
+#include "surfaceFields.H"
 #include "fvScalarMatrix.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -42,51 +43,87 @@ namespace Foam
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<>
-void Foam::ggiFvPatchField<scalar>::manipulateMatrix(fvMatrix<scalar>& matrix)
+void Foam::ggiFvPatchField<scalar>::patchFlux
+(
+    GeometricField<scalar, fvsPatchField, surfaceMesh>& flux,
+    const fvMatrix<scalar>& matrix
+) const
 {
-    Info<< "IN MANIPULATE MATRIX!" << endl;
+    // Since we have adjusted the internal/boundary coefficients in the
+    // manipulateMatrix member function below, we must not use
+    // patchNeighbourField for reconstructing the flux. We only need to use
+    // interpolated shadow field. VV, 6/Mar/2018.
+
+    // Get patch ID
+    const label patchI = this->patch().index();
+
+    // Get internal field
+    const scalarField& iField = this->internalField();
+
+    // Get shadow face-cells and assemble shadow field
+    const unallocLabelList& sfc = ggiPatch_.shadow().faceCells();
+
+    scalarField sField(sfc.size());
+    forAll (sField, i)
+    {
+        sField[i] = iField[sfc[i]];
+    }
+
+    // Interpolate shadow to this side. Note: must not bridge since internal
+    // coeffs and boundary coeffs take it into account
+    scalarField neighbourField(ggiPatch_.interpolate(sField));
 
-    // Conservative treatment for bridged overlap
     if (ggiPatch_.bridgeOverlap())
     {
-        // Get this patch and shadow patch index
+        // Only set fully uncovered faces. Partially covered faces taken into
-        const label patchI = this->patch().index();
+        // account by manipulating value and gradient matrix coefficients. Note:
-        const label sPatchI = ggiPatch_.shadowIndex();
+        // mirror field is the same as patch internal field
+        ggiPatch_.setUncoveredFaces
+        (
+            this->patchInternalField()(),
+            neighbourField
+        );
+    }
 
-        // Get all matrix coefficients
+    // Calculate the flux with correct neighbour field (fully uncovered faces
-        scalarField& thisIC = matrix.internalCoeffs()[patchI];
+    // bridged, while partially uncovered faces taken into account by
-        scalarField& thisBC = matrix.boundaryCoeffs()[patchI];
+    // manipulating value and gradient matrix coefficients in order to ensure
-        scalarField& shadowIC = matrix.internalCoeffs()[sPatchI];
+    // conservation for both convection and diffusion part across partially
-        scalarField& shadowBC = matrix.boundaryCoeffs()[sPatchI];
+    // overlapping faces). VV, 14/Mar/2018.
+    flux.boundaryField()[patchI] =
+        matrix.internalCoeffs()[patchI]*this->patchInternalField()
+      - matrix.boundaryCoeffs()[patchI]*neighbourField;
 
-//        if (!ggiPatch_.master())
+    // Scale the flux on slave patch to ensure global conservation across this
+    // partially overlapping GGI patch. The slight disbalance can happen since
+    // we interpolate the matrix coeffs and field separately, and we should do
+    // it together to ensure conservation. Current code design does not easily
+    // allow this, so this is a meaningful temporary solution when we have
+    // partially overlapping faces. VV, 14/Mar/2018.
+    if (ggiPatch_.bridgeOverlap() && !ggiPatch_.master())
     {
-            // This is master, interpolate shadow on this side
+        // This is slave patch, master already updated the fluxes and we can use
-            const scalarField shadowInterpolatedIC = ggiPatch_.interpolate(shadowIC);
+        // that info to scale the fluxes on this side
-            const scalarField shadowInterpolatedBC = ggiPatch_.interpolate(shadowBC);
 
-            // Set new internal coeffs using the data from the other side and
+        // Get the total flux through master patch
-            // taking into account partially and fully uncovered faces
+        const scalar masterFlux =
-            const scalarField origIC = thisIC;
+            gSum(flux.boundaryField()[ggiPatch_.shadowIndex()]);
-            thisIC = shadowInterpolatedBC;
-            ggiPatch_.scaleForPartialCoverage(origIC, thisIC);
-            ggiPatch_.scaleForPartialCoverage(origIC, thisIC);
 
-            // Set new boundary coeffs using the data from the other side and
+        // Get the total flux through slave patch
-            // taking into account partially and fully uncovered faces
+        const scalar slaveFlux = gSum(flux.boundaryField()[patchI]);
-            const scalarField origBC = thisBC;
-            thisBC = shadowInterpolatedIC;
-            ggiPatch_.scaleForPartialCoverage(origBC, thisBC);
-            ggiPatch_.scaleForPartialCoverage(origBC, thisBC);
 
-            Info<< "This new internal coeffs: " << thisIC << nl
+        // Calculate the scaling factor. Note: negative sign since master and
-                << "This old internal coeffs: " << origIC << nl
+        // slave fluxes have opposite sign
-                << "This shadow boundary coeffs" << shadowBC << nl << endl;
+        const scalar scalingFactor = -masterFlux/(slaveFlux + SMALL);
 
-            Info<< "This new boundary coeffs: " << thisBC << nl
+        // Scale the slave flux to ensure global patch conservation
-                << "This old boundary coeffs: " << origBC << nl
+        flux.boundaryField()[patchI] *= scalingFactor;
-                << "This shadow internale coeffs" << shadowIC << nl << endl;
 
+        if (debug)
+        {
+            Info<< "Scaling flux on patch: " << ggiPatch_.name()
+                << " with " << scalingFactor
+                << " to ensure conservation." << endl;
         }
     }
 }

src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchScalarField.H

@@ -25,7 +25,7 @@ Author
     Vuko Vukcevic, Wikki Ltd.  All rights reserved
 
 Description
-    Specilalisation of manipulateMatrix member function for scalars needed to
+    Specilalisation of patchFlux member function for scalars needed to
     ensure conservation across GGI patches that are partially covered if the
     bridge overlap switch is on.
 
@@ -47,10 +47,11 @@ namespace Foam
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<>
-void Foam::ggiFvPatchField<scalar>::manipulateMatrix
+void Foam::ggiFvPatchField<scalar>::patchFlux
 (
-    fvMatrix<scalar>& matrix
+    GeometricField<scalar, fvsPatchField, surfaceMesh>& flux,
-);
+    const fvMatrix<scalar>& matrix
+) const;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C

@@ -156,6 +156,14 @@ void Foam::ggiFvPatch::makeCorrVecs(vectorField& cv) const
     // Non-orthogonality correction on a ggi interface
     // MB, 7/April/2009
 
+    // No non-orthognal correction if the bridge overlap is switched on to
+    // ensure conservative interpolation for partially overlapping faces
+    if (bridgeOverlap())
+    {
+        cv = vector::zero;
+    }
+    else
+    {
         // Calculate correction vectors on coupled patches
         const scalarField& patchDeltaCoeffs = deltaCoeffs();
 
@@ -165,6 +173,7 @@ void Foam::ggiFvPatch::makeCorrVecs(vectorField& cv) const
         // If non-orthogonality is over 90 deg, kill correction vector
         // HJ, 6/Jan/2011
         cv = pos(patchDeltas & n)*(n - patchDeltas*patchDeltaCoeffs);
+    }
 }