diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index c975296d1..3bd8bcaa0 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -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
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C
new file mode 100644
index 000000000..ff7c0cddd
--- /dev/null
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.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 .
+
+\*---------------------------------------------------------------------------*/
+
+#include "PBiCGStab.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(PBiCGStab, 0);
+
+ lduMatrix::solver::addasymMatrixConstructorToTable
+ addPBiCGStabAsymMatrixConstructorToTable_;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::PBiCGStab::PBiCGStab
+(
+ 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
+ )
+{}
+
+
+// * * * * * * * * * * * * * * * 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 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; cellprecondition(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; cellprecondition(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.
+
+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& interfaceBouCoeffs,
+ const FieldField& 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
+
+// ************************************************************************* //