ENH: lduMatrix: new matrix solvers: PPCG,PPCR

PPCG is pipelined version of PCG, PPCR is conjugate
residual version.

src/OpenFOAM/Make/files

@@ -383,6 +383,8 @@ $(lduMatrix)/solvers/smoothSolver/smoothSolver.C
 $(lduMatrix)/solvers/PCG/PCG.C
 $(lduMatrix)/solvers/PBiCG/PBiCG.C
 $(lduMatrix)/solvers/PBiCGStab/PBiCGStab.C
+$(lduMatrix)/solvers/PPCG/PPCG.C
+$(lduMatrix)/solvers/PPCR/PPCR.C
 
 $(lduMatrix)/smoothers/GaussSeidel/GaussSeidelSmoother.C
 $(lduMatrix)/smoothers/symGaussSeidel/symGaussSeidelSmoother.C

src/OpenFOAM/db/IOstreams/Pstreams/PstreamReduceOps.H

@@ -156,6 +156,21 @@ void reduce
     NotImplemented;
 }
 
+// Non-blocking version of reduce. Sets request.
+template<class T, class BinaryOp>
+void reduce
+(
+    T Value[],
+    const int size,
+    const BinaryOp& bop,
+    const int tag,
+    const label comm,
+    label& request
+)
+{
+    NotImplemented;
+}
+
 
 // Insist there are specialisations for the common reductions of scalar(s)
 void reduce
@@ -199,6 +214,16 @@ void reduce
     label& request
 );
 
+void reduce
+(
+    scalar Value[],
+    const int size,
+    const sumOp<scalar>& bop,
+    const int tag,
+    const label comm,
+    label& request
+);
+
 
 #if defined(WM_SPDP)
 void reduce
@@ -241,6 +266,16 @@ void reduce
     const label comm,
     label& request
 );
+
+void reduce
+(
+    solveScalar Value[],
+    const int size,
+    const sumOp<solveScalar>& bop,
+    const int tag,
+    const label comm,
+    label& request
+);
 #endif
 
 

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H

@@ -708,7 +708,8 @@ public:
                 const lduInterfaceFieldPtrsList& interfaces,
                 const solveScalarField& psiif,
                 solveScalarField& result,
-                const direction cmpt
+                const direction cmpt,
+                const label startRequest // starting request (for non-blocking)
             ) const;
 
             //- Set the residual field using an IOField on the object registry

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixATmul.C

@@ -56,6 +56,8 @@ void Foam::lduMatrix::Amul
     const scalar* const __restrict__ upperPtr = upper().begin();
     const scalar* const __restrict__ lowerPtr = lower().begin();
 
+    const label startRequest = Pstream::nRequests();
+
     // Initialise the update of interfaced interfaces
     initMatrixInterfaces
     (
@@ -90,7 +92,8 @@ void Foam::lduMatrix::Amul
         interfaces,
         psi,
         Apsi,
-        cmpt
+        cmpt,
+        startRequest
     );
 
     tpsi.clear();
@@ -119,6 +122,8 @@ void Foam::lduMatrix::Tmul
     const scalar* const __restrict__ lowerPtr = lower().begin();
     const scalar* const __restrict__ upperPtr = upper().begin();
 
+    const label startRequest = Pstream::nRequests();
+
     // Initialise the update of interfaced interfaces
     initMatrixInterfaces
     (
@@ -151,7 +156,8 @@ void Foam::lduMatrix::Tmul
         interfaces,
         psi,
         Tpsi,
-        cmpt
+        cmpt,
+        startRequest
     );
 
     tpsi.clear();
@@ -240,6 +246,8 @@ void Foam::lduMatrix::residual
     // To compensate for this, it is necessary to turn the
     // sign of the contribution.
 
+    const label startRequest = Pstream::nRequests();
+
     // Initialise the update of interfaced interfaces
     initMatrixInterfaces
     (
@@ -274,7 +282,8 @@ void Foam::lduMatrix::residual
         interfaces,
         psi,
         rA,
-        cmpt
+        cmpt,
+        startRequest
     );
 }
 

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C

