GAMGSolverInterpolate: Correct interpolate for conservation

such that the sum of the values in the fines cells is equal to the value in the coarse cell

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H

@@ -230,7 +230,8 @@ class GAMGSolver
             const lduMatrix& m,
             const FieldField<Field, scalar>& interfaceBouCoeffs,
             const lduInterfaceFieldPtrsList& interfaces,
-            const scalarField& source,
+            const labelList& restrictAddressing,
+            const scalarField& psiC,
             const direction cmpt
         ) const;
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverInterpolate.C

@@ -34,7 +34,8 @@ void Foam::GAMGSolver::interpolate
     const lduMatrix& m,
     const FieldField<Field, scalar>& interfaceBouCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
-    const scalarField& source,
+    const labelList& restrictAddressing,
+    const scalarField& psiC,
     const direction cmpt
 ) const
 {
@@ -80,6 +81,30 @@ void Foam::GAMGSolver::interpolate
     {
         psiPtr[celli] = -ApsiPtr[celli]/(diagPtr[celli]);
     }
+
+    register const label nCCells = psiC.size();
+    scalarField corrC(nCCells, 0);
+    scalarField diagC(nCCells, 0);
+
+    for (register label celli=0; celli<nCells; celli++)
+    {
+        corrC[restrictAddressing[celli]] += diagPtr[celli]*psiPtr[celli];
+        diagC[restrictAddressing[celli]] += diagPtr[celli];
+        //corrC[restrictAddressing[celli]] += psiPtr[celli]/diagPtr[celli];
+        //diagC[restrictAddressing[celli]] += 1.0/diagPtr[celli];
+        //corrC[restrictAddressing[celli]] += psiPtr[celli];
+        //diagC[restrictAddressing[celli]] += 1.0;
+    }
+
+    for (register label ccelli=0; ccelli<nCCells; ccelli++)
+    {
+        corrC[ccelli] = psiC[ccelli] - corrC[ccelli]/diagC[ccelli];
+    }
+
+    for (register label celli=0; celli<nCells; celli++)
+    {
+        psiPtr[celli] += corrC[restrictAddressing[celli]];
+    }
 }
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C

@@ -313,7 +313,7 @@ void Foam::GAMGSolver::Vcycle
             scalarField& ACfRef =
                 const_cast<scalarField&>(ACf.operator const scalarField&());
 
-            if (interpolateCorrection_)
+            if (interpolateCorrection_ && leveli < coarsestLevel - 2)
             {
                 interpolate
                 (
@@ -322,7 +322,8 @@ void Foam::GAMGSolver::Vcycle
                     matrixLevels_[leveli],
                     interfaceLevelsBouCoeffs_[leveli],
                     interfaceLevels_[leveli],
-                    coarseSources[leveli],
+                    agglomeration_.restrictAddressing(leveli + 1),
+                    coarseCorrFields[leveli + 1],
                     cmpt
                 );
             }
@@ -382,7 +383,8 @@ void Foam::GAMGSolver::Vcycle
             matrix_,
             interfaceBouCoeffs_,
             interfaces_,
-            finestResidual,
+            agglomeration_.restrictAddressing(0),
+            coarseCorrFields[0],
             cmpt
         );
     }