From 7c8d3866cb5fb2d4fe399d674c33d74a516dc1a0 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Fri, 3 Feb 2017 21:53:53 +0000
Subject: [PATCH] PBiCG: Suggest changing to the more robust PBiCGStab solver

if convergence is not achieved within the maximum number of iterations.

Sometimes, particularly running in parallel, PBiCG fails to converge or diverges
without warning or obvious cause leaving a solution field containing significant
errors which can cause divergence of the application.  PBiCGStab is more robust
and does not suffer from the problems encountered with PBiCG.
---
 .../matrices/lduMatrix/lduMatrix/lduMatrix.C         |  5 ++++-
 .../matrices/lduMatrix/lduMatrix/lduMatrix.H         |  5 ++++-
 .../matrices/lduMatrix/lduMatrix/lduMatrixSolver.C   |  8 ++++----
 .../matrices/lduMatrix/solvers/PBiCG/PBiCG.C         | 12 +++++++++++-
 4 files changed, 23 insertions(+), 7 deletions(-)

diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C
index 8de739562..373a1a569 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,9 @@ namespace Foam
 }
 
 
+const Foam::label Foam::lduMatrix::solver::defaultMaxIter_ = 1000;
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 Foam::lduMatrix::lduMatrix(const lduMesh& mesh)
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
index 575916f63..30ba325a5 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -106,6 +106,9 @@ public:
             //- Dictionary of controls
             dictionary controlDict_;
 
+            //- Default maximum number of iterations in the solver
+            static const label defaultMaxIter_;
+
             //- Maximum number of iterations in the solver
             label maxIter_;
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C
index 1f6bd82d1..1466fa816 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -157,10 +157,10 @@ Foam::lduMatrix::solver::solver
 
 void Foam::lduMatrix::solver::readControls()
 {
-    maxIter_   = controlDict_.lookupOrDefault<label>("maxIter", 1000);
-    minIter_   = controlDict_.lookupOrDefault<label>("minIter", 0);
+    maxIter_ = controlDict_.lookupOrDefault<label>("maxIter", defaultMaxIter_);
+    minIter_ = controlDict_.lookupOrDefault<label>("minIter", 0);
     tolerance_ = controlDict_.lookupOrDefault<scalar>("tolerance", 1e-6);
-    relTol_    = controlDict_.lookupOrDefault<scalar>("relTol", 0);
+    relTol_ = controlDict_.lookupOrDefault<scalar>("relTol", 0);
 }
 
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C
index 3d43c64a0..61f117c12 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -204,6 +204,16 @@ Foam::solverPerformance Foam::PBiCG::solve
         );
     }
 
+    // Recommend PBiCGStab if PBiCG fails to converge
+    if (solverPerf.nIterations() > max(defaultMaxIter_, maxIter_))
+    {
+        FatalErrorInFunction
+            << "PBiCG has failed to converge within the maximum number"
+               " of iterations " << max(defaultMaxIter_, maxIter_) << nl
+            << "    Please try the more robust PBiCGStab solver."
+            << exit(FatalError);
+    }
+
     return solverPerf;
 }