Merge branch 'development' of ssh://git.code.sf.net/u/hjasak/foam-extend-4.0 into development

src/lduSolvers/lduPrecon/CholeskyPrecon/CholeskyPrecon.C

@@ -65,13 +65,12 @@ void Foam::CholeskyPrecon::calcPreconDiag()
 
                 // Get interface coefficiens
                 const scalarField& bouCoeffs = coupleBouCoeffs_[patchI];
-                const scalarField& intCoeffs = coupleIntCoeffs_[patchI];
 
                 forAll (fc, coeffI)
                 {
+                    // Note change of sign for boundary coeffs
                     preconDiag_[fc[coeffI]] +=
-                        bouCoeffs[coeffI]*intCoeffs[coeffI]/
+                        sqr(bouCoeffs[coeffI])/preconDiag_[fc[coeffI]];
-                        preconDiag_[fc[coeffI]];
                 }
             }
         }
@@ -174,39 +173,7 @@ void Foam::CholeskyPrecon::precondition
             << abort(FatalError);
     }
 
-    // In order to properly execute parallel preconditioning, re-use
+    // Diagonal block
-    // x to zero and execute coupled boundary update
-    // HJ, 19/Jun/2017
-
-    x = 0;
-
-    // Coupled boundary update
-    {
-        matrix_.initMatrixInterfaces
-        (
-            coupleBouCoeffs_,
-            interfaces_,
-            x,
-            x,               // put result into x
-            cmpt,
-            false
-        );
-
-        matrix_.updateMatrixInterfaces
-        (
-            coupleBouCoeffs_,
-            interfaces_,
-            x,
-            x,               // put result into x
-            cmpt,
-            false
-        );
-    }
-
-    // Multiply with inverse diag to precondition
-    x *= preconDiag_;
-
-    // Diagonal block: note +=
     {
         scalar* __restrict__ xPtr = x.begin();
 
@@ -218,7 +185,7 @@ void Foam::CholeskyPrecon::precondition
 
         for (register label rowI = 0; rowI < nRows; rowI++)
         {
-            xPtr[rowI] += bPtr[rowI]*preconDiagPtr[rowI];
+            xPtr[rowI] = bPtr[rowI]*preconDiagPtr[rowI];
         }
     }
 

src/lduSolvers/lduPrecon/ILU0/ILU0.C

@@ -65,12 +65,18 @@ void Foam::ILU0::calcPreconDiag()
 
                 // Get interface coefficiens
                 const scalarField& bouCoeffs = coupleBouCoeffs_[patchI];
-                const scalarField& intCoeffs = coupleIntCoeffs_[patchI];
 
+                // Note:
+                // In order to do the preconditioning correctly, the lower
+                // triangular coefficient on the coupled interface should
+                // also be available.  Since it is not, we will re-create it
+                // assuming the negSumDiag rule
+                
                 forAll (fc, coeffI)
                 {
-                    preconDiag_[fc[coeffI]] +=
+                    // Note consistent sign for boundary coeffs
-                        bouCoeffs[coeffI]*intCoeffs[coeffI]/
+                    preconDiag_[fc[coeffI]] -=
+                        bouCoeffs[coeffI]*(1 - bouCoeffs[coeffI])/
                         preconDiag_[fc[coeffI]];
                 }
             }
@@ -179,33 +185,31 @@ void Foam::ILU0::precondition
     // x to zero and execute coupled boundary update first
     // HJ, 19/Jun/2017
 
-    x = 0;
+    // // Coupled boundary update
+    // {
+    //     matrix_.initMatrixInterfaces
+    //     (
+    //         coupleBouCoeffs_,
+    //         interfaces_,
+    //         x,
+    //         x,               // put result into x
+    //         cmpt,
+    //         false
+    //     );
 
