GAMG: Add optional correction interpolate after prolongation

src/OpenFOAM/Make/files

@@ -284,7 +284,8 @@ $(lduInterfaceFields)/cyclicLduInterfaceField/cyclicLduInterfaceField.C
 GAMG = $(lduMatrix)/solvers/GAMG
 $(GAMG)/GAMGSolver.C
 $(GAMG)/GAMGSolverAgglomerateMatrix.C
-$(GAMG)/GAMGSolverScalingFactor.C
+$(GAMG)/GAMGSolverInterpolate.C
+$(GAMG)/GAMGSolverScale.C
 $(GAMG)/GAMGSolverSolve.C
 
 GAMGInterfaces = $(GAMG)/interfaces

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

@@ -71,6 +71,7 @@ Foam::GAMGSolver::GAMGSolver
     postSweepsLevelMultiplier_(1),
     maxPostSweeps_(4),
     nFinestSweeps_(2),
+    interpolateCorrection_(false),
     scaleCorrection_(matrix.symmetric()),
     directSolveCoarsest_(false),
     agglomeration_(GAMGAgglomeration::New(matrix_, controlDict_)),
@@ -174,6 +175,7 @@ void Foam::GAMGSolver::readControls()
     );
     controlDict_.readIfPresent("maxPostSweeps", maxPostSweeps_);
     controlDict_.readIfPresent("nFinestSweeps", nFinestSweeps_);
+    controlDict_.readIfPresent("interpolateCorrection", interpolateCorrection_);
     controlDict_.readIfPresent("scaleCorrection", scaleCorrection_);
     controlDict_.readIfPresent("directSolveCoarsest", directSolveCoarsest_);
 }

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

@@ -46,6 +46,8 @@ SourceFiles
     GAMGSolverCalcAgglomeration.C
     GAMGSolverMakeCoarseMatrix.C
     GAMGSolverOperations.C
+    GAMGSolverInterpolate.C
+    GAMGSolverScale.C
     GAMGSolverSolve.C
 
 \*---------------------------------------------------------------------------*/
@@ -97,6 +99,10 @@ class GAMGSolver
         //- Number of smoothing sweeps on finest mesh
         label nFinestSweeps_;
 
+        //- Choose if the corrections should be interpolated after injection.
+        //  By default corrections are not interpolated.
+        bool interpolateCorrection_;
+
         //- Choose if the corrections should be scaled.
         //  By default corrections for symmetric matrices are scaled
         //  but not for asymmetric matrices.
@@ -154,32 +160,34 @@ class GAMGSolver
         //- Agglomerate coarse matrix
         void agglomerateMatrix(const label fineLevelIndex);
 
-        //- Calculate and return the scaling factor from Acf, coarseSource
+        //-  Interpolate the correction after injected prolongation
+        void interpolate
+        (
+            scalarField& psi,
+            scalarField& Apsi,
+            const lduMatrix& m,
+            const FieldField<Field, scalar>& interfaceBouCoeffs,
+            const lduInterfaceFieldPtrsList& interfaces,
+            const scalarField& source,
+            const direction cmpt
+        ) const;
+
+        //- Calculate and apply the scaling factor from Acf, coarseSource
         //  and coarseField.
         //  At the same time do a Jacobi iteration on the coarseField using
         //  the Acf provided after the coarseField values are used for the
         //  scaling factor.
-        scalar scalingFactor
+        void scale
         (
             scalarField& field,
-            const scalarField& source,
-            const scalarField& Acf,
-            const scalarField& D
-        ) const;
-
-        //- Calculate Acf and calculate and return the scaling factor.
-        scalar scalingFactor
-        (
             scalarField& Acf,
             const lduMatrix& A,
-            scalarField& field,
             const FieldField<Field, scalar>& interfaceLevelBouCoeffs,
             const lduInterfaceFieldPtrsList& interfaceLevel,
             const scalarField& source,
             const direction cmpt
         ) const;
 
