Improved normalisation of reported residuals for block-coupled systems

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAMGSolver/BlockAMGSolver.C

@@ -72,10 +72,10 @@ Foam::BlockAMGSolver<Type>::solve
         this->fieldName()
     );
 
-    scalar norm = this->normFactor(x, b);
+    Type norm = this->normFactor(x, b);
 
     // Calculate initial residual
-    solverPerf.initialResidual() = gSum(cmptMag(amg_.residual(x, b)))/norm;
+    solverPerf.initialResidual() = cmptDivide(gSum(cmptMag(amg_.residual(x, b))),norm);
     solverPerf.finalResidual() = solverPerf.initialResidual();
 
     // Stop solver on divergence
@@ -88,8 +88,8 @@ Foam::BlockAMGSolver<Type>::solve
         {
             amg_.cycle(x, b);
 
-            solverPerf.finalResidual() =
+              solverPerf.finalResidual() =
-                gSum(cmptMag(amg_.residual(x, b)))/norm;
+                  cmptDivide(gSum(cmptMag(amg_.residual(x, b))), norm);
 
             solverPerf.nIterations()++;
 

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C

@@ -76,7 +76,7 @@ Foam::BlockBiCGStabSolver<Type>::solve
         this->fieldName()
     );
 
-    scalar norm = this->normFactor(x, b);
+    Type norm = this->normFactor(x, b);
 
     // Multiplication helper
     typename BlockCoeff<Type>::multiply mult;
@@ -87,7 +87,8 @@ Foam::BlockBiCGStabSolver<Type>::solve
     matrix.Amul(p, x);
     Field<Type> r(b - p);
 
-    solverPerf.initialResidual() = gSum(cmptMag(r))/norm;
+    // NOTE: Normalisation of residual per component! TU, Feb 2019
+    solverPerf.initialResidual() = cmptDivide(gSum(cmptMag(r)),norm);
     solverPerf.finalResidual() = solverPerf.initialResidual();
 
     // Check convergence, solve if not converged
@@ -160,7 +161,7 @@ Foam::BlockBiCGStabSolver<Type>::solve
                 r[i] = s[i] - omega*t[i];
             }
 
-            solverPerf.finalResidual() = gSum(cmptMag(r))/norm;
+            solverPerf.finalResidual() = cmptDivide(gSum(cmptMag(r)), norm);
             solverPerf.nIterations()++;
         } while (!this->stop(solverPerf));
     }

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C

@@ -75,7 +75,7 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
         this->fieldName()
     );
 
-    scalar norm = this->normFactor(x, b);
+    Type norm = this->normFactor(x, b);
 
     Field<Type> wA(x.size());
 
@@ -83,7 +83,8 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
     matrix.Amul(wA, x);
     Field<Type> rA(b - wA);
 
-    solverPerf.initialResidual() = gSum(cmptMag(rA))/norm;
+    // NOTE: Normalisation of residual per component! TU, Feb 2019
+    solverPerf.initialResidual() = cmptDivide(gSum(cmptMag(rA)),norm);
     solverPerf.finalResidual() = solverPerf.initialResidual();
 
     // Check convergence, solve if not converged
@@ -120,7 +121,7 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
             wApA = gSumProd(wA, pA);
 
             // Check for singularity
-            if (solverPerf.checkSingularity(mag(wApA)/norm))
+            if (solverPerf.checkSingularity(mag(wApA)/mag(norm)))
             {
                 break;
             }
@@ -138,7 +139,7 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
                 rA[i] -= alpha*wA[i];
             }
 
-            solverPerf.finalResidual() = gSum(cmptMag(rA))/norm;
+            solverPerf.finalResidual() = cmptDivide(gSum(cmptMag(rA)), norm);
             solverPerf.nIterations()++;
         } while (!this->stop(solverPerf));
     }

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C

@@ -111,7 +111,7 @@ Foam::BlockGMRESSolver<Type>::solve
         this->fieldName()
     );
 
-    scalar norm = this->normFactor(x, b);
+    Type norm = this->normFactor(x, b);
 
     Field<Type> wA(x.size());
 
@@ -119,7 +119,7 @@ Foam::BlockGMRESSolver<Type>::solve
     matrix.Amul(wA, x);
     Field<Type> rA(b - wA);
 
-    solverPerf.initialResidual() = gSum(cmptMag(rA))/norm;
+    solverPerf.initialResidual() = cmptDivide(gSum(cmptMag(rA)),norm);
     solverPerf.finalResidual() = solverPerf.initialResidual();
 
     // Check convergence, solve if not converged
@@ -234,7 +234,7 @@ Foam::BlockGMRESSolver<Type>::solve
                 rA[raI] = b[raI] - wA[raI];
             }
 
-            solverPerf.finalResidual() = gSum(cmptMag(rA))/norm;
+            solverPerf.finalResidual() = cmptDivide(gSum(cmptMag(rA)), norm);
             solverPerf.nIterations()++;
         } while (!this->stop(solverPerf));
     }

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C

