From 2bc20518c4139544cd055aabe177ad78b4f72342 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Fri, 3 Feb 2017 18:36:40 +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.H         |  5 ++++-
 .../matrices/LduMatrix/LduMatrix/LduMatrixSolver.C   |  4 ++--
 .../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 +++++++++++-
 6 files changed, 29 insertions(+), 10 deletions(-)

diff --git a/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H b/src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H
index 1da27f106b..f96d583b1f 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
@@ -122,6 +122,9 @@ public:
             //- Dictionary of controls
             dictionary controlDict_;
 
+            //- Default maximum number of iterations in the solver
+            static const label defaultMaxIter_ = 1000;
+
             //- 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 651be6cee2..dd908b91f0 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-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -131,7 +131,7 @@ Foam::LduMatrix<Type, DType, LUType>::solver::solver
 
     controlDict_(solverDict),
 
-    maxIter_(1000),
+    maxIter_(defaultMaxIter_),
     minIter_(0),
     tolerance_(1e-6*pTraits<Type>::one),
     relTol_(Zero)
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C
index 11cdf31485..e9eb25143a 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 575916f633..30ba325a5a 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 d554375cae..bc1e134741 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-2016 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 54088f550f..b4e630571f 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-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -208,6 +208,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;
 }