Commit 0857f479 authored by Henry Weller's avatar Henry Weller
Browse files

PBiCGStab: New preconditioned bi-conjugate gradient stabilized solver for asymmetric lduMatrices

using a run-time selectable preconditioner

References:
    Van der Vorst, H. A. (1992).
    Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG
    for the solution of nonsymmetric linear systems.
    SIAM Journal on scientific and Statistical Computing, 13(2), 631-644.

    Barrett, R., Berry, M. W., Chan, T. F., Demmel, J., Donato, J.,
    Dongarra, J., Eijkhout, V., Pozo, R., Romine, C. & Van der Vorst, H.
    (1994).
    Templates for the solution of linear systems:
    building blocks for iterative methods
    (Vol. 43). Siam.

See also: https://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method

Tests have shown that PBiCGStab with the DILU preconditioner is more
robust, reliable and shows faster convergence (~2x) than PBiCG with
DILU, in particular in parallel where PBiCG occasionally diverges.

This remarkable improvement over PBiCG prompted the update of all
tutorial cases currently using PBiCG to use PBiCGStab instead.  If any
issues arise with this update please report on Mantis: http://bugs.openfoam.org
parent b19fa431
......@@ -271,6 +271,7 @@ $(lduMatrix)/solvers/diagonalSolver/diagonalSolver.C
$(lduMatrix)/solvers/smoothSolver/smoothSolver.C
$(lduMatrix)/solvers/PCG/PCG.C
$(lduMatrix)/solvers/PBiCG/PBiCG.C
$(lduMatrix)/solvers/PBiCGStab/PBiCGStab.C
$(lduMatrix)/smoothers/GaussSeidel/GaussSeidelSmoother.C
$(lduMatrix)/smoothers/symGaussSeidel/symGaussSeidelSmoother.C
......
......@@ -76,37 +76,25 @@ Foam::solverPerformance Foam::PBiCG::solve
fieldName_
);
label nCells = psi.size();
const label nCells = psi.size();
scalar* __restrict__ psiPtr = psi.begin();
scalarField pA(nCells);
scalar* __restrict__ pAPtr = pA.begin();
scalarField pT(nCells, 0.0);
scalar* __restrict__ pTPtr = pT.begin();
scalarField wA(nCells);
scalar* __restrict__ wAPtr = wA.begin();
scalarField wT(nCells);
scalar* __restrict__ wTPtr = wT.begin();
scalar wArT = solverPerf.great_;
scalar wArTold = wArT;
// --- Calculate A.psi and T.psi
// --- Calculate A.psi
matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
matrix_.Tmul(wT, psi, interfaceIntCoeffs_, interfaces_, cmpt);
// --- Calculate initial residual and transpose residual fields
// --- Calculate initial residual field
scalarField rA(source - wA);
scalarField rT(source - wT);
scalar* __restrict__ rAPtr = rA.begin();
scalar* __restrict__ rTPtr = rT.begin();
// --- Calculate normalisation factor
scalar normFactor = this->normFactor(psi, source, wA, pA);
const scalar normFactor = this->normFactor(psi, source, wA, pA);
if (lduMatrix::debug >= 2)
{
......@@ -126,6 +114,22 @@ Foam::solverPerformance Foam::PBiCG::solve
|| !solverPerf.checkConvergence(tolerance_, relTol_)
)
{
scalarField pT(nCells, 0);
scalar* __restrict__ pTPtr = pT.begin();
scalarField wT(nCells);
scalar* __restrict__ wTPtr = wT.begin();
// --- Calculate T.psi
matrix_.Tmul(wT, psi, interfaceIntCoeffs_, interfaces_, cmpt);
// --- Calculate initial transpose residual field
scalarField rT(source - wT);
scalar* __restrict__ rTPtr = rT.begin();
// --- Initial value not used
scalar wArT = 0;
// --- Select and construct the preconditioner
autoPtr<lduMatrix::preconditioner> preconPtr =
lduMatrix::preconditioner::New
......@@ -138,7 +142,7 @@ Foam::solverPerformance Foam::PBiCG::solve
do
{
// --- Store previous wArT
wArTold = wArT;
const scalar wArTold = wArT;
// --- Precondition residuals
preconPtr->precondition(wA, rA, cmpt);
......@@ -157,7 +161,7 @@ Foam::solverPerformance Foam::PBiCG::solve
}
else
{
scalar beta = wArT/wArTold;
const scalar beta = wArT/wArTold;
for (label cell=0; cell<nCells; cell++)
{
......@@ -171,7 +175,7 @@ Foam::solverPerformance Foam::PBiCG::solve
matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);
matrix_.Tmul(wT, pT, interfaceIntCoeffs_, interfaces_, cmpt);
scalar wApT = gSumProd(wA, pT, matrix().mesh().comm());
const scalar wApT = gSumProd(wA, pT, matrix().mesh().comm());
// --- Test for singularity
if (solverPerf.checkSingularity(mag(wApT)/normFactor))
......@@ -182,7 +186,7 @@ Foam::solverPerformance Foam::PBiCG::solve
// --- Update solution and residual:
scalar alpha = wArT/wApT;
const scalar alpha = wArT/wApT;
for (label cell=0; cell<nCells; cell++)
{
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 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 "PBiCGStab.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(PBiCGStab, 0);
lduMatrix::solver::addasymMatrixConstructorToTable<PBiCGStab>
addPBiCGStabAsymMatrixConstructorToTable_;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::PBiCGStab::PBiCGStab
(
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::PBiCGStab::solve
(
scalarField& psi,
const scalarField& source,
const direction cmpt
) const
{
// --- Setup class containing solver performance data
solverPerformance solverPerf
(
lduMatrix::preconditioner::getName(controlDict_) + typeName,
fieldName_
);
const label nCells = psi.size();
scalar* __restrict__ psiPtr = psi.begin();
scalarField pA(nCells);
scalar* __restrict__ pAPtr = pA.begin();
scalarField yA(nCells);
scalar* __restrict__ yAPtr = yA.begin();
// --- Calculate A.psi
matrix_.Amul(yA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
// --- Calculate initial residual field
scalarField rA(source - yA);
scalar* __restrict__ rAPtr = rA.begin();
// --- Calculate normalisation factor
const scalar normFactor = this->normFactor(psi, source, yA, pA);
if (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_)
)
{
scalarField AyA(nCells);
scalar* __restrict__ AyAPtr = AyA.begin();
scalarField sA(nCells);
scalar* __restrict__ sAPtr = sA.begin();
scalarField zA(nCells);
scalar* __restrict__ zAPtr = zA.begin();
scalarField tA(nCells);
scalar* __restrict__ tAPtr = tA.begin();
// --- Store initial residual
const scalarField rA0(rA);
// --- Initial values not used
scalar rA0rA = 0;
scalar alpha = 0;
scalar omega = 0;
// --- Select and construct the preconditioner
autoPtr<lduMatrix::preconditioner> preconPtr =
lduMatrix::preconditioner::New
(
*this,
controlDict_
);
// --- Solver iteration
do
{
// --- Store previous rA0rA
const scalar rA0rAold = rA0rA;
rA0rA = gSumProd(rA0, rA, matrix().mesh().comm());
// --- Test for singularity
if (solverPerf.checkSingularity(mag(rA0rA)))
{
break;
}
// --- Update pA
if (solverPerf.nIterations() == 0)
{
for (label cell=0; cell<nCells; cell++)
{
pAPtr[cell] = rAPtr[cell];
}
}
else
{
// --- Test for singularity
if (solverPerf.checkSingularity(mag(omega)))
{
break;
}
const scalar beta = (rA0rA/rA0rAold)*(alpha/omega);
for (label cell=0; cell<nCells; cell++)
{
pAPtr[cell] =
rAPtr[cell] + beta*(pAPtr[cell] - omega*AyAPtr[cell]);
}
}
// --- Precondition pA
preconPtr->precondition(yA, pA, cmpt);
// --- Calculate AyA
matrix_.Amul(AyA, yA, interfaceBouCoeffs_, interfaces_, cmpt);
const scalar rA0AyA = gSumProd(rA0, AyA, matrix().mesh().comm());
alpha = rA0rA/rA0AyA;
// --- Calculate sA
for (label cell=0; cell<nCells; cell++)
{
sAPtr[cell] = rAPtr[cell] - alpha*AyAPtr[cell];
}
// --- Test sA for convergence
solverPerf.finalResidual() =
gSumMag(sA, matrix().mesh().comm())/normFactor;
if (solverPerf.checkConvergence(tolerance_, relTol_))
{
for (label cell=0; cell<nCells; cell++)
{
psiPtr[cell] += alpha*yAPtr[cell];
}
solverPerf.nIterations()++;
return solverPerf;
}
// --- Precondition sA
preconPtr->precondition(zA, sA, cmpt);
// --- Calculate tA
matrix_.Amul(tA, zA, interfaceBouCoeffs_, interfaces_, cmpt);
const scalar tAtA = gSumSqr(tA, matrix().mesh().comm());
// --- Calculate omega from tA and sA
// (cheaper than using zA with preconditioned tA)
omega = gSumProd(tA, sA)/tAtA;
// --- Update solution and residual
for (label cell=0; cell<nCells; cell++)
{
psiPtr[cell] += alpha*yAPtr[cell] + omega*zAPtr[cell];
rAPtr[cell] = sAPtr[cell] - omega*tAPtr[cell];
}
solverPerf.finalResidual() =
gSumMag(rA, matrix().mesh().comm())
/normFactor;
} while
(
(
solverPerf.nIterations()++ < maxIter_
&& !solverPerf.checkConvergence(tolerance_, relTol_)
)
|| solverPerf.nIterations() < minIter_
);
}
return solverPerf;
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 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/>.
Class
Foam::PBiCGStab
Description
Preconditioned bi-conjugate gradient stabilized solver for asymmetric
lduMatrices using a run-time selectable preconditiioner.
References:
\verbatim
Van der Vorst, H. A. (1992).
Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG
for the solution of nonsymmetric linear systems.
SIAM Journal on scientific and Statistical Computing, 13(2), 631-644.
Barrett, R., Berry, M. W., Chan, T. F., Demmel, J., Donato, J.,
Dongarra, J., Eijkhout, V., Pozo, R., Romine, C. & Van der Vorst, H.
(1994).
Templates for the solution of linear systems:
building blocks for iterative methods
(Vol. 43). Siam.
\endverbatim
SourceFiles
PBiCGStab.C
\*---------------------------------------------------------------------------*/
#ifndef PBiCGStab_H
#define PBiCGStab_H
#include "lduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class PBiCGStab Declaration
\*---------------------------------------------------------------------------*/
class PBiCGStab
:
public lduMatrix::solver
{
// Private Member Functions
//- Disallow default bitwise copy construct
PBiCGStab(const PBiCGStab&);
//- Disallow default bitwise assignment
void operator=(const PBiCGStab&);
public:
//- Runtime type information
TypeName("PBiCGStab");
// Constructors
//- Construct from matrix components and solver data stream
PBiCGStab
(
const word& fieldName,
const lduMatrix& matrix,
const FieldField<Field, scalar>& interfaceBouCoeffs,
const FieldField<Field, scalar>& interfaceIntCoeffs,
const lduInterfaceFieldPtrsList& interfaces,
const dictionary& solverControls
);
//- Destructor
virtual ~PBiCGStab()
{}
// Member Functions
//- Solve the matrix with this solver
virtual solverPerformance solve
(
scalarField& psi,
const scalarField& source,
const direction cmpt=0
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
......@@ -19,7 +19,7 @@ solvers
{
T
{
solver PBiCG;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-06;
relTol 0;
......
......@@ -49,7 +49,7 @@ solvers
"(b|Xi|ft|ha|hau|k|epsilon)"
{
solver PBiCG;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-06;
relTol 0.1;
......@@ -57,7 +57,7 @@ solvers
"(b|Xi|ft|ha|hau|k|epsilon)Final"
{
solver PBiCG;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-06;
relTol 0;
......
......@@ -34,7 +34,7 @@ solvers
"(U|b|Su|Xi|ha|hau|k|epsilon)"
{
solver PBiCG;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
......@@ -42,7 +42,7 @@ solvers
"(U|b|Su|Xi|ha|hau|k|epsilon)Final"
{
solver PBiCG;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-05;
relTol 0;
......
......@@ -19,7 +19,7 @@ solvers
{
Yi
{
solver PBiCG;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-12;
relTol 0;
......
......@@ -19,7 +19,7 @@ solvers
{
Yi
{
solver PBiCG;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-12;
relTol 0;
......
......@@ -19,7 +19,7 @@ solvers
{
Yi
{
solver PBiCG;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-12;
relTol 0;
......
......@@ -19,7 +19,7 @@ solvers
{
Yi
{
solver PBiCG;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-12;
relTol 0;
......
......@@ -19,7 +19,7 @@ solvers
{
Yi
{
solver PBiCG;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-12;
relTol 0;
......
......@@ -49,7 +49,7 @@ solvers
"(U|Xi|hau|eau|ft|b|ha|ea|k|epsilon)"
{
solver PBiCG;
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
......
......@@ -19,21 +19,21 @@ solvers
{
hf
{
solver PBiCG;
solver PBiCGStab;
preconditioner DILU;
tolerance 0;
relTol 1e-3;
}
"(Uf|deltaf\*rhof)"
{
solver PBiCG;
solver PBiCGStab;
preconditioner DILU;