GAMG: Add correction interpolation after prolongation

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverInterpolate.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"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::GAMGSolver::interpolate
+(
+    scalarField& psi,
+    scalarField& Apsi,
+    const lduMatrix& m,
+    const FieldField<Field, scalar>& interfaceBouCoeffs,
+    const lduInterfaceFieldPtrsList& interfaces,
+    const scalarField& source,
+    const direction cmpt
+) const
+{
+    scalar* __restrict__ psiPtr = psi.begin();
+
+    const label* const __restrict__ uPtr = m.lduAddr().upperAddr().begin();
+    const label* const __restrict__ lPtr = m.lduAddr().lowerAddr().begin();
+
+    const scalar* const __restrict__ diagPtr = m.diag().begin();
+    const scalar* const __restrict__ upperPtr = m.upper().begin();
+    const scalar* const __restrict__ lowerPtr = m.lower().begin();
+
+    Apsi = 0;
+    scalar* __restrict__ ApsiPtr = Apsi.begin();
+
+    m.initMatrixInterfaces
+    (
+        interfaceBouCoeffs,
+        interfaces,
+        psi,
+        Apsi,
+        cmpt
+    );
+
+    register const label nFaces = m.upper().size();
+    for (register label face=0; face<nFaces; face++)
+    {
+        ApsiPtr[uPtr[face]] += lowerPtr[face]*psiPtr[lPtr[face]];
+        ApsiPtr[lPtr[face]] += upperPtr[face]*psiPtr[uPtr[face]];
+    }
+
+    m.updateMatrixInterfaces
+    (
+        interfaceBouCoeffs,
+        interfaces,
+        psi,
+        Apsi,
+        cmpt
+    );
+
+    register const label nCells = m.diag().size();
+    for (register label celli=0; celli<nCells; celli++)
+    {
+        psiPtr[celli] = -ApsiPtr[celli]/(diagPtr[celli]);
+    }
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverScale.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  * * * * * * * * * * * //
+
+void Foam::GAMGSolver::scale
+(
+    scalarField& field,
+    scalarField& Acf,
+    const lduMatrix& A,
+    const FieldField<Field, scalar>& interfaceLevelBouCoeffs,
+    const lduInterfaceFieldPtrsList& interfaceLevel,
+    const scalarField& source,
+    const direction cmpt
+) const
+{
+    A.Amul
+    (
+        Acf,
+        field,
+        interfaceLevelBouCoeffs,
+        interfaceLevel,
+        cmpt
+    );
+
+    scalar scalingFactorNum = 0.0;
+    scalar scalingFactorDenom = 0.0;
+
+    forAll(field, i)
+    {
+        scalingFactorNum += source[i]*field[i];
+        scalingFactorDenom += Acf[i]*field[i];
+    }
+
+    vector2D scalingVector(scalingFactorNum, scalingFactorDenom);
+    reduce(scalingVector, sumOp<vector2D>());
+    scalar sf = scalingVector.x()/stabilise(scalingVector.y(), VSMALL);
+
+    if (debug >= 2)
+    {
+        Pout<< sf << " ";
+    }
+
+    const scalarField& D = A.diag();
+
+    forAll(field, i)
+    {
+        field[i] = sf*field[i] + (source[i] - sf*Acf[i])/D[i];
+    }
+}
+
+
+// ************************************************************************* //