From ccefe8bb902133160f28c53190e58da90ad2b563 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 9 Nov 2011 12:31:57 +0000
Subject: [PATCH] ENH: fvMatrix: debug printing of diagonal dominance

---
 .../fvMatrices/fvMatrix/fvMatrix.C            | 34 +++++++++++++++++++
 .../fvMatrices/fvMatrix/fvMatrix.H            |  1 +
 2 files changed, 35 insertions(+)

diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 33922b5019..fd78acb546 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -587,7 +587,41 @@ void Foam::fvMatrix<Type>::relax(const scalar alpha)
         }
     }
 
+
+    if (debug)
+    {
+        // Calculate amount of non-dominance.
+        label nNon = 0;
+        scalar maxNon = 0.0;
+        scalar sumNon = 0.0;
+        forAll(D, celli)
+        {
+            scalar d = (sumOff[celli] - D[celli])/D[celli];
+
+            if (d > 0)
+            {
+                nNon++;
+                maxNon = max(maxNon, d);
+                sumNon += d;
+            }
+        }
+
+        reduce(nNon, sumOp<label>());
+        reduce(maxNon, maxOp<scalar>());
+        reduce(sumNon, sumOp<scalar>());
+        sumNon /= returnReduce(D.size(), sumOp<label>());
+
+        InfoIn("fvMatrix<Type>::relax(const scalar alpha)")
+            << "Matrix dominance test for " << psi_.name() << nl
+            << "    number of non-dominant cells   : " << nNon << nl
+            << "    maximum relative non-dominance : " << maxNon << nl
+            << "    average relative non-dominance : " << sumNon << nl
+            << endl;
+    }
+
+
     // Ensure the matrix is diagonally dominant...
+    // (assumes that the central coefficient is positive)
     max(D, D, sumOff);
 
     // ... then relax
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
index 196a48285e..25021653f3 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
@@ -364,6 +364,7 @@ public:
             //  alpha = 1 : diagonally equal
             //  alpha < 1 : diagonally dominant
             //  alpha = 0 : do nothing
+            //  Note: Requires positive diagonal.
             void relax(const scalar alpha);
 
             //- Relax matrix (for steady-state solution).