-
         //- Initialise the data structures for the V-cycle
         void initVcycle
         (

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

-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-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 "GAMGSolver.H"
-#include "vector2D.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::scalar Foam::GAMGSolver::scalingFactor
-(
-    scalarField& field,
-    const scalarField& source,
-    const scalarField& Acf,
-    const scalarField& D
-) const
-{
-    scalar scalingFactorNum = 0.0;
-    scalar scalingFactorDenom = 0.0;
-
-    forAll(field, i)
-    {
-        scalingFactorNum += source[i]*field[i];
-        scalingFactorDenom += Acf[i]*field[i];
-
-        // While the matrix-multiply done for the scaling it is
-        // possible to perform a point-Jacobi smoothing operation cheaply
-        field[i] += (source[i] - Acf[i])/D[i];
-    }
-
-    vector2D scalingVector(scalingFactorNum, scalingFactorDenom);
-    reduce(scalingVector, sumOp<vector2D>());
-    return scalingVector.x()/stabilise(scalingVector.y(), VSMALL);
-}
-
-
-Foam::scalar Foam::GAMGSolver::scalingFactor
-(
-    scalarField& Acf,
-    const lduMatrix& A,
-    scalarField& field,
-    const FieldField<Field, scalar>& interfaceLevelBouCoeffs,
-    const lduInterfaceFieldPtrsList& interfaceLevel,
-    const scalarField& source,
-    const direction cmpt
-) const
-{
-    A.Amul
-    (
-        Acf,
-        field,
-        interfaceLevelBouCoeffs,
-        interfaceLevel,
-        cmpt
-    );
-
-    return scalingFactor
-    (
-        field,
-        source,
-        Acf,
-        A.diag()
-    );
-}
-
-
-// ************************************************************************* //

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

@@ -28,6 +28,7 @@ License
 #include "BICCG.H"
 #include "SubField.H"
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::solverPerformance Foam::GAMGSolver::solve
@@ -173,23 +174,16 @@ void Foam::GAMGSolver::Vcycle
             // but not on the coarsest level because it evaluates to 1
             if (scaleCorrection_ && leveli < coarsestLevel - 1)
             {
-                scalar sf = scalingFactor
+                scale
                 (
+                    coarseCorrFields[leveli],
                     const_cast<scalarField&>(ACf.operator const scalarField&()),
                     matrixLevels_[leveli],
-                    coarseCorrFields[leveli],
                     interfaceLevelsBouCoeffs_[leveli],
                     interfaceLevels_[leveli],
                     coarseSources[leveli],
                     cmpt
                 );
-
-                if (debug >= 2)
-                {
-                    Pout<< sf << " ";
-                }
-
-                coarseCorrFields[leveli] *= sf;
             }
 
             // Correct the residual with the new solution
@@ -246,7 +240,7 @@ void Foam::GAMGSolver::Vcycle
             coarseCorrFields[leveli].size()
         );
 
-        // Only store the preSmoothedCoarseCorrField is pre-smoothing is used
+        // Only store the preSmoothedCoarseCorrField if pre-smoothing is used
         if (nPreSweeps_)
         {
             preSmoothedCoarseCorrField.assign(coarseCorrFields[leveli]);
@@ -259,38 +253,47 @@ void Foam::GAMGSolver::Vcycle
             leveli + 1
         );
 
-        // Scale coarse-grid correction field
-        // but not on the coarsest level because it evaluates to 1
-        if (scaleCorrection_ && leveli < coarsestLevel - 1)
+        // Create A.psi for this coarse level as a sub-field of Apsi
+        scalarField::subField ACf
+        (
+            Apsi,
+            coarseCorrFields[leveli].size()
+        );
+
+        scalarField& ACfRef =
+            const_cast<scalarField&>(ACf.operator const scalarField&());
+
+        if (interpolateCorrection_)
         {
-            // Create A.psi for this coarse level as a sub-field of Apsi
-            scalarField::subField ACf
+            interpolate
             (
-                Apsi,
-                coarseCorrFields[leveli].size()
+                coarseCorrFields[leveli],
+                ACfRef,
+                matrixLevels_[leveli],
+                interfaceLevelsBouCoeffs_[leveli],
+                interfaceLevels_[leveli],
+                coarseSources[leveli],
+                cmpt
             );
+        }
 
-            scalar sf = scalingFactor
+        // Scale coarse-grid correction field
+        // but not on the coarsest level because it evaluates to 1
+        if (scaleCorrection_ && leveli < coarsestLevel - 1)
+        {
+            scale
             (
-                const_cast<scalarField&>(ACf.operator const scalarField&()),
-                matrixLevels_[leveli],
                 coarseCorrFields[leveli],
+                ACfRef,
+                matrixLevels_[leveli],
                 interfaceLevelsBouCoeffs_[leveli],
                 interfaceLevels_[leveli],
                 coarseSources[leveli],
                 cmpt
             );
-
-
-            if (debug >= 2)
-            {
-                Pout<< sf << " ";
-            }
-
-            coarseCorrFields[leveli] *= sf;
         }
 
-        // Only add the preSmoothedCoarseCorrField is pre-smoothing is used
+        // Only add the preSmoothedCoarseCorrField if pre-smoothing is used
         if (nPreSweeps_)
         {
             coarseCorrFields[leveli] += preSmoothedCoarseCorrField;
@@ -317,36 +320,38 @@ void Foam::GAMGSolver::Vcycle
         0
     );
 
-    if (scaleCorrection_)
+    if (interpolateCorrection_)
     {
-        // Calculate finest level scaling factor
-        scalar fsf = scalingFactor
+        interpolate
         (
+            finestCorrection,
             Apsi,
             matrix_,
-            finestCorrection,
             interfaceBouCoeffs_,
             interfaces_,
             finestResidual,
             cmpt
         );
+    }
 
-        if (debug >= 2)
-        {
-            Pout<< fsf << endl;
-        }
-
-        forAll(psi, i)
-        {
-            psi[i] += fsf*finestCorrection[i];
-        }
+    if (scaleCorrection_)
+    {
+        // Scale the finest level correction
+        scale
+        (
+            finestCorrection,
+            Apsi,
+            matrix_,
+            interfaceBouCoeffs_,
+            interfaces_,
+            finestResidual,
+            cmpt
+        );
     }
-    else
+
+    forAll(psi, i)
     {
-        forAll(psi, i)
-        {
-            psi[i] += finestCorrection[i];
-        }
+        psi[i] += finestCorrection[i];
     }
 
     smoothers[0].smooth