research/OpenFOAM-4.x
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

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

Browse source 
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 commit is contained in:
Henry Weller 2016-09-12 15:19:24 +01:00
parent 168b29e2cf
commit 1ca4bbf1d2
3 changed files with 376 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
src/OpenFOAM/Make/files
View file

@ -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

252
src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C Normal file
View file

@ -0,0 +1,252 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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;
}
// ************************************************************************* //

123
src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.H Normal file
View file

@ -0,0 +1,123 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 preconditioner.
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
// ************************************************************************* //