-    // Coupled boundary update
+    //     matrix_.updateMatrixInterfaces
-    {
+    //     (
-        matrix_.initMatrixInterfaces
+    //         coupleBouCoeffs_,
-        (
+    //         interfaces_,
-            coupleBouCoeffs_,
+    //         x,
-            interfaces_,
+    //         x,               // put result into x
-            x,
+    //         cmpt,
-            x,               // put result into x
+    //         false
-            cmpt,
+    //     );
-            false
+    // }
-        );
 
-        matrix_.updateMatrixInterfaces
+    // // Multiply with inverse diag to precondition
-        (
+    // x *= preconDiag_;
-            coupleBouCoeffs_,
-            interfaces_,
-            x,
-            x,               // put result into x
-            cmpt,
-            false
-        );
-    }
-
-    // Multiply with inverse diag to precondition
-    x *= preconDiag_;
 
     // Diagonal block: note +=
     forAll (x, i)
@@ -270,38 +274,10 @@ void Foam::ILU0::preconditionT
     // x to zero and execute coupled boundary update first
     // HJ, 19/Jun/2017
 
-    x = 0;
+    // Diagonal block
-
+    forAll (x, i)
-    // Coupled boundary update
     {
-        matrix_.initMatrixInterfaces
+        x[i] = b[i]*preconDiag_[i];
-        (
-            coupleIntCoeffs_, // Note: transpose coupled coeffs
-            interfaces_,
-            x,
-            x,                // put result into x
-            cmpt,
-            false
-        );
-
-        matrix_.updateMatrixInterfaces
-        (
-            coupleIntCoeffs_, // Note: transpose coupled coeffs
-            interfaces_,
-            x,
-            x,                // put result into x
-            cmpt,
-            false
-        );
-    }
-
-    // Multiply with inverse diag to precondition
-    x *= preconDiag_;
-
-    // Diagonal block: note +=
-    forAll(x, i)
-    {
-        x[i] += b[i]*preconDiag_[i];
     }
 
     if (matrix_.asymmetric())

src/lduSolvers/lduSmoother/iluSmoother/iluSmoother.C

@@ -25,7 +25,7 @@ Class
     iluSmoother
 
 Description
-    Symmetric Gauss-Seidel smoother
+    ILU smoother
 
 Author
     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
@@ -33,6 +33,8 @@ Author
 \*---------------------------------------------------------------------------*/
 
 #include "iluSmoother.H"
+#include "CholeskyPrecon.H"
+#include "ILU0.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -66,16 +68,37 @@ Foam::iluSmoother::iluSmoother
         coupleIntCoeffs,
         interfaces
     ),
-    precon_
+    preconPtr_(),
-    (
-        matrix,
-        coupleBouCoeffs,
-        coupleIntCoeffs,
-        interfaces
-    ),
     xCorr_(matrix.lduAddr().size()),
     residual_(matrix.lduAddr().size())
-{}
+{
+    if (matrix.symmetric())
+    {
+        preconPtr_.set
+        (
+            new CholeskyPrecon
+            (
+                matrix,
+                coupleBouCoeffs,
+                coupleIntCoeffs,
+                interfaces
+            )
+        );
+    }
+    else if (matrix.asymmetric())
+    {
+        preconPtr_.set
+        (
+            new ILU0
+            (
+                matrix,
+                coupleBouCoeffs,
+                coupleIntCoeffs,
+                interfaces
+            )
+        );
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -106,7 +129,10 @@ void Foam::iluSmoother::smooth
             residual_[i] = b[i] - residual_[i];
         }
 
-        precon_.precondition(xCorr_, residual_, cmpt);
+        if (!matrix_.diagonal())
+        {
+            preconPtr_->precondition(xCorr_, residual_, cmpt);
+        }
 
         // Add correction to x
         x += xCorr_;

src/lduSolvers/lduSmoother/iluSmoother/iluSmoother.H

@@ -38,7 +38,7 @@ SourceFiles
 #ifndef iluSmoother_H
 #define iluSmoother_H
 
-#include "CholeskyPrecon.H"
+#include "lduMatrix.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,7 +56,7 @@ class iluSmoother
     // Private Data
 
         //- Cholesky preconditioner
-        CholeskyPrecon precon_;
+        autoPtr<lduPreconditioner> preconPtr_;
 
         //- Correction array
         mutable scalarField xCorr_;
@@ -92,10 +92,9 @@ public:
         );
 
 
-    // Destructor
+    //- Destructor
-
+    virtual ~iluSmoother()
-        virtual ~iluSmoother()
+    {}
-        {}
 
 
     // Member Functions