@@ -72,7 +72,7 @@ Foam::BlockGaussSeidelSolver<Type>::solve
         this->fieldName()
     );
 
-    scalar norm = this->normFactor(x, b);
+    Type norm = this->normFactor(x, b);
 
     Field<Type> wA(x.size());
 
@@ -80,7 +80,7 @@ Foam::BlockGaussSeidelSolver<Type>::solve
     matrix.Amul(wA, x);
     wA -= b;
 
-    solverPerf.initialResidual() = gSum(cmptMag(wA))/norm;
+    solverPerf.initialResidual() = cmptDivide(gSum(cmptMag(wA)),norm);
     solverPerf.finalResidual() = solverPerf.initialResidual();
 
     // Check convergence, solve if not converged
@@ -103,7 +103,7 @@ Foam::BlockGaussSeidelSolver<Type>::solve
             matrix.Amul(wA, x);
             wA -= b;
 
-            solverPerf.finalResidual() = gSum(cmptMag(wA))/norm;
+            solverPerf.finalResidual() = cmptDivide(gSum(cmptMag(wA)), norm);
         } while (!this->stop(solverPerf));
     }
 

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockILU/BlockILUSolver.C

@@ -72,7 +72,7 @@ Foam::BlockILUSolver<Type>::solve
         this->fieldName()
     );
 
-    scalar norm = this->normFactor(x, b);
+    Type norm = this->normFactor(x, b);
 
     Field<Type> residual(x.size());
 
@@ -84,7 +84,7 @@ Foam::BlockILUSolver<Type>::solve
         residual[i] = b[i] - residual[i];
     }
 
-    solverPerf.initialResidual() = gSum(cmptMag(residual))/norm;
+    solverPerf.initialResidual() = cmptDivide(gSum(cmptMag(residual)),norm);
     solverPerf.finalResidual() = solverPerf.initialResidual();
 
     // Check convergence, solve if not converged
@@ -114,7 +114,7 @@ Foam::BlockILUSolver<Type>::solve
                 residual[i] = b[i] - residual[i];
             }
 
-            solverPerf.finalResidual() = gSum(cmptMag(residual))/norm;
+            solverPerf.finalResidual() = cmptDivide(gSum(cmptMag(residual)), norm);
         } while (!this->stop(solverPerf));
     }
 

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.C

@@ -50,7 +50,7 @@ Foam::BlockIterativeSolver<Type>::BlockIterativeSolver
 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 
 template<class Type>
-Foam::scalar Foam::BlockIterativeSolver<Type>::normFactor
+Type Foam::BlockIterativeSolver<Type>::normFactor
 (
     Field<Type>& x,
     const Field<Type>& b
@@ -77,7 +77,13 @@ Foam::scalar Foam::BlockIterativeSolver<Type>::normFactor
         Field<Type>(nRows, xRef)
     );
 
-    scalar normFactor = gSum(mag(wA - pA) + mag(b - pA)) + this->small_;
+    // The components of the normalisation factor cannot be equal to 0 since it
+    // appears in the denominator: replace zero with small_
+    Type small = pTraits<Type>::one*this->small_;
+
+    // Norm factor for the block matrix is not a scalar since each component
+    // should be normalised with a value corresponding to a particular variable
+    Type normFactor = gSum(cmptMag(wA - pA) + cmptMag(b - pA)) + small;
 
     if (blockLduMatrix::debug >= 2)
     {

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.H

@@ -84,7 +84,7 @@ protected:
     // Protected Member Functions
 
         //- Return normalisation factor
-        scalar normFactor
+        Type normFactor
         (
             Field<Type>& x,
             const Field<Type>& b
@@ -107,10 +107,9 @@ public:
         );
 
 
-    // Destructor
+    //- Destructor
-
+    virtual ~BlockIterativeSolver()
-        virtual ~BlockIterativeSolver()
+    {}
-        {}
 
 
     // Member Functions

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformance.C

@@ -39,14 +39,9 @@ bool Foam::BlockSolverPerformance<Type>::checkConvergence
     {
         Info<< solverName_
             << ":  Iteration " << nIterations_
-            << " residual = " << finalResidual_;
+            << " residual = " << finalResidual_
-
+            << " mag = " << mag(finalResidual_)
-        if (pTraits<Type>::rank > 0)
+            << " tol = "
-        {
-            Info<< " mag = " << mag(finalResidual_);
-        }
-
-        Info<< " tol = "
             << Foam::max(Tolerance, RelTolerance*mag(initialResidual_))
             << endl;
     }

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockLduSolver/BlockSolverPerformances.C

@@ -41,6 +41,7 @@ defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceVector, 1);
 defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceSphericalTensor, 1);
 defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceSymmTensor, 1);
 defineNamedTemplateTypeNameAndDebug(BlockSolverPerformanceTensor, 1);
+
 };