Skip to content
Permalink
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.
  • Loading branch information...
Henry Weller
Henry Weller committed Sep 12, 2016
1 parent 168b29e commit 1ca4bbf1d288075a76940d58fdf61cbfc49298ee
@@ -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
@@ -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;
}


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

// ************************************************************************* //

0 comments on commit 1ca4bbf

Please sign in to comment.
You can’t perform that action at this time.