@@ -106,7 +106,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
     const lduInterfaceFieldPtrsList& interfaces,
     const solveScalarField& psiif,
     solveScalarField& result,
-    const direction cmpt
+    const direction cmpt,
+    const label startRequest
 ) const
 {
     if (Pstream::defaultCommsType == Pstream::commsTypes::blocking)
@@ -183,7 +184,7 @@ void Foam::lduMatrix::updateMatrixInterfaces
             else
             {
                 // Block for all requests and remove storage
-                UPstream::waitRequests();
+                UPstream::waitRequests(startRequest);
             }
         }
 

src/OpenFOAM/matrices/lduMatrix/smoothers/GaussSeidel/GaussSeidelSmoother.C

@@ -115,6 +115,8 @@ void Foam::GaussSeidelSmoother::smooth
     {
         bPrime = source;
 
+        const label startRequest = Pstream::nRequests();
+
         matrix_.initMatrixInterfaces
         (
             false,
@@ -132,7 +134,8 @@ void Foam::GaussSeidelSmoother::smooth
             interfaces_,
             psi,
             bPrime,
-            cmpt
+            cmpt,
+            startRequest
         );
 
         solveScalar psii;

src/OpenFOAM/matrices/lduMatrix/smoothers/nonBlockingGaussSeidel/nonBlockingGaussSeidelSmoother.C

@@ -142,6 +142,8 @@ void Foam::nonBlockingGaussSeidelSmoother::smooth
     {
         bPrime = source;
 
+        const label startRequest = Pstream::nRequests();
+
         matrix_.initMatrixInterfaces
         (
             false,
@@ -190,7 +192,8 @@ void Foam::nonBlockingGaussSeidelSmoother::smooth
             interfaces_,
             psi,
             bPrime,
-            cmpt
+            cmpt,
+            startRequest
         );
 
         // Update rest of the cells

src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C

@@ -115,6 +115,8 @@ void Foam::symGaussSeidelSmoother::smooth
     {
         bPrime = source;
 
+        const label startRequest = Pstream::nRequests();
+
         matrix_.initMatrixInterfaces
         (
             false,
@@ -132,7 +134,8 @@ void Foam::symGaussSeidelSmoother::smooth
             interfaces_,
             psi,
             bPrime,
-            cmpt
+            cmpt,
+            startRequest
         );
 
         solveScalar psii;

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

@@ -52,6 +52,8 @@ void Foam::GAMGSolver::interpolate
     Apsi = 0;
     solveScalar* __restrict__ ApsiPtr = Apsi.begin();
 
+    const label startRequest = Pstream::nRequests();
+
     m.initMatrixInterfaces
     (
         true,
@@ -76,7 +78,8 @@ void Foam::GAMGSolver::interpolate
         interfaces,
         psi,
         Apsi,
-        cmpt
+        cmpt,
+        startRequest
     );
 
     const label nCells = m.diag().size();

src/OpenFOAM/matrices/lduMatrix/solvers/PPCG/PPCG.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019-2020 M. Janssens
+-------------------------------------------------------------------------------
+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 "PPCG.H"
+#include "PrecisionAdaptor.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(PPCG, 0);
+
+    lduMatrix::solver::addsymMatrixConstructorToTable<PPCG>
+        addPPCGSymMatrixConstructorToTable_;
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::PPCG::gSumMagProd
+(
+    FixedList<solveScalar, 3>& globalSum,
+    const solveScalarField& a,
+    const solveScalarField& b,
+    const solveScalarField& c,
+    const solveScalarField& sumMag,
+    label& outstandingRequest,
+    const label comm
+) const
+{
+    const label nCells = a.size();
+
+    globalSum = 0.0;
+    for (label cell=0; cell<nCells; cell++)
+    {
+        globalSum[0] += a[cell]*b[cell];    // sumProd(a, b)
+        globalSum[1] += a[cell]*c[cell];    // sumProd(a, c)
+        globalSum[2] += mag(sumMag[cell]);
+    }
+
+    if (Pstream::parRun())
+    {
+        Foam::reduce
+        (
+            globalSum.begin(),
+            globalSum.size(),
+            sumOp<solveScalar>(),
+            Pstream::msgType(),
+            comm,
+            outstandingRequest
+        );
+    }
+}
+
+
+Foam::solverPerformance Foam::PPCG::scalarSolve
+(
+    solveScalarField& psi,
+    const solveScalarField& source,
+    const direction cmpt,
+    const bool cgMode
+) const
+{
+    // --- Setup class containing solver performance data
+    solverPerformance solverPerf
+    (
+        lduMatrix::preconditioner::getName(controlDict_) + type(),
+        fieldName_
+    );
+
+    const label comm = matrix().mesh().comm();
+    const label nCells = psi.size();
+    solveScalarField w(nCells);
+
+    // --- Calculate A.psi
+    matrix_.Amul(w, psi, interfaceBouCoeffs_, interfaces_, cmpt);
+
+    // --- Calculate initial residual field
+    solveScalarField r(source - w);
+
+    // --- Calculate normalisation factor
+    solveScalarField p(nCells);
+    const solveScalar normFactor = this->normFactor(psi, source, w, p);
+
+    if (lduMatrix::debug >= 2)
+    {
+        Info<< "   Normalisation factor = " << normFactor << endl;
+    }
+
+    // --- Select and construct the preconditioner
+    autoPtr<lduMatrix::preconditioner> preconPtr =
+    lduMatrix::preconditioner::New
+    (
+        *this,
+        controlDict_
+    );
+
+    // --- Precondition residual (= u0)
+    solveScalarField u(nCells);
+    preconPtr->precondition(u, r, cmpt);
+
+    // --- Calculate A*u - reuse w
+    matrix_.Amul(w, u, interfaceBouCoeffs_, interfaces_, cmpt);
+
+
+    // State
+    solveScalarField s(nCells);
+    solveScalarField q(nCells);
+    solveScalarField z(nCells);
+
+    solveScalarField m(nCells);
+
+    FixedList<solveScalar, 3> globalSum;
+    label outstandingRequest = -1;
+    if (cgMode)
+    {
+        // --- Start global reductions for inner products
+        gSumMagProd(globalSum, u, r, w, r, outstandingRequest, comm);
+
+        // --- Precondition residual
+        preconPtr->precondition(m, w, cmpt);
+    }
+    else
+    {
+        // --- Precondition residual
+        preconPtr->precondition(m, w, cmpt);
+
+        // --- Start global reductions for inner products
+        gSumMagProd(globalSum, w, u, m, r, outstandingRequest, comm);
+    }
+
+    // --- Calculate A*m
+    solveScalarField n(nCells);
+    matrix_.Amul(n, m, interfaceBouCoeffs_, interfaces_, cmpt);
+
+    solveScalar alpha = 0.0;
+    solveScalar gamma = 0.0;
+
+    // --- Solver iteration
+    for
+    (
+        solverPerf.nIterations() = 0;
+        solverPerf.nIterations() < maxIter_;
+        solverPerf.nIterations()++
+    )
+    {
+        // Make sure gamma,delta are available
+        if (Pstream::parRun())
+        {
+            Pstream::waitRequest(outstandingRequest);
+            outstandingRequest = -1;
+        }
+
+        const solveScalar gammaOld = gamma;
+        gamma = globalSum[0];
+        const solveScalar delta = globalSum[1];
+
+        solverPerf.finalResidual() = globalSum[2]/normFactor;
+        if (solverPerf.nIterations() == 0)
+        {
+            solverPerf.initialResidual() = solverPerf.finalResidual();
+        }
+
+        // Check convergence (bypass if not enough iterations yet)
+        if
+        (
+            (minIter_ <= 0 || solverPerf.nIterations() >= minIter_)
+         && solverPerf.checkConvergence(tolerance_, relTol_)
+        )
+        {
+            break;
+        }
+
+
+        if (solverPerf.nIterations() == 0)
+        {
+            alpha = gamma/delta;
+            z = n;
+            q = m;
+            s = w;
+            p = u;
+        }
+        else
+        {
+            const solveScalar beta = gamma/gammaOld;
+            alpha = gamma/(delta-beta*gamma/alpha);
+
+            for (label cell=0; cell<nCells; cell++)
+            {
+                z[cell] = n[cell] + beta*z[cell];
+                q[cell] = m[cell] + beta*q[cell];
+                s[cell] = w[cell] + beta*s[cell];
+                p[cell] = u[cell] + beta*p[cell];
+            }
+        }
+
+        for (label cell=0; cell<nCells; cell++)
+        {
+            psi[cell] += alpha*p[cell];
+            r[cell] -= alpha*s[cell];
+            u[cell] -= alpha*q[cell];
+            w[cell] -= alpha*z[cell];
+        }
+
+        if (cgMode)
+        {
+            // --- Start global reductions for inner products
+            gSumMagProd(globalSum, u, r, w, r, outstandingRequest, comm);
+
+            // --- Precondition residual
+            preconPtr->precondition(m, w, cmpt);
+        }
+        else
+        {
+            // --- Precondition residual
+            preconPtr->precondition(m, w, cmpt);
+
+            // --- Start global reductions for inner products
+            gSumMagProd(globalSum, w, u, m, r, outstandingRequest, comm);
+        }
+
+        // --- Calculate A*m
+        matrix_.Amul(n, m, interfaceBouCoeffs_, interfaces_, cmpt);
+    }
+
+    return solverPerf;
+}
+
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::PPCG::PPCG
+(
+    const word& fieldName,
+    const lduMatrix& matrix,
+    const FieldField<Field, scalar>& interfaceBouCoeffs,
+    const FieldField<Field, scalar>& interfaceIntCoeffs,
+    const lduInterfaceFieldPtrsList& interfaces,
+    const dictionary& solverControls
+)
+:
+    lduMatrix::solver
+    (
+        fieldName,
+        matrix,
+        interfaceBouCoeffs,
+        interfaceIntCoeffs,
+        interfaces,
+        solverControls
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::solverPerformance Foam::PPCG::solve
+(
+    scalarField& psi_s,
+    const scalarField& source,
+    const direction cmpt
+) const
+{
+    PrecisionAdaptor<solveScalar, scalar> tpsi(psi_s);
+    return scalarSolve
+    (
+        tpsi.ref(),
+        ConstPrecisionAdaptor<solveScalar, scalar>(source)(),
+        cmpt,
+        true   // operate in conjugate-gradient mode
+    );
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/PPCG/PPCG.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019-2020 M. Janssens
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::PPCG
+
+Description
+    Preconditioned pipelined conjugate gradient solver for symmetric
+    lduMatrices using a run-time selectable preconditioner.
+
+    Reference:
+    \verbatim
+        P. Ghysels, W. Vanroose.
+        "Hiding global synchronization latency in the
+         preconditioned Conjugate Gradient algorithm"
+    \endverbatim
+    and implementation details from
+    \verbatim
+        Paul Eller, William Gropp
+        "Scalable Non-blocking Preconditioned Conjugate Gradient Methods"
+    \endverbatim
+
+SourceFiles
+    PPCG.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PPCG_H
+#define PPCG_H
+
+#include "lduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class PPCG Declaration
+\*---------------------------------------------------------------------------*/
+
+class PPCG
+:
+    public lduMatrix::solver
+{
+    // Private Member Functions
+
+        //- Non-blocking version of sum(a*b), sum(a*c), sum(mag(sumMag))
+        void gSumMagProd
+        (
+            FixedList<solveScalar, 3>& globalSum,
+            const solveScalarField& a,
+            const solveScalarField& b,
+            const solveScalarField& c,
+            const solveScalarField& sumMag,
+            label& outstandingRequest,
+