From 1354e6971c0277642744d0f5e72c20501d18cca7 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 19 Jun 2017 17:18:50 +0100 Subject: [PATCH 1/2] Update iluSmoother for symmetric and asymmetric matrices --- .../lduSmoother/iluSmoother/iluSmoother.C | 46 +++++++++++++++---- .../lduSmoother/iluSmoother/iluSmoother.H | 11 ++--- 2 files changed, 41 insertions(+), 16 deletions(-) diff --git a/src/lduSolvers/lduSmoother/iluSmoother/iluSmoother.C b/src/lduSolvers/lduSmoother/iluSmoother/iluSmoother.C index 570cf6000..bd48c4482 100644 --- a/src/lduSolvers/lduSmoother/iluSmoother/iluSmoother.C +++ b/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_ - ( - matrix, - coupleBouCoeffs, - coupleIntCoeffs, - interfaces - ), + preconPtr_(), 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_; diff --git a/src/lduSolvers/lduSmoother/iluSmoother/iluSmoother.H b/src/lduSolvers/lduSmoother/iluSmoother/iluSmoother.H index 63313358f..040945a3d 100644 --- a/src/lduSolvers/lduSmoother/iluSmoother/iluSmoother.H +++ b/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 preconPtr_; //- Correction array mutable scalarField xCorr_; @@ -92,10 +92,9 @@ public: ); - // Destructor - - virtual ~iluSmoother() - {} + //- Destructor + virtual ~iluSmoother() + {} // Member Functions From d17c293e6d6e092bd7066e1963d6cc77b38cbeb0 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Mon, 19 Jun 2017 17:19:17 +0100 Subject: [PATCH 2/2] Cholesky and ILU0 precon on coupled patches: bug fix --- .../lduPrecon/CholeskyPrecon/CholeskyPrecon.C | 41 +------- src/lduSolvers/lduPrecon/ILU0/ILU0.C | 94 +++++++------------ 2 files changed, 39 insertions(+), 96 deletions(-) diff --git a/src/lduSolvers/lduPrecon/CholeskyPrecon/CholeskyPrecon.C b/src/lduSolvers/lduPrecon/CholeskyPrecon/CholeskyPrecon.C index 3386f7c59..40c90d034 100644 --- a/src/lduSolvers/lduPrecon/CholeskyPrecon/CholeskyPrecon.C +++ b/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]/ - preconDiag_[fc[coeffI]]; + sqr(bouCoeffs[coeffI])/preconDiag_[fc[coeffI]]; } } } @@ -174,39 +173,7 @@ void Foam::CholeskyPrecon::precondition << abort(FatalError); } - // In order to properly execute parallel preconditioning, re-use - // 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 += + // Diagonal block { 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]; } } diff --git a/src/lduSolvers/lduPrecon/ILU0/ILU0.C b/src/lduSolvers/lduPrecon/ILU0/ILU0.C index cc0c87d9c..3bac83b89 100644 --- a/src/lduSolvers/lduPrecon/ILU0/ILU0.C +++ b/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]] += - bouCoeffs[coeffI]*intCoeffs[coeffI]/ + // Note consistent sign for boundary coeffs + 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_.initMatrixInterfaces - ( - coupleBouCoeffs_, - interfaces_, - x, - x, // put result into x - cmpt, - false - ); + // matrix_.updateMatrixInterfaces + // ( + // 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_; + // // 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; - - // Coupled boundary update + // Diagonal block + forAll (x, i) { - matrix_.initMatrixInterfaces - ( - 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]; + x[i] = b[i]*preconDiag_[i]; } if (matrix_.asymmetric())