From e435bf40f6cf0298eba86ef368b24121aa963099 Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Sun, 27 May 2018 16:25:14 +0200
Subject: [PATCH] Change top-level AMG solver settings: minIter = 0 and reset
 top x in each call

---
 .../BlockAMG/coarseBlockAMGLevel.C            | 25 ++++++++++---------
 src/lduSolvers/amg/coarseAmgLevel.C           | 13 +++++-----
 2 files changed, 19 insertions(+), 19 deletions(-)

diff --git a/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.C b/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.C
index 790cfd0dd..6016125d9 100644
--- a/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.C
+++ b/src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.C
@@ -39,6 +39,7 @@ Author
 #include "BlockSolverPerformance.H"
 #include "BlockCGSolver.H"
 #include "BlockBiCGStabSolver.H"
+#include "BlockGMRESSolver.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -56,8 +57,8 @@ Foam::coarseBlockAMGLevel<Type>::coarseBlockAMGLevel
 :
     addrPtr_(addrPtr),
     matrixPtr_(matrixPtr),
-    x_(matrixPtr_->diag().size(),pTraits<Type>::zero),
-    b_(matrixPtr_->diag().size(),pTraits<Type>::zero),
+    x_(matrixPtr_->diag().size(), pTraits<Type>::zero),
+    b_(matrixPtr_->diag().size(), pTraits<Type>::zero),
     dict_(dict),
     coarseningPtr_
     (
@@ -215,15 +216,15 @@ void Foam::coarseBlockAMGLevel<Type>::solve
     // Create artificial dictionary for top-level solution
     dictionary topLevelDict;
     topLevelDict.add("nDirections", "5");
-    topLevelDict.add("minIter", 1);
+    topLevelDict.add("preconditioner", "ILUC0");
+    topLevelDict.add("minIter", 0);
     topLevelDict.add("maxIter", 500);
     topLevelDict.add("tolerance", tolerance);
     topLevelDict.add("relTol", relTol);
 
     // Avoid issues with round-off on strict tolerance setup
     // HJ, 27/Jun/2013
-    // Create multiplication function object
-    typename BlockCoeff<Type>::multiply mult;
+    x = pTraits<Type>::zero;
 
     // Switch of debug in top-level direct solve
     label oldDebug = blockLduMatrix::debug();
@@ -243,9 +244,6 @@ void Foam::coarseBlockAMGLevel<Type>::solve
 
     if (matrixPtr_->symmetric())
     {
-        // Note: must change preconditioner to C0.  HJ. 10/Oct/2017
-        topLevelDict.add("preconditioner", "ILUC0");
-
         coarseSolverPerf = BlockCGSolver<Type>
         (
             "topLevelCorr",
@@ -255,8 +253,6 @@ void Foam::coarseBlockAMGLevel<Type>::solve
     }
     else
     {
-        topLevelDict.add("preconditioner", "ILUC0");
-
         coarseSolverPerf =
         BlockBiCGStabSolver<Type>
         (
@@ -267,10 +263,15 @@ void Foam::coarseBlockAMGLevel<Type>::solve
     }
 
     // Check for convergence
+
+//    blockLduMatrix::debug = 1;
+//    coarseSolverPerf.print();
+//    Info<< "b: " << gSum(b) << " corr (min, max): (" << gMin(x) << ", " << gMax(x) << endl;
+
     const scalar magInitialRes = mag(coarseSolverPerf.initialResidual());
     const scalar magFinalRes = mag(coarseSolverPerf.finalResidual());
 
-    if (magFinalRes > magInitialRes && magInitialRes > SMALL)
+    if (magFinalRes > magInitialRes && magInitialRes > 1e-12)
     {
         if (blockLduMatrix::debug)
         {
@@ -284,7 +285,7 @@ void Foam::coarseBlockAMGLevel<Type>::solve
     // Restore debug
     blockLduMatrix::debug = oldDebug;
 
-    if (blockLduMatrix::debug >= 2)
+    if (blockLduMatrix::debug >= 3)
     {
         coarseSolverPerf.print();
     }
diff --git a/src/lduSolvers/amg/coarseAmgLevel.C b/src/lduSolvers/amg/coarseAmgLevel.C
index 5e25384dc..8a79047fb 100644
--- a/src/lduSolvers/amg/coarseAmgLevel.C
+++ b/src/lduSolvers/amg/coarseAmgLevel.C
@@ -194,11 +194,15 @@ void Foam::coarseAmgLevel::solve
     lduSolverPerformance coarseSolverPerf;
 
     dictionary topLevelDict;
-    topLevelDict.add("minIter", 1);
+    topLevelDict.add("preconditioner", "ILUC0");
+    topLevelDict.add("minIter", 0);
     topLevelDict.add("maxIter", 500);
     topLevelDict.add("tolerance", tolerance);
     topLevelDict.add("relTol", relTol);
 
+    // Top-level round-off error control.  HJ, 28/May/2018
+    x = 0;
+
     // Switch off debug in top-level direct solve
     label oldDebug = blockLduMatrix::debug();
 
@@ -217,9 +221,6 @@ void Foam::coarseAmgLevel::solve
 
     if (matrixPtr_->matrix().symmetric())
     {
-        // Note: must change preconditioner to C0.  HJ. 10/Oct/2017
-        topLevelDict.add("preconditioner", "ILUC0");
-
         coarseSolverPerf = cgSolver
         (
             "topLevelCorr",
@@ -232,8 +233,6 @@ void Foam::coarseAmgLevel::solve
     }
     else
     {
-        topLevelDict.add("preconditioner", "ILUC0");
-
         coarseSolverPerf = bicgStabSolver
         (
             "topLevelCorr",
@@ -249,7 +248,7 @@ void Foam::coarseAmgLevel::solve
     const scalar magInitialRes = mag(coarseSolverPerf.initialResidual());
     const scalar magFinalRes = mag(coarseSolverPerf.finalResidual());
 
-    if (magFinalRes > magInitialRes && magInitialRes > SMALL)
+    if (magFinalRes > magInitialRes && magInitialRes > 1e-12)
     {
         if (blockLduMatrix::debug)
         {