From 4f7aa70ad6b1bfdf85262f33f561d1b23c751fd2 Mon Sep 17 00:00:00 2001 From: Yu Ankun Date: Thu, 20 Apr 2023 13:56:38 +0300 Subject: [PATCH] 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 | 1 + .../matrices/lduMatrix/solvers/FPCG/FPCG.C | 269 ++++++++++++++++++ .../matrices/lduMatrix/solvers/FPCG/FPCG.H | 134 +++++++++ .../matrices/lduMatrix/solvers/PCG/PCG.C | 4 +- 4 files changed, 406 insertions(+), 2 deletions(-) create mode 100644 src/OpenFOAM/matrices/lduMatrix/solvers/FPCG/FPCG.C create mode 100644 src/OpenFOAM/matrices/lduMatrix/solvers/FPCG/FPCG.H diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index 7753533552..9eb84d9243 100644 --- a/src/OpenFOAM/Make/files +++ b/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 diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/FPCG/FPCG.C b/src/OpenFOAM/matrices/lduMatrix/solvers/FPCG/FPCG.C new file mode 100644 index 0000000000..36bcdc99d0 --- /dev/null +++ b/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 . + +\*---------------------------------------------------------------------------*/ + +#include "FPCG.H" +#include "PrecisionAdaptor.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(FPCG, 0); + + lduMatrix::solver::addsymMatrixConstructorToTable + addFPCGSymMatrixConstructorToTable_; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::FPCG::FPCG +( + const word& fieldName, + const lduMatrix& matrix, + const FieldField& interfaceBouCoeffs, + const FieldField& interfaceIntCoeffs, + const lduInterfaceFieldPtrsList& interfaces, + const dictionary& solverControls +) +: + lduMatrix::solver + ( + fieldName, + matrix, + interfaceBouCoeffs, + interfaceIntCoeffs, + interfaces, + solverControls + ) +{} + + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +void Foam::FPCG::gSumMagProd +( + FixedList& globalSum, + const solveScalarField& a, + const solveScalarField& b, + const label comm +) +{ + const label nCells = a.size(); + + globalSum = 0.0; + for (label cell=0; cell(), + 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(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 preconPtr = + lduMatrix::preconditioner::New + ( + *this, + controlDict_ + ); + + FixedList 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(rA)(), + fieldName_, + false + ); + + return solverPerf; +} + + + +Foam::solverPerformance Foam::FPCG::solve +( + scalarField& psi_s, + const scalarField& source, + const direction cmpt +) const +{ + PrecisionAdaptor tpsi(psi_s); + return scalarSolve + ( + tpsi.ref(), + ConstPrecisionAdaptor(source)(), + cmpt + ); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/FPCG/FPCG.H b/src/OpenFOAM/matrices/lduMatrix/solvers/FPCG/FPCG.H new file mode 100644 index 0000000000..16daa7b795 --- /dev/null +++ b/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 . + +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& 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& interfaceBouCoeffs, + const FieldField& 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 + +// ************************************************************************* // diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C b/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C index ce33493e17..98023fbf83 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/PCG/PCG.C +++ b/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