ENH: new FPCG solver (faster PCG)

- this should be faster than the regular PCG on larger systems
  since it combines two global reductions

src/OpenFOAM/Make/files

@@ -431,6 +431,7 @@ $(lduMatrix)/solvers/smoothSolver/smoothSolver.C
 $(lduMatrix)/solvers/PCG/PCG.C
 $(lduMatrix)/solvers/PBiCG/PBiCG.C
 $(lduMatrix)/solvers/PBiCGStab/PBiCGStab.C
+$(lduMatrix)/solvers/FPCG/FPCG.C
 $(lduMatrix)/solvers/PPCG/PPCG.C
 $(lduMatrix)/solvers/PPCR/PPCR.C
 

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

@@ -0,0 +1,269 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2023 Huawei (Yu Ankun)
+-------------------------------------------------------------------------------
+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 "FPCG.H"
+#include "PrecisionAdaptor.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(FPCG, 0);
+
+    lduMatrix::solver::addsymMatrixConstructorToTable<FPCG>
+        addFPCGSymMatrixConstructorToTable_;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::FPCG::FPCG
+(
+    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
+    )
+{}
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+void Foam::FPCG::gSumMagProd
+(
+    FixedList<solveScalar, 2>& globalSum,
+    const solveScalarField& a,
+    const solveScalarField& b,
+    const label comm
+)
+{
+    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] += mag(b[cell]);       // sumMag(b)
+    }
+
+    Foam::reduce
+    (
+        globalSum.data(),
+        globalSum.size(),
+        sumOp<solveScalar>(),
+        UPstream::msgType(),  // (ignored): direct MPI call
+        comm
+    );
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::solverPerformance Foam::FPCG::scalarSolve
+(
+    solveScalarField& psi,
+    const solveScalarField& source,
+    const direction cmpt
+) const
+{
+    // --- Setup class containing solver performance data
+    solverPerformance solverPerf
+    (
+        lduMatrix::preconditioner::getName(controlDict_) + typeName,
+        fieldName_
+    );
+
+    label nCells = psi.size();
+
+    solveScalar* __restrict__ psiPtr = psi.begin();
+
+    solveScalarField pA(nCells);
+    solveScalar* __restrict__ pAPtr = pA.begin();
+
+    solveScalarField wA(nCells);
+    solveScalar* __restrict__ wAPtr = wA.begin();
+
+    solveScalar wArA = solverPerf.great_;
+    solveScalar wArAold = wArA;
+
+    // --- Calculate A.psi
+    matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
+
+    // --- Calculate initial residual field
+    solveScalarField rA(source - wA);
+    solveScalar* __restrict__ rAPtr = rA.begin();
+
+    matrix().setResidualField
+    (
+        ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),
+        fieldName_,
+        true
+    );
+
+    // --- Calculate normalisation factor
+    solveScalar normFactor = this->normFactor(psi, source, wA, pA);
+
+    if ((log_ >= 2) || (lduMatrix::debug >= 2))
+    {
+        Info<< "   Normalisation factor = " << normFactor << endl;
+    }
+
+    // --- Calculate normalised residual norm
+    solverPerf.initialResidual() =
+        gSumMag(rA, matrix().mesh().comm())
+       /normFactor;
+    solverPerf.finalResidual() = solverPerf.initialResidual();
+
+    // --- Check convergence, solve if not converged
+    if
+    (
+        minIter_ > 0
+     || !solverPerf.checkConvergence(tolerance_, relTol_, log_)
+    )
+    {
+        // --- Select and construct the preconditioner
+        autoPtr<lduMatrix::preconditioner> preconPtr =
+            lduMatrix::preconditioner::New
+            (
+                *this,
+                controlDict_
+            );
+
+        FixedList<solveScalar, 2> globalSum;
+
+        // --- Solver iteration
+        do
+        {
+            // --- Store previous wArA
+            wArAold = wArA;
+
+            // --- Precondition residual
+            preconPtr->precondition(wA, rA, cmpt);
+
+            // --- Update search directions and calculate residual:
+            gSumMagProd(globalSum, wA, rA, matrix().mesh().comm());
+
+            wArA = globalSum[0];
+
+            solverPerf.finalResidual() = globalSum[1]/normFactor;
+
+            // Check convergence (bypass if not enough iterations yet)
+            if
+            (
+                (minIter_ <= 0 || solverPerf.nIterations() >= minIter_)
+             && solverPerf.checkConvergence(tolerance_, relTol_, log_)
+            )
+            {
+                break;
+            }
+
+            if (solverPerf.nIterations() == 0)
+            {
+                for (label cell=0; cell<nCells; cell++)
+                {
+                    pAPtr[cell] = wAPtr[cell];
+                }
+            }
+            else
+            {
+                const solveScalar beta = wArA/wArAold;
+
+                for (label cell=0; cell<nCells; cell++)
+                {
+                    pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
+                }
+            }
+
+
+            // --- Update preconditioned residual
+            matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);
+
+            solveScalar wApA = gSumProd(wA, pA, matrix().mesh().comm());
+
+            // --- Test for singularity
+            if (solverPerf.checkSingularity(mag(wApA)/normFactor)) break;
+
+
+            // --- Update solution and residual:
+
+            const solveScalar alpha = wArA/wApA;
+
+            for (label cell=0; cell<nCells; cell++)
+            {
+                psiPtr[cell] += alpha*pAPtr[cell];
+                rAPtr[cell] -= alpha*wAPtr[cell];
+            }
+
+        } while
+        (
+            ++solverPerf.nIterations() < maxIter_
+        );
+    }
+
+    matrix().setResidualField
+    (
+        ConstPrecisionAdaptor<scalar, solveScalar>(rA)(),
+        fieldName_,
+        false
+    );
+
+    return solverPerf;
+}
+
+
+
+Foam::solverPerformance Foam::FPCG::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
+    );
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2012 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2023 Huawei (Yu Ankun)
+-------------------------------------------------------------------------------
+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::FPCG
+
+Group
+    grpLduMatrixSolvers
+
+Description
+    A 'faster' preconditioned conjugate gradient solver for
+    symmetric lduMatrices using a run-time selectable preconditioner.
+
+    This is termed \em "faster" than the regular PCG since it combines
+    global reductions.
+
+SourceFiles
+    FPCG.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_FPCG_H
+#define Foam_FPCG_H
+
+#include "lduMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                            Class FPCG Declaration
+\*---------------------------------------------------------------------------*/
+
+class FPCG
+:
+    public lduMatrix::solver
+{
+    // Private Member Functions
+
+        //- Blocking version of sum(a*b), sum(mag(b))
+        static void gSumMagProd
+        (
+            FixedList<solveScalar, 2>& globalSum,
+            const solveScalarField& a,
+            const solveScalarField& b,
+            const label comm
+        );
+
+        //- No copy construct
+        FPCG(const FPCG&) = delete;
+
+        //- No copy assignment
+        void operator=(const FPCG&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("FPCG");
+
+
+    // Constructors
+
+        //- Construct from matrix components and solver controls
+        FPCG
+        (
+            const word& fieldName,
+            const lduMatrix& matrix,
+            const FieldField<Field, scalar>& interfaceBouCoeffs,
+            const FieldField<Field, scalar>& interfaceIntCoeffs,
+            const lduInterfaceFieldPtrsList& interfaces,
+            const dictionary& solverControls
+        );
+
+
+    //- Destructor
+    virtual ~FPCG() = default;
+
+
+    // Member Functions
+
+        //- Solve the matrix with this solver
+        virtual solverPerformance scalarSolve
+        (
+            solveScalarField& psi,
+            const solveScalarField& source,
+            const direction cmpt=0
+        ) const;
+
+        //- Solve the matrix with this solver
+        virtual solverPerformance solve
+        (
+            scalarField& psi,
+            const scalarField& source,
+            const direction cmpt=0
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -157,7 +157,7 @@ Foam::solverPerformance Foam::PCG::scalarSolve
             }
             else
             {
-                solveScalar beta = wArA/wArAold;
+                const solveScalar beta = wArA/wArAold;
 
                 for (label cell=0; cell<nCells; cell++)
                 {
@@ -177,7 +177,7 @@ Foam::solverPerformance Foam::PCG::scalarSolve
 
             // --- Update solution and residual:
 
-            solveScalar alpha = wArA/wApA;
+            const solveScalar alpha = wArA/wApA;
 
             for (label cell=0; cell<nCells; cell++)
             {