From eb2dcad003401d7b7bc7e867622748c24a3300ae Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 19 Aug 2015 09:14:33 +0100 Subject: [PATCH 1/7] Fixed backslash-newline in options --- applications/utilities/surface/surfaceAutoPatch/Make/options | 4 ++-- applications/utilities/surface/surfaceClean/Make/options | 2 +- .../utilities/surface/surfaceFeatureConvert/Make/options | 2 +- .../utilities/surface/surfaceFeatureExtract/Make/options | 4 ++-- .../utilities/surface/surfaceMeshConvert/Make/options | 4 +++- .../utilities/surface/surfaceMeshExport/Make/options | 4 +++- .../utilities/surface/surfaceMeshImport/Make/options | 4 +++- .../utilities/surface/surfaceMeshTriangulate/Make/options | 5 ++--- applications/utilities/surface/surfaceOrient/Make/options | 2 +- .../utilities/surface/surfacePointMerge/Make/options | 2 +- .../utilities/surface/surfaceRedistributePar/Make/options | 2 +- .../utilities/surface/surfaceRefineRedGreen/Make/options | 4 ++-- .../utilities/surface/surfaceSplitByPatch/Make/options | 1 - .../utilities/surface/surfaceSplitNonManifolds/Make/options | 4 ++-- 14 files changed, 24 insertions(+), 20 deletions(-) diff --git a/applications/utilities/surface/surfaceAutoPatch/Make/options b/applications/utilities/surface/surfaceAutoPatch/Make/options index b83e5e8ab..54c035b8f 100644 --- a/applications/utilities/surface/surfaceAutoPatch/Make/options +++ b/applications/utilities/surface/surfaceAutoPatch/Make/options @@ -1,5 +1,5 @@ EXE_INC = \ - -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude EXE_LIBS = \ - -lmeshTools \ + -lmeshTools diff --git a/applications/utilities/surface/surfaceClean/Make/options b/applications/utilities/surface/surfaceClean/Make/options index 18ff94699..54c035b8f 100644 --- a/applications/utilities/surface/surfaceClean/Make/options +++ b/applications/utilities/surface/surfaceClean/Make/options @@ -2,4 +2,4 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude EXE_LIBS = \ - -lmeshTools \ + -lmeshTools diff --git a/applications/utilities/surface/surfaceFeatureConvert/Make/options b/applications/utilities/surface/surfaceFeatureConvert/Make/options index c2415625a..7418ca2b9 100644 --- a/applications/utilities/surface/surfaceFeatureConvert/Make/options +++ b/applications/utilities/surface/surfaceFeatureConvert/Make/options @@ -4,4 +4,4 @@ EXE_INC = \ EXE_LIBS = \ -lmeshTools \ - -ledgeMesh \ + -ledgeMesh diff --git a/applications/utilities/surface/surfaceFeatureExtract/Make/options b/applications/utilities/surface/surfaceFeatureExtract/Make/options index 45ebbd803..a10600ec1 100644 --- a/applications/utilities/surface/surfaceFeatureExtract/Make/options +++ b/applications/utilities/surface/surfaceFeatureExtract/Make/options @@ -1,6 +1,6 @@ EXE_INC = \ -I$(LIB_SRC)/cfdTools/general/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude EXE_LIBS = \ - -lmeshTools \ + -lmeshTools diff --git a/applications/utilities/surface/surfaceMeshConvert/Make/options b/applications/utilities/surface/surfaceMeshConvert/Make/options index 42b30c865..02c293cee 100644 --- a/applications/utilities/surface/surfaceMeshConvert/Make/options +++ b/applications/utilities/surface/surfaceMeshConvert/Make/options @@ -2,4 +2,6 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/surfMesh/lnInclude -EXE_LIBS = -lmeshTools -lsurfMesh +EXE_LIBS = \ + -lmeshTools \ + -lsurfMesh diff --git a/applications/utilities/surface/surfaceMeshExport/Make/options b/applications/utilities/surface/surfaceMeshExport/Make/options index 42b30c865..02c293cee 100644 --- a/applications/utilities/surface/surfaceMeshExport/Make/options +++ b/applications/utilities/surface/surfaceMeshExport/Make/options @@ -2,4 +2,6 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/surfMesh/lnInclude -EXE_LIBS = -lmeshTools -lsurfMesh +EXE_LIBS = \ + -lmeshTools \ + -lsurfMesh diff --git a/applications/utilities/surface/surfaceMeshImport/Make/options b/applications/utilities/surface/surfaceMeshImport/Make/options index 42b30c865..02c293cee 100644 --- a/applications/utilities/surface/surfaceMeshImport/Make/options +++ b/applications/utilities/surface/surfaceMeshImport/Make/options @@ -2,4 +2,6 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/surfMesh/lnInclude -EXE_LIBS = -lmeshTools -lsurfMesh +EXE_LIBS = \ + -lmeshTools \ + -lsurfMesh diff --git a/applications/utilities/surface/surfaceMeshTriangulate/Make/options b/applications/utilities/surface/surfaceMeshTriangulate/Make/options index 990b27070..54c035b8f 100644 --- a/applications/utilities/surface/surfaceMeshTriangulate/Make/options +++ b/applications/utilities/surface/surfaceMeshTriangulate/Make/options @@ -1,6 +1,5 @@ EXE_INC = \ - -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude EXE_LIBS = \ - -lmeshTools \ - + -lmeshTools diff --git a/applications/utilities/surface/surfaceOrient/Make/options b/applications/utilities/surface/surfaceOrient/Make/options index 18ff94699..54c035b8f 100644 --- a/applications/utilities/surface/surfaceOrient/Make/options +++ b/applications/utilities/surface/surfaceOrient/Make/options @@ -2,4 +2,4 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude EXE_LIBS = \ - -lmeshTools \ + -lmeshTools diff --git a/applications/utilities/surface/surfacePointMerge/Make/options b/applications/utilities/surface/surfacePointMerge/Make/options index b83e5e8ab..f63a1a87c 100644 --- a/applications/utilities/surface/surfacePointMerge/Make/options +++ b/applications/utilities/surface/surfacePointMerge/Make/options @@ -2,4 +2,4 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude \ EXE_LIBS = \ - -lmeshTools \ + -lmeshTools diff --git a/applications/utilities/surface/surfaceRedistributePar/Make/options b/applications/utilities/surface/surfaceRedistributePar/Make/options index 18ff94699..54c035b8f 100644 --- a/applications/utilities/surface/surfaceRedistributePar/Make/options +++ b/applications/utilities/surface/surfaceRedistributePar/Make/options @@ -2,4 +2,4 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude EXE_LIBS = \ - -lmeshTools \ + -lmeshTools diff --git a/applications/utilities/surface/surfaceRefineRedGreen/Make/options b/applications/utilities/surface/surfaceRefineRedGreen/Make/options index b83e5e8ab..54c035b8f 100644 --- a/applications/utilities/surface/surfaceRefineRedGreen/Make/options +++ b/applications/utilities/surface/surfaceRefineRedGreen/Make/options @@ -1,5 +1,5 @@ EXE_INC = \ - -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude EXE_LIBS = \ - -lmeshTools \ + -lmeshTools diff --git a/applications/utilities/surface/surfaceSplitByPatch/Make/options b/applications/utilities/surface/surfaceSplitByPatch/Make/options index 8940e1ec4..54c035b8f 100644 --- a/applications/utilities/surface/surfaceSplitByPatch/Make/options +++ b/applications/utilities/surface/surfaceSplitByPatch/Make/options @@ -3,4 +3,3 @@ EXE_INC = \ EXE_LIBS = \ -lmeshTools - diff --git a/applications/utilities/surface/surfaceSplitNonManifolds/Make/options b/applications/utilities/surface/surfaceSplitNonManifolds/Make/options index b83e5e8ab..54c035b8f 100644 --- a/applications/utilities/surface/surfaceSplitNonManifolds/Make/options +++ b/applications/utilities/surface/surfaceSplitNonManifolds/Make/options @@ -1,5 +1,5 @@ EXE_INC = \ - -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude EXE_LIBS = \ - -lmeshTools \ + -lmeshTools From bf0316a4a741df3302265d712873a2744eb8d93c Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 20 Aug 2015 13:53:03 +0100 Subject: [PATCH 2/7] Use Householder inverse for reconstruct --- src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C b/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C index bc6b8dbd9..435d8adea 100644 --- a/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C +++ b/src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C @@ -68,8 +68,8 @@ reconstruct IOobject::NO_READ, IOobject::NO_WRITE ), - inv(surfaceSum(sqr(mesh.Sf())/mesh.magSf())) - & surfaceSum((mesh.Sf()/mesh.magSf())*ssf), + mesh, + ssf.dimensions()/dimArea, zeroGradientFvPatchField::typeName ) ); @@ -87,9 +87,14 @@ reconstruct GeometricField fluxTimesNormal = surfaceSum((mesh.Sf()/mesh.magSf())*ssf); + // Note: hinv inverse must be used to stabilise the inverse on bad meshes + // HJ, 19/Aug/2015 reconField.internalField() = ( - inv(surfaceSum(sqr(mesh.Sf())/mesh.magSf())().internalField()) + hinv + ( + surfaceSum(sqr(mesh.Sf())/mesh.magSf())().internalField() + ) & fluxTimesNormal.internalField() ); From 8dca6a3763dafc0ebb633fc75b5a82f1675e144c Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 20 Aug 2015 21:11:23 +0100 Subject: [PATCH 3/7] Formatting --- src/foam/primitives/VectorN/DiagTensorNI.H | 278 +++++++++++++++--- .../primitives/VectorN/SphericalTensorNI.H | 96 ++++-- src/foam/primitives/VectorN/TensorNI.H | 88 +++--- .../VectorN/expandContract/ExpandTensorN.H | 6 +- 4 files changed, 357 insertions(+), 111 deletions(-) diff --git a/src/foam/primitives/VectorN/DiagTensorNI.H b/src/foam/primitives/VectorN/DiagTensorNI.H index f521b5557..53d2864d7 100644 --- a/src/foam/primitives/VectorN/DiagTensorNI.H +++ b/src/foam/primitives/VectorN/DiagTensorNI.H @@ -40,6 +40,7 @@ const DiagTensorN DiagTensorN::zero(0); template const DiagTensorN DiagTensorN::one(1); + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // Construct null @@ -63,7 +64,12 @@ inline DiagTensorN::DiagTensorN template inline DiagTensorN::DiagTensorN(const Cmpt& tx) { - VectorSpaceOps::nComponents,0>::eqOpS(*this, tx, eqOp()); + VectorSpaceOps::nComponents,0>::eqOpS + ( + *this, + tx, + eqOp() + ); } @@ -94,80 +100,154 @@ inline DiagTensorN DiagTensorN::T() const //- Assign to a SphericalTensorN template -inline void DiagTensorN::operator=(const SphericalTensorN& st) +inline void DiagTensorN::operator= +( + const SphericalTensorN& st +) { const Cmpt& s = st.v_[0]; - VectorSpaceOps::nComponents,0>::eqOpS(*this, s, eqOp()); + VectorSpaceOps::nComponents,0>::eqOpS + ( + *this, + s, + eqOp() + ); } // * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * // //- Addition of DiagTensorN and DiagTensorN template -inline DiagTensorN -operator+(const DiagTensorN& dt1, const DiagTensorN& dt2) +inline DiagTensorN +operator+ +( + const DiagTensorN& dt1, + const DiagTensorN& dt2 +) { DiagTensorN res; - VectorSpaceOps::nComponents,0>::op(res, dt1, dt2, plusOp()); + VectorSpaceOps::nComponents,0>::op + ( + res, + dt1, + dt2, + plusOp() + ); + return res; } //- Addition of DiagTensorN and SphericalTensorN template -inline DiagTensorN -operator+(const DiagTensorN& dt1, const SphericalTensorN& st2) +inline DiagTensorN +operator+ +( + const DiagTensorN& dt1, + const SphericalTensorN& st2 +) { const Cmpt& s = st2.v_[0]; DiagTensorN res; - VectorSpaceOps::nComponents,0>::opVS(res, dt1, s, plusOp()); + VectorSpaceOps::nComponents,0>::opVS + ( + res, + dt1, + s, + plusOp() + ); + return res; } //- Addition of SphericalTensorN and DiagTensorN template -inline DiagTensorN -operator+(const SphericalTensorN& st1, const DiagTensorN& dt2) +inline DiagTensorN +operator+ +( + const SphericalTensorN& st1, + const DiagTensorN& dt2 +) { const Cmpt& s = st1.v_[0]; DiagTensorN res; - VectorSpaceOps::nComponents,0>::opSV(res, s, dt2, plusOp()); + VectorSpaceOps::nComponents,0>::opSV + ( + res, + s, + dt2, + plusOp() + ); + return res; } //- Subtraction of DiagTensorN and DiagTensorN template -inline DiagTensorN -operator-(const DiagTensorN& dt1, const DiagTensorN& dt2) +inline DiagTensorN +operator- +( + const DiagTensorN& dt1, + const DiagTensorN& dt2 +) { DiagTensorN res; - VectorSpaceOps::nComponents,0>::op(res, dt1, dt2, minusOp()); + VectorSpaceOps::nComponents,0>::op + ( + res, + dt1, + dt2, + minusOp() + ); + return res; } //- Subtraction of DiagTensorN and SphericalTensorN template -inline DiagTensorN -operator-(const DiagTensorN& dt1, const SphericalTensorN& st2) +inline DiagTensorN +operator- +( + const DiagTensorN& dt1, + const SphericalTensorN& st2 +) { const Cmpt& s = st2.v_[0]; DiagTensorN res; - VectorSpaceOps::nComponents,0>::opVS(res, dt1, s, minusOp()); + VectorSpaceOps::nComponents,0>::opVS + ( + res, + dt1, + s, + minusOp() + ); + return res; } //- Subtraction of SphericalTensorN and DiagTensorN template -inline DiagTensorN -operator-(const SphericalTensorN& st1, const DiagTensorN& dt2) +inline DiagTensorN +operator- +( + const SphericalTensorN& st1, + const DiagTensorN& dt2 +) { const Cmpt& s = st1.v_[0]; DiagTensorN res; - VectorSpaceOps::nComponents,0>::opSV(res, s, dt2, minusOp()); + VectorSpaceOps::nComponents,0>::opSV + ( + res, + s, + dt2, + minusOp() + ); + return res; } @@ -176,10 +256,21 @@ operator-(const SphericalTensorN& st1, const DiagTensorN inline typename innerProduct, DiagTensorN >::type -operator&(const DiagTensorN& dt1, const DiagTensorN& dt2) +operator& +( + const DiagTensorN& dt1, + const DiagTensorN& dt2 +) { DiagTensorN res; - VectorSpaceOps::nComponents,0>::op(res, dt1, dt2, multiplyOp()); + VectorSpaceOps::nComponents,0>::op + ( + res, + dt1, + dt2, + multiplyOp() + ); + return res; } @@ -188,11 +279,22 @@ operator&(const DiagTensorN& dt1, const DiagTensorN& template inline typename innerProduct, DiagTensorN >::type -operator&(const SphericalTensorN& st1, const DiagTensorN& dt2) +operator& +( + const SphericalTensorN& st1, + const DiagTensorN& dt2 +) { const Cmpt& s = st1.v_[0]; DiagTensorN res; - VectorSpaceOps::nComponents,0>::opSV(res, s, dt2, multiplyOp()); + VectorSpaceOps::nComponents,0>::opSV + ( + res, + s, + dt2, + multiplyOp() + ); + return res; } @@ -201,11 +303,22 @@ operator&(const SphericalTensorN& st1, const DiagTensorN inline typename innerProduct, SphericalTensorN >::type -operator&(const DiagTensorN& dt1, const SphericalTensorN& st2) +operator& +( + const DiagTensorN& dt1, + const SphericalTensorN& st2 +) { const Cmpt& s = st2.v_[0]; DiagTensorN res; - VectorSpaceOps::nComponents,0>::opVS(res, dt1, s, multiplyOp()); + VectorSpaceOps::nComponents,0>::opVS + ( + res, + dt1, + s, + multiplyOp() + ); + return res; } @@ -214,10 +327,21 @@ operator&(const DiagTensorN& dt1, const SphericalTensorN inline typename innerProduct, VectorN >::type -operator&(const DiagTensorN& dt, const VectorN& v) +operator& +( + const DiagTensorN& dt, + const VectorN& v +) { VectorN res; - VectorSpaceOps::nComponents,0>::opVV(res, dt, v, multiplyOp()); + VectorSpaceOps::nComponents,0>::opVV + ( + res, + dt, + v, + multiplyOp() + ); + return res; } @@ -227,10 +351,21 @@ operator&(const DiagTensorN& dt, const VectorN& v) template inline typename innerProduct, DiagTensorN >::type -operator&(const VectorN& v, const DiagTensorN& dt) +operator& +( + const VectorN& v, + const DiagTensorN& dt +) { VectorN res; - VectorSpaceOps::nComponents,0>::opVV(res, v, dt, multiplyOp()); + VectorSpaceOps::nComponents,0>::opVV + ( + res, + v, + dt, + multiplyOp() + ); + return res; } @@ -241,53 +376,99 @@ inline DiagTensorN operator/(const scalar s, const DiagTensorN& dt) { DiagTensorN res; - VectorSpaceOps::nComponents,0>::opSV(res, s, dt, divideOp3()); + VectorSpaceOps::nComponents,0>::opSV + ( + res, + s, + dt, + divideOp3() + ); + return res; } //- Inner Product of a VectorN by an inverse diagonalTensor template -inline VectorN -operator/(const VectorN& v, const DiagTensorN& dt) +inline VectorN +operator/(const VectorN& v, const DiagTensorN& dt) { VectorN res(v); - VectorSpaceOps::nComponents,0>::eqOp(res, dt, divideEqOp()); + VectorSpaceOps::nComponents,0>::eqOp + ( + res, + dt, + divideEqOp() + ); + return res; } //- Inner Product of a DiagTensorN and an inverse DiagTensorN template -inline DiagTensorN -operator/(const DiagTensorN& dt1, const DiagTensorN& dt2) +inline DiagTensorN +operator/ +( + const DiagTensorN& dt1, + const DiagTensorN& dt2 +) { DiagTensorN res; - VectorSpaceOps::nComponents,0>::op(res, dt1, dt2, divideOp()); + VectorSpaceOps::nComponents,0>::op + ( + res, + dt1, + dt2, + divideOp() + ); + return res; } //- Inner Product of a SphericalTensorN and an inverse DiagTensorN template -inline DiagTensorN -operator/(const SphericalTensorN& st1, const DiagTensorN& dt2) +inline DiagTensorN +operator/ +( + const SphericalTensorN& st1, + const DiagTensorN& dt2 +) { const Cmpt& s = st1.v_[0]; DiagTensorN res; - VectorSpaceOps::nComponents,0>::opSV(res, s, dt2, divideOp()); + VectorSpaceOps::nComponents,0>::opSV + ( + res, + s, + dt2, + divideOp() + ); + return res; } //- Inner Product of a DiagTensorN and an inverse SphericalTensorN template -inline DiagTensorN -operator/(const DiagTensorN& dt1, const SphericalTensorN& st2) +inline DiagTensorN +operator/ +( + const DiagTensorN& dt1, + const SphericalTensorN& st2 +) { const Cmpt& s = st2.v_[0]; DiagTensorN res; - VectorSpaceOps::nComponents,0>::opVS(res, dt1, s, divideOp()); + VectorSpaceOps::nComponents,0>::opVS + ( + res, + dt1, + s, + divideOp() + ); + return res; } @@ -297,7 +478,14 @@ template inline DiagTensorN inv(const DiagTensorN& dt) { DiagTensorN res; - VectorSpaceOps::nComponents,0>::opSV(res, 1.0, dt, divideOp()); + VectorSpaceOps::nComponents,0>::opSV + ( + res, + 1.0, + dt, + divideOp() + ); + return res; } diff --git a/src/foam/primitives/VectorN/SphericalTensorNI.H b/src/foam/primitives/VectorN/SphericalTensorNI.H index 13d27595c..8dd7bd2d1 100644 --- a/src/foam/primitives/VectorN/SphericalTensorNI.H +++ b/src/foam/primitives/VectorN/SphericalTensorNI.H @@ -79,14 +79,16 @@ inline SphericalTensorN::SphericalTensorN(Istream& is) //- Return diagonal tensor diagonal template -inline SphericalTensorN SphericalTensorN::diag() const +inline SphericalTensorN +SphericalTensorN::diag() const { return *this; } //- Return spherical tensor transpose template -inline SphericalTensorN SphericalTensorN::T() const +inline SphericalTensorN +SphericalTensorN::T() const { return *this; } @@ -111,8 +113,12 @@ inline SphericalTensorN transform //- Addition of SphericalTensorN and SphericalTensorN template -inline SphericalTensorN -operator+(const SphericalTensorN& st1, const SphericalTensorN& st2) +inline SphericalTensorN +operator+ +( + const SphericalTensorN& st1, + const SphericalTensorN& st2 +) { return SphericalTensorN(st1.v_[0] + st2.v_[0]); } @@ -120,8 +126,12 @@ operator+(const SphericalTensorN& st1, const SphericalTensorN -inline SphericalTensorN -operator-(const SphericalTensorN& st1, const SphericalTensorN& st2) +inline SphericalTensorN +operator- +( + const SphericalTensorN& st1, + const SphericalTensorN& st2 +) { return SphericalTensorN(st1.v_[0] - st2.v_[0]); } @@ -130,8 +140,16 @@ operator-(const SphericalTensorN& st1, const SphericalTensorN inline typename -innerProduct, SphericalTensorN >::type -operator&(const SphericalTensorN& st1, const SphericalTensorN& st2) +innerProduct +< + SphericalTensorN, + SphericalTensorN +>::type +operator& +( + const SphericalTensorN& st1, + const SphericalTensorN& st2 +) { return SphericalTensorN(st1.v_[0]*st2.v_[0]); } @@ -141,11 +159,22 @@ operator&(const SphericalTensorN& st1, const SphericalTensorN inline typename innerProduct, VectorN >::type -operator&(const SphericalTensorN& st, const VectorN& v) +operator& +( + const SphericalTensorN& st, + const VectorN& v +) { const Cmpt& s = st.v_[0]; VectorN res; - VectorSpaceOps::nComponents,0>::opSV(res, s, v, multiplyOp()); + VectorSpaceOps::nComponents,0>::opSV + ( + res, + s, + v, + multiplyOp() + ); + return res; } @@ -154,11 +183,22 @@ operator&(const SphericalTensorN& st, const VectorN& template inline typename innerProduct, SphericalTensorN >::type -operator&(const VectorN& v, const SphericalTensorN& st) +operator& +( + const VectorN& v, + const SphericalTensorN& st +) { const Cmpt& s = st.v_[0]; VectorN res; - VectorSpaceOps::nComponents,0>::opVS(res, v, s, multiplyOp()); + VectorSpaceOps::nComponents,0>::opVS + ( + res, + v, + s, + multiplyOp() + ); + return res; } @@ -183,8 +223,12 @@ operator/(const scalar s, const SphericalTensorN& st) //- Inner Product of a VectorN by an inverse SphericalTensorN template -inline VectorN -operator/(const VectorN& v, const SphericalTensorN& st) +inline VectorN +operator/ +( + const VectorN& v, + const SphericalTensorN& st +) { return v/st.v_[0]; } @@ -192,8 +236,12 @@ operator/(const VectorN& v, const SphericalTensorN& st //- Inner Product of a SphericalTensorN and an inverse SphericalTensorN template -inline SphericalTensorN -operator/(const SphericalTensorN& st1, const SphericalTensorN& st2) +inline SphericalTensorN +operator/ +( + const SphericalTensorN& st1, + const SphericalTensorN& st2 +) { return SphericalTensorN(st1.v_[0]/st2.v_[0]); } @@ -201,7 +249,10 @@ operator/(const SphericalTensorN& st1, const SphericalTensorN -inline SphericalTensorN inv(const SphericalTensorN& st) +inline SphericalTensorN inv +( + const SphericalTensorN& st +) { return SphericalTensorN(pTraits::one/st.v_[0]); } @@ -209,7 +260,10 @@ inline SphericalTensorN inv(const SphericalTensorN& //- Return tensor diagonal template -inline SphericalTensorN diag(const SphericalTensorN& st) +inline SphericalTensorN diag +( + const SphericalTensorN& st +) { return st; } @@ -258,7 +312,11 @@ public: }; template -class innerProduct, SphericalTensorN > +class innerProduct +< + SphericalTensorN, + SphericalTensorN +> { public: diff --git a/src/foam/primitives/VectorN/TensorNI.H b/src/foam/primitives/VectorN/TensorNI.H index 3514c12c0..764024d1c 100644 --- a/src/foam/primitives/VectorN/TensorNI.H +++ b/src/foam/primitives/VectorN/TensorNI.H @@ -483,11 +483,11 @@ operator*(const VectorN& v1, const VectorN& v2) //- Addition of TensorN and TensorN template -inline TensorN -operator+(const TensorN& t1, const TensorN& t2) +inline TensorN +operator+(const TensorN& t1, const TensorN& t2) { - TensorN res; - VectorSpaceOps::nComponents,0>::op + TensorN res; + VectorSpaceOps::nComponents,0>::op ( res, t1, @@ -501,8 +501,8 @@ operator+(const TensorN& t1, const TensorN& t2) //- Addition of TensorN and DiagTensorN template -inline TensorN -operator+(const TensorN& t1, const DiagTensorN& dt2) +inline TensorN +operator+(const TensorN& t1, const DiagTensorN& dt2) { TensorN result(t1); @@ -519,8 +519,8 @@ operator+(const TensorN& t1, const DiagTensorN& dt2) //- Addition of DiagTensorN and TensorN template -inline TensorN -operator+(const DiagTensorN& dt1, const TensorN& t2) +inline TensorN +operator+(const DiagTensorN& dt1, const TensorN& t2) { TensorN result(t2); @@ -537,11 +537,11 @@ operator+(const DiagTensorN& dt1, const TensorN& t2) //- Addition of TensorN and SphericalTensorN template -inline TensorN +inline TensorN operator+ ( - const TensorN& t1, - const SphericalTensorN& st2 + const TensorN& t1, + const SphericalTensorN& st2 ) { TensorN result(t1); @@ -560,11 +560,11 @@ operator+ //- Addition of SphericalTensorN and TensorN template -inline TensorN +inline TensorN operator+ ( - const SphericalTensorN& st1, - const TensorN& t2 + const SphericalTensorN& st1, + const TensorN& t2 ) { TensorN result(t2); @@ -583,11 +583,11 @@ operator+ //- Subtraction of TensorN and TensorN template -inline TensorN -operator-(const TensorN& t1, const TensorN& t2) +inline TensorN +operator-(const TensorN& t1, const TensorN& t2) { - TensorN res; - VectorSpaceOps::nComponents,0>::op + TensorN res; + VectorSpaceOps::nComponents,0>::op ( res, t1, @@ -601,8 +601,8 @@ operator-(const TensorN& t1, const TensorN& t2) //- Subtraction of TensorN and DiagTensorN template -inline TensorN -operator-(const TensorN& t1, const DiagTensorN& dt2) +inline TensorN +operator-(const TensorN& t1, const DiagTensorN& dt2) { TensorN result(t1); @@ -619,8 +619,8 @@ operator-(const TensorN& t1, const DiagTensorN& dt2) //- Subtraction of DiagTensorN and TensorN template -inline TensorN -operator-(const DiagTensorN& dt1, const TensorN& t2) +inline TensorN +operator-(const DiagTensorN& dt1, const TensorN& t2) { TensorN result(-t2); @@ -637,11 +637,11 @@ operator-(const DiagTensorN& dt1, const TensorN& t2) //- Subtraction of TensorN and SphericalTensorN template -inline TensorN +inline TensorN operator- ( - const TensorN& t1, - const SphericalTensorN& st2 + const TensorN& t1, + const SphericalTensorN& st2 ) { TensorN result(t1); @@ -660,11 +660,11 @@ operator- //- Subtraction of SphericalTensorN and TensorN template -inline TensorN +inline TensorN operator- ( - const SphericalTensorN& st1, - const TensorN& t2 + const SphericalTensorN& st1, + const TensorN& t2 ) { TensorN result(-t2); @@ -683,24 +683,24 @@ operator- //- Division of a scalar by a TensorN template -inline TensorN -operator/(const scalar s, const TensorN& t) +inline TensorN +operator/(const scalar s, const TensorN& t) { return s*inv(t); } //- Inner Product of a VectorN by an inverse TensorN template -inline VectorN -operator/(const VectorN& v, const TensorN& t) +inline VectorN +operator/(const VectorN& v, const TensorN& t) { return v & inv(t); } //- Inner Product of a TensorN by an inverse TensorN template -inline TensorN -operator/(const TensorN& t1, const TensorN& t2) +inline TensorN +operator/(const TensorN& t1, const TensorN& t2) { return t1 & inv(t2); } @@ -708,8 +708,8 @@ operator/(const TensorN& t1, const TensorN& t2) //- Inner Product of a DiagTensorN and an inverse TensorN template -inline TensorN -operator/(const DiagTensorN& dt1, const TensorN& t2) +inline TensorN +operator/(const DiagTensorN& dt1, const TensorN& t2) { return dt1 & inv(t2); } @@ -717,8 +717,8 @@ operator/(const DiagTensorN& dt1, const TensorN& t2) //- Inner Product of a TensorN and an inverse DiagTensorN template -inline TensorN -operator/(const TensorN& t1, const DiagTensorN& dt2) +inline TensorN +operator/(const TensorN& t1, const DiagTensorN& dt2) { return t1 & inv(dt2); } @@ -726,11 +726,11 @@ operator/(const TensorN& t1, const DiagTensorN& dt2) //- Inner Product of a SphericalTensorN and an inverse TensorN template -inline TensorN +inline TensorN operator/ ( - const SphericalTensorN& st1, - const TensorN& t2 + const SphericalTensorN& st1, + const TensorN& t2 ) { return st1.v_[0] * inv(t2); @@ -739,11 +739,11 @@ operator/ //- Inner Product of a TensorN and an inverse SphericalTensorN template -inline TensorN +inline TensorN operator/ ( - const TensorN& t1, - const SphericalTensorN& st2 + const TensorN& t1, + const SphericalTensorN& st2 ) { TensorN result; diff --git a/src/foam/primitives/VectorN/expandContract/ExpandTensorN.H b/src/foam/primitives/VectorN/expandContract/ExpandTensorN.H index 71f50f570..d78304882 100644 --- a/src/foam/primitives/VectorN/expandContract/ExpandTensorN.H +++ b/src/foam/primitives/VectorN/expandContract/ExpandTensorN.H @@ -40,15 +40,15 @@ Author #define UNARY_TEMPLATE_FUNCTION(ReturnType, Type, Func) \ template \ -inline void Func(ReturnType&, const Type&); +inline void Func(ReturnType&, const Type&); #define UNARY_TEMPLATE_FUNCTION_VS(ReturnType, Func) \ template \ -inline void Func(ReturnType&, const Cmpt&); +inline void Func(ReturnType&, const Cmpt&); #define UNARY_TEMPLATE_FUNCTION_SV(Type, Func) \ template \ -inline void Func(Cmpt&, const Type&); +inline void Func(Cmpt&, const Type&); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // From bb9cbfbbec2cf9dc27e5f46d87ee5299cdc75bc7 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 20 Aug 2015 21:12:00 +0100 Subject: [PATCH 4/7] Formatting --- .../fields/CoeffField/CoeffFieldFunctions.C | 24 +++++++++---------- .../CoeffField/DecoupledCoeffFieldFunctions.H | 5 +++- 2 files changed, 16 insertions(+), 13 deletions(-) diff --git a/src/foam/fields/CoeffField/CoeffFieldFunctions.C b/src/foam/fields/CoeffField/CoeffFieldFunctions.C index 1baf1bb96..b25bb602c 100644 --- a/src/foam/fields/CoeffField/CoeffFieldFunctions.C +++ b/src/foam/fields/CoeffField/CoeffFieldFunctions.C @@ -143,7 +143,7 @@ void multiply #define UNARY_OPERATOR(op, opFunc) \ \ template \ -void opFunc \ +void opFunc \ ( \ CoeffField& f, \ const CoeffField& f1 \ @@ -176,7 +176,7 @@ void opFunc \ } \ \ template \ -tmp > operator op \ +tmp > operator op \ ( \ const CoeffField& f1 \ ) \ @@ -187,7 +187,7 @@ tmp > operator op \ } \ \ template \ -tmp > operator op \ +tmp > operator op \ ( \ const tmp >& tf1 \ ) \ @@ -205,7 +205,7 @@ UNARY_OPERATOR(-, negate) #define BINARY_OPERATOR_FF(Type1, Type2, op, opFunc) \ \ template \ -tmp > operator op \ +tmp > operator op \ ( \ const CoeffField& f1, \ const Type2& f2 \ @@ -218,7 +218,7 @@ tmp > operator op \ \ \ template \ -tmp > operator op \ +tmp > operator op \ ( \ const CoeffField& f1, \ const Field& f2 \ @@ -231,7 +231,7 @@ tmp > operator op \ \ \ template \ -tmp > operator op \ +tmp > operator op \ ( \ const Field& f1, \ const CoeffField& f2 \ @@ -244,7 +244,7 @@ tmp > operator op \ #define BINARY_OPERATOR_FTR(Type1, Type2, op, opFunc) \ template \ -tmp > operator op \ +tmp > operator op \ ( \ const CoeffField& f1, \ const tmp >& tf2 \ @@ -257,7 +257,7 @@ tmp > operator op \ #define BINARY_OPERATOR_FT(Type1, Type2, op, opFunc) \ template \ -tmp > operator op \ +tmp > operator op \ ( \ const Field& f1, \ const tmp >& tf2 \ @@ -270,7 +270,7 @@ tmp > operator op \ #define BINARY_OPERATOR_TRF(Type1, Type2, op, opFunc) \ template \ -tmp > operator op \ +tmp > operator op \ ( \ const tmp >& tf1, \ const Field& f2 \ @@ -283,7 +283,7 @@ tmp > operator op \ #define BINARY_OPERATOR_TF(Type1, Type2, op, opFunc) \ template \ -tmp > operator op \ +tmp > operator op \ ( \ const tmp >& tf1, \ const Field& f2 \ @@ -296,7 +296,7 @@ tmp > operator op \ #define BINARY_OPERATOR_TRT(Type1, Type2, op, opFunc) \ template \ -tmp > operator op \ +tmp > operator op \ ( \ const tmp >& tf1, \ const tmp >& tf2 \ @@ -310,7 +310,7 @@ tmp > operator op \ #define BINARY_OPERATOR_TTR(Type1, Type2, op, opFunc) \ template \ -tmp > operator op \ +tmp > operator op \ ( \ const tmp >& tf1, \ const tmp >& tf2 \ diff --git a/src/foam/fields/CoeffField/DecoupledCoeffFieldFunctions.H b/src/foam/fields/CoeffField/DecoupledCoeffFieldFunctions.H index 3f749284d..095d66d86 100644 --- a/src/foam/fields/CoeffField/DecoupledCoeffFieldFunctions.H +++ b/src/foam/fields/CoeffField/DecoupledCoeffFieldFunctions.H @@ -37,7 +37,10 @@ namespace Foam /* * * * * * * * * * * * * * * * Global functions * * * * * * * * * * * * * */ template -tmp > inv(const DecoupledCoeffField& f); +tmp > inv +( + const DecoupledCoeffField& f +); template From fb9862668610ec7f990a9d8912f176276d953314 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 20 Aug 2015 21:12:15 +0100 Subject: [PATCH 5/7] Formatting --- .../GGIInterpolation/GGIInterpolationQuickRejectTests.C | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/foam/interpolations/GGIInterpolation/GGIInterpolationQuickRejectTests.C b/src/foam/interpolations/GGIInterpolation/GGIInterpolationQuickRejectTests.C index 8626171cb..f7de56cba 100644 --- a/src/foam/interpolations/GGIInterpolation/GGIInterpolationQuickRejectTests.C +++ b/src/foam/interpolations/GGIInterpolation/GGIInterpolationQuickRejectTests.C @@ -181,7 +181,7 @@ void GGIInterpolation::findNeighbours3D scalar tmpValue = Foam::magSqr(bbSlave.max() - bbSlave.min())/4.0; - // We will compare squared distances, so save the sqrt() if value > 1.0. + // We compare squared distances, so save the sqrt() if value > 1.0. if (tmpValue < 1.0) { // Take the sqrt, otherwise, we underestimate the radius @@ -566,7 +566,6 @@ void GGIInterpolation::findNeighboursBBOctree } - // Projects a list of points onto a plane located at planeOrig, // oriented along planeNormal. Return the projected points in a // pointField, and the normal distance of each points from the From 1685dc456557c56f705bf039b25968f4ef2ffe36 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 21 Aug 2015 08:16:02 +0100 Subject: [PATCH 6/7] Clean-up of block Cholesky precon --- .../BlockCholeskyPrecon/BlockCholeskyPrecon.C | 21 +- .../BlockCholeskyPreconDecoupled.C | 5 + .../BlockGaussSeidelPrecon.C | 270 ++++++++++++++---- .../BlockGaussSeidelPrecon.H | 30 +- .../BlockGaussSeidelPreconDecoupled.C | 67 +++-- .../scalarBlockGaussSeidelPrecon.C | 35 +-- .../scalarBlockGaussSeidelPrecon.H | 20 +- .../tensorBlockGaussSeidelPrecon.C | 38 +-- .../tensorBlockGaussSeidelPrecon.H | 21 +- 9 files changed, 289 insertions(+), 218 deletions(-) diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C index f3f7d4ed8..fb0c2eca5 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C @@ -220,6 +220,9 @@ void Foam::BlockCholeskyPrecon::calcPreconDiag() } } } + + // Invert the diagonal + preconDiag_ = inv(preconDiag_); } @@ -250,12 +253,6 @@ void Foam::BlockCholeskyPrecon::diagMultiply upper[coeffI] ); } - - // Invert the diagonal for future use - forAll (dDiag, i) - { - dDiag[i] = mult.inverse(dDiag[i]); - } } @@ -286,12 +283,6 @@ void Foam::BlockCholeskyPrecon::diagMultiplyCoeffT upper[coeffI] ); } - - // Invert the diagonal for future use - forAll (dDiag, i) - { - dDiag[i] = mult.inverse(dDiag[i]); - } } @@ -323,12 +314,6 @@ void Foam::BlockCholeskyPrecon::diagMultiply upper[coeffI] ); } - - // Invert the diagonal for future use - forAll (dDiag, i) - { - dDiag[i] = mult.inverse(dDiag[i]); - } } diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPreconDecoupled.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPreconDecoupled.C index e426b1d2b..3ca795564 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPreconDecoupled.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPreconDecoupled.C @@ -126,6 +126,11 @@ void Foam::BlockCholeskyPrecon::calcDecoupledPreconDiag() } } } + + // Invert the diagonal + // Note: since square diag type does not exist, simple inversion + // covers all cases + preconDiag_ = inv(preconDiag_); } diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C index bf6d2f4f6..374387af8 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C @@ -34,6 +34,79 @@ Author // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +template +void Foam::BlockGaussSeidelPrecon::calcInvDiag() +{ + typedef CoeffField TypeCoeffField; + typedef typename TypeCoeffField::linearTypeField linearTypeField; + typedef typename TypeCoeffField::linearType valueType; + + typedef typename TypeCoeffField::squareTypeField squareTypeField; + + // Note: Cannot use inv since the square coefficient requires + // special treatment. HJ, 20/Aug/2015 + + // Get reference to diagonal + const TypeCoeffField& d = this->matrix_.diag(); + + if (d.activeType() == blockCoeffBase::SCALAR) + { + invDiag_.asScalar() = 1/d.asScalar(); + } + else if (d.activeType() == blockCoeffBase::LINEAR) + { + invDiag_.asLinear() = + cmptDivide + ( + linearTypeField(d.size(), pTraits::one), + d.asLinear() + ); + } + else if (d.activeType() == blockCoeffBase::SQUARE) + { + // For square diagonal invert diagonal only and store the rest + // info LUDiag coefficient. This avoids full inverse of the + // diagonal matrix. HJ, 20/Aug/2015 + Info<< "Square diag inverse" << endl; + + // Get reference to active diag + const squareTypeField& activeDiag = d.asSquare(); + + // Get reference to LU: remove diagonal from active diag + squareTypeField& luDiag = LUDiag_.asSquare(); + + linearTypeField lf(activeDiag.size()); + + // Take out the diagonal from the diag as a linear type + contractLinear(lf, activeDiag); + + // Expand diagonal only to full square type and store into luDiag + expandLinear(luDiag, lf); + + // Keep only off-diagonal in ldDiag. + // Note change of sign to avoid multiplication with -1 when moving + // to the other side. HJ, 20/Aug/2015 + luDiag -= activeDiag; + + // Inverse is the inverse of diagonal only + invDiag_.asLinear() = + cmptDivide + ( + linearTypeField(lf.size(), pTraits::one), + lf + ); + } + else + { + FatalErrorIn + ( + "void BlockGaussSeidelPrecon::calcInvDiag()" + ) << "Problem with coefficient type morphing." + << abort(FatalError); + } +} + + // Block sweep, symmetric matrix template template @@ -46,7 +119,8 @@ void Foam::BlockGaussSeidelPrecon::BlockSweep ) const { const unallocLabelList& u = this->matrix_.lduAddr().upperAddr(); - const unallocLabelList& ownStart = this->matrix_.lduAddr().ownerStartAddr(); + const unallocLabelList& ownStart = + this->matrix_.lduAddr().ownerStartAddr(); const label nRows = ownStart.size() - 1; @@ -155,7 +229,8 @@ void Foam::BlockGaussSeidelPrecon::BlockSweep ) const { const unallocLabelList& u = this->matrix_.lduAddr().upperAddr(); - const unallocLabelList& ownStart = this->matrix_.lduAddr().ownerStartAddr(); + const unallocLabelList& ownStart = + this->matrix_.lduAddr().ownerStartAddr(); const label nRows = ownStart.size() - 1; @@ -247,6 +322,43 @@ void Foam::BlockGaussSeidelPrecon::BlockSweep } +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::BlockGaussSeidelPrecon::BlockGaussSeidelPrecon +( + const BlockLduMatrix& matrix +) +: + BlockLduPrecon(matrix), + invDiag_(matrix.lduAddr().size()), + LUDiag_(matrix.lduAddr().size()), + bPlusLU_(), + bPrime_(matrix.lduAddr().size()), + nSweeps_(1) +{ + calcInvDiag(); +} + + +template +Foam::BlockGaussSeidelPrecon::BlockGaussSeidelPrecon +( + const BlockLduMatrix& matrix, + const dictionary& dict +) +: + BlockLduPrecon(matrix), + invDiag_(matrix.lduAddr().size()), + LUDiag_(matrix.lduAddr().size()), + bPlusLU_(), + bPrime_(matrix.lduAddr().size()), + nSweeps_(readLabel(dict.lookup("nSweeps"))) +{ + calcInvDiag(); +} + + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template @@ -260,14 +372,11 @@ void Foam::BlockGaussSeidelPrecon::precondition if (this->matrix_.diagonal()) { - TypeCoeffField dDCoeff = inv(this->matrix_.diag()); - - multiply(x, dDCoeff, b); + multiply(x, invDiag_, b); } else if (this->matrix_.symmetric()) { - TypeCoeffField dDCoeff = inv(this->matrix_.diag()); - + const TypeCoeffField& DiagCoeff = this->matrix_.diag(); const TypeCoeffField& UpperCoeff = this->matrix_.upper(); // Note @@ -279,16 +388,16 @@ void Foam::BlockGaussSeidelPrecon::precondition // to be enforced without the per-element if-condition, which // makes for ugly code. HJ, 19/May/2005 - //Note: Assuming lower and upper triangle have the same active type + // Note: Assuming lower and upper triangle have the same active type - if (dDCoeff.activeType() == blockCoeffBase::SCALAR) + if (DiagCoeff.activeType() == blockCoeffBase::SCALAR) { if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) { BlockSweep ( x, - dDCoeff.asScalar(), + invDiag_.asScalar(), UpperCoeff.asScalar(), b ); @@ -298,7 +407,7 @@ void Foam::BlockGaussSeidelPrecon::precondition BlockSweep ( x, - dDCoeff.asScalar(), + invDiag_.asScalar(), UpperCoeff.asLinear(), b ); @@ -308,20 +417,20 @@ void Foam::BlockGaussSeidelPrecon::precondition BlockSweep ( x, - dDCoeff.asScalar(), + invDiag_.asScalar(), UpperCoeff.asSquare(), b ); } } - else if (dDCoeff.activeType() == blockCoeffBase::LINEAR) + else if (DiagCoeff.activeType() == blockCoeffBase::LINEAR) { if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) { BlockSweep ( x, - dDCoeff.asLinear(), + invDiag_.asLinear(), UpperCoeff.asScalar(), b ); @@ -331,7 +440,7 @@ void Foam::BlockGaussSeidelPrecon::precondition BlockSweep ( x, - dDCoeff.asLinear(), + invDiag_.asLinear(), UpperCoeff.asLinear(), b ); @@ -341,42 +450,57 @@ void Foam::BlockGaussSeidelPrecon::precondition BlockSweep ( x, - dDCoeff.asLinear(), + invDiag_.asLinear(), UpperCoeff.asSquare(), b ); } } - else if (dDCoeff.activeType() == blockCoeffBase::SQUARE) + else if (DiagCoeff.activeType() == blockCoeffBase::SQUARE) { + // Add diag coupling to b + if (bPlusLU_.empty()) + { + bPlusLU_.setSize(b.size(), pTraits::zero); + } + + // Multiply overwrites bPlusLU_: no need to initialise + // Change of sign accounted via change of sign of bPlusLU_ + // HJ, 20/Aug/2015 + multiply(bPlusLU_, LUDiag_, x); + bPlusLU_ += b; + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) { + // Note linear diag inversed due to decoupling BlockSweep ( x, - dDCoeff.asSquare(), + invDiag_.asLinear(), UpperCoeff.asScalar(), - b + bPlusLU_ ); } else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) { + // Note linear diag inversed due to decoupling BlockSweep ( x, - dDCoeff.asSquare(), + invDiag_.asLinear(), UpperCoeff.asLinear(), - b + bPlusLU_ ); } else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) { + // Note linear diag inversed due to decoupling BlockSweep ( x, - dDCoeff.asSquare(), + invDiag_.asLinear(), UpperCoeff.asSquare(), - b + bPlusLU_ ); } } @@ -395,8 +519,7 @@ void Foam::BlockGaussSeidelPrecon::precondition } else if (this->matrix_.asymmetric()) { - TypeCoeffField dDCoeff = inv(this->matrix_.diag()); - + const TypeCoeffField& DiagCoeff = this->matrix_.diag(); const TypeCoeffField& LowerCoeff = this->matrix_.lower(); const TypeCoeffField& UpperCoeff = this->matrix_.upper(); @@ -411,14 +534,14 @@ void Foam::BlockGaussSeidelPrecon::precondition //Note: Assuming lower and upper triangle have the same active type - if (dDCoeff.activeType() == blockCoeffBase::SCALAR) + if (DiagCoeff.activeType() == blockCoeffBase::SCALAR) { if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) { BlockSweep ( x, - dDCoeff.asScalar(), + invDiag_.asScalar(), LowerCoeff.asScalar(), UpperCoeff.asScalar(), b @@ -429,7 +552,7 @@ void Foam::BlockGaussSeidelPrecon::precondition BlockSweep ( x, - dDCoeff.asScalar(), + invDiag_.asScalar(), LowerCoeff.asLinear(), UpperCoeff.asLinear(), b @@ -440,21 +563,21 @@ void Foam::BlockGaussSeidelPrecon::precondition BlockSweep ( x, - dDCoeff.asScalar(), + invDiag_.asScalar(), LowerCoeff.asSquare(), UpperCoeff.asSquare(), b ); } } - else if (dDCoeff.activeType() == blockCoeffBase::LINEAR) + else if (DiagCoeff.activeType() == blockCoeffBase::LINEAR) { if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) { BlockSweep ( x, - dDCoeff.asLinear(), + invDiag_.asLinear(), LowerCoeff.asScalar(), UpperCoeff.asScalar(), b @@ -465,7 +588,7 @@ void Foam::BlockGaussSeidelPrecon::precondition BlockSweep ( x, - dDCoeff.asLinear(), + invDiag_.asLinear(), LowerCoeff.asLinear(), UpperCoeff.asLinear(), b @@ -476,46 +599,61 @@ void Foam::BlockGaussSeidelPrecon::precondition BlockSweep ( x, - dDCoeff.asLinear(), + invDiag_.asLinear(), LowerCoeff.asSquare(), UpperCoeff.asSquare(), b ); } } - else if (dDCoeff.activeType() == blockCoeffBase::SQUARE) + else if (DiagCoeff.activeType() == blockCoeffBase::SQUARE) { + // Add diag coupling to b + if (bPlusLU_.empty()) + { + bPlusLU_.setSize(b.size(), pTraits::zero); + } + + // Multiply overwrites bPlusLU_: no need to initialise + // Change of sign accounted via change of sign of bPlusLU_ + // HJ, 20/Aug/2015 + multiply(bPlusLU_, LUDiag_, x); + bPlusLU_ += b; + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) { + // Note linear diag inversed due to decoupling BlockSweep ( x, - dDCoeff.asSquare(), + invDiag_.asLinear(), LowerCoeff.asScalar(), UpperCoeff.asScalar(), - b + bPlusLU_ ); } else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) { + // Note linear diag inversed due to decoupling BlockSweep ( x, - dDCoeff.asSquare(), + invDiag_.asLinear(), LowerCoeff.asLinear(), UpperCoeff.asLinear(), - b + bPlusLU_ ); } else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) { + // Note linear diag inversed due to decoupling BlockSweep ( x, - dDCoeff.asSquare(), + invDiag_.asLinear(), LowerCoeff.asSquare(), UpperCoeff.asSquare(), - b + bPlusLU_ ); } } @@ -558,14 +696,11 @@ void Foam::BlockGaussSeidelPrecon::preconditionT if (this->matrix_.diagonal()) { - TypeCoeffField dDCoeff = inv(this->matrix_.diag()); - - multiply(xT, dDCoeff, bT); + multiply(xT, invDiag_, bT); } else if (this->matrix_.symmetric() || this->matrix_.asymmetric()) { - TypeCoeffField dDCoeff = inv(this->matrix_.diag()); - + const TypeCoeffField& DiagCoeff = this->matrix_.lower(); const TypeCoeffField& LowerCoeff = this->matrix_.lower(); const TypeCoeffField& UpperCoeff = this->matrix_.upper(); @@ -580,7 +715,7 @@ void Foam::BlockGaussSeidelPrecon::preconditionT //Note: Assuming lower and upper triangle have the same active type - if (dDCoeff.activeType() == blockCoeffBase::SCALAR) + if (DiagCoeff.activeType() == blockCoeffBase::SCALAR) { if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) { @@ -588,7 +723,7 @@ void Foam::BlockGaussSeidelPrecon::preconditionT BlockSweep ( xT, - dDCoeff.asScalar(), + invDiag_.asScalar(), UpperCoeff.asScalar(), LowerCoeff.asScalar(), bT @@ -600,7 +735,7 @@ void Foam::BlockGaussSeidelPrecon::preconditionT BlockSweep ( xT, - dDCoeff.asScalar(), + invDiag_.asScalar(), UpperCoeff.asLinear(), LowerCoeff.asLinear(), bT @@ -612,14 +747,14 @@ void Foam::BlockGaussSeidelPrecon::preconditionT BlockSweep ( xT, - dDCoeff.asScalar(), + invDiag_.asScalar(), UpperCoeff.asSquare(), LowerCoeff.asSquare(), bT ); } } - else if (dDCoeff.activeType() == blockCoeffBase::LINEAR) + else if (DiagCoeff.activeType() == blockCoeffBase::LINEAR) { if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) { @@ -627,7 +762,7 @@ void Foam::BlockGaussSeidelPrecon::preconditionT BlockSweep ( xT, - dDCoeff.asLinear(), + invDiag_.asLinear(), UpperCoeff.asScalar(), LowerCoeff.asScalar(), bT @@ -639,7 +774,7 @@ void Foam::BlockGaussSeidelPrecon::preconditionT BlockSweep ( xT, - dDCoeff.asLinear(), + invDiag_.asLinear(), UpperCoeff.asLinear(), LowerCoeff.asLinear(), bT @@ -651,49 +786,64 @@ void Foam::BlockGaussSeidelPrecon::preconditionT BlockSweep ( xT, - dDCoeff.asLinear(), + invDiag_.asLinear(), UpperCoeff.asSquare(), LowerCoeff.asSquare(), bT ); } } - else if (dDCoeff.activeType() == blockCoeffBase::SQUARE) + else if (DiagCoeff.activeType() == blockCoeffBase::SQUARE) { + // Add diag coupling to b + if (bPlusLU_.empty()) + { + bPlusLU_.setSize(bT.size(), pTraits::zero); + } + + // Multiply overwrites bPlusLU_: no need to initialise + // Change of sign accounted via change of sign of bPlusLU_ + // HJ, 20/Aug/2015 + multiply(bPlusLU_, LUDiag_, xT); + bPlusLU_ += bT; + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) { // Transpose multiplication - swap lower and upper coeff arrays + // Note linear diag inversed due to decoupling BlockSweep ( xT, - dDCoeff.asSquare(), + invDiag_.asLinear(), UpperCoeff.asScalar(), LowerCoeff.asScalar(), - bT + bPlusLU_ ); } else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) { // Transpose multiplication - swap lower and upper coeff arrays + // Note linear diag inversed due to decoupling BlockSweep ( xT, - dDCoeff.asSquare(), + invDiag_.asLinear(), UpperCoeff.asLinear(), LowerCoeff.asLinear(), - bT + bPlusLU_ ); } else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) { // Transpose multiplication - swap lower and upper coeff arrays + // Note linear diag inversed due to decoupling BlockSweep ( xT, - dDCoeff.asSquare(), + invDiag_.asLinear(), UpperCoeff.asSquare(), LowerCoeff.asSquare(), - bT + bPlusLU_ ); } } diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H index 8ba166789..1a3efe91a 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H @@ -58,6 +58,16 @@ class BlockGaussSeidelPrecon { // Private Data + //- Inverse of diagonal diagonal + CoeffField invDiag_; + + //- Off-diag part of diagonal + CoeffField LUDiag_; + + //- Temporary space for updated decoupled source + // Initialised with zero size and resized on first use + mutable Field bPlusLU_; + //- Temporary space for solution intermediate mutable Field bPrime_; @@ -74,6 +84,9 @@ class BlockGaussSeidelPrecon void operator=(const BlockGaussSeidelPrecon&); + //- Calculate inverse diagonal + void calcInvDiag(); + // Block Gauss-Seidel sweep, symetric matrix template void BlockSweep @@ -98,6 +111,9 @@ class BlockGaussSeidelPrecon // Decoupled operations, used in template specialisation + //- Calculate inverse diagonal, decoupled version + void calcDecoupledInvDiag(); + //- Execute preconditioning, decoupled version void decoupledPrecondition ( @@ -126,24 +142,14 @@ public: BlockGaussSeidelPrecon ( const BlockLduMatrix& matrix - ) - : - BlockLduPrecon(matrix), - bPrime_(matrix.lduAddr().size()), - nSweeps_(1) - {} + ); //- Construct from components BlockGaussSeidelPrecon ( const BlockLduMatrix& matrix, const dictionary& dict - ) - : - BlockLduPrecon(matrix), - bPrime_(matrix.lduAddr().size()), - nSweeps_(readLabel(dict.lookup("nSweeps"))) - {} + ); // Destructor diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPreconDecoupled.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPreconDecoupled.C index 3b72a19d6..6b8b6cc62 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPreconDecoupled.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPreconDecoupled.C @@ -34,6 +34,19 @@ Author // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +template +void Foam::BlockGaussSeidelPrecon::calcDecoupledInvDiag() +{ + // Get reference to diagonal and obtain inverse by casting + typedef CoeffField TypeCoeffField; + + const TypeCoeffField& d = this->matrix_.diag(); + const DecoupledCoeffField& dd = d; + + invDiag_ = CoeffField(inv(dd)()); +} + + template void Foam::BlockGaussSeidelPrecon::decoupledPrecondition ( @@ -45,14 +58,17 @@ void Foam::BlockGaussSeidelPrecon::decoupledPrecondition if (this->matrix_.diagonal()) { - TypeCoeffField dDCoeff = inv(this->matrix_.diag()); - - multiply(x, dDCoeff, b); + if (invDiag_.activeType() == blockCoeffBase::SCALAR) + { + x = invDiag_.asScalar()*b; + } + else if (invDiag_.activeType() == blockCoeffBase::LINEAR) + { + x = cmptMultiply(invDiag_.asLinear(), b); + } } else if (this->matrix_.symmetric() || this->matrix_.asymmetric()) { - TypeCoeffField dDCoeff = inv(this->matrix_.diag()); - const TypeCoeffField& LowerCoeff = this->matrix_.lower(); const TypeCoeffField& UpperCoeff = this->matrix_.upper(); @@ -65,16 +81,16 @@ void Foam::BlockGaussSeidelPrecon::decoupledPrecondition // to be enforced without the per-element if-condition, which // makes for ugly code. HJ, 19/May/2005 - //Note: Assuming lower and upper triangle have the same active type + // Note: Assuming lower and upper triangle have the same active type - if (dDCoeff.activeType() == blockCoeffBase::SCALAR) + if (invDiag_.activeType() == blockCoeffBase::SCALAR) { if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) { BlockSweep ( x, - dDCoeff.asScalar(), + invDiag_.asScalar(), LowerCoeff.asScalar(), UpperCoeff.asScalar(), b @@ -85,21 +101,21 @@ void Foam::BlockGaussSeidelPrecon::decoupledPrecondition BlockSweep ( x, - dDCoeff.asScalar(), + invDiag_.asScalar(), LowerCoeff.asLinear(), UpperCoeff.asLinear(), b ); } } - else if (dDCoeff.activeType() == blockCoeffBase::LINEAR) + else if (invDiag_.activeType() == blockCoeffBase::LINEAR) { if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) { BlockSweep ( x, - dDCoeff.asLinear(), + invDiag_.asLinear(), LowerCoeff.asScalar(), UpperCoeff.asScalar(), b @@ -110,7 +126,7 @@ void Foam::BlockGaussSeidelPrecon::decoupledPrecondition BlockSweep ( x, - dDCoeff.asLinear(), + invDiag_.asLinear(), LowerCoeff.asLinear(), UpperCoeff.asLinear(), b @@ -156,14 +172,17 @@ void Foam::BlockGaussSeidelPrecon::decoupledPreconditionT if (this->matrix_.diagonal()) { - TypeCoeffField dDCoeff = inv(this->matrix_.diag()); - - multiply(xT, dDCoeff, bT); + if (invDiag_.activeType() == blockCoeffBase::SCALAR) + { + xT = invDiag_.asScalar()*bT; + } + else if (invDiag_.activeType() == blockCoeffBase::LINEAR) + { + xT = cmptMultiply(invDiag_.asLinear(), bT); + } } else if (this->matrix_.symmetric() || this->matrix_.asymmetric()) { - TypeCoeffField dDCoeff = inv(this->matrix_.diag()); - const TypeCoeffField& LowerCoeff = this->matrix_.lower(); const TypeCoeffField& UpperCoeff = this->matrix_.upper(); @@ -176,9 +195,9 @@ void Foam::BlockGaussSeidelPrecon::decoupledPreconditionT // to be enforced without the per-element if-condition, which // makes for ugly code. HJ, 19/May/2005 - //Note: Assuming lower and upper triangle have the same active type + // Note: Assuming lower and upper triangle have the same active type - if (dDCoeff.activeType() == blockCoeffBase::SCALAR) + if (invDiag_.activeType() == blockCoeffBase::SCALAR) { if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) { @@ -186,7 +205,7 @@ void Foam::BlockGaussSeidelPrecon::decoupledPreconditionT BlockSweep ( xT, - dDCoeff.asScalar(), + invDiag_.asScalar(), UpperCoeff.asScalar(), LowerCoeff.asScalar(), bT @@ -198,14 +217,14 @@ void Foam::BlockGaussSeidelPrecon::decoupledPreconditionT BlockSweep ( xT, - dDCoeff.asScalar(), + invDiag_.asScalar(), UpperCoeff.asLinear(), LowerCoeff.asLinear(), bT ); } } - else if (dDCoeff.activeType() == blockCoeffBase::LINEAR) + else if (invDiag_.activeType() == blockCoeffBase::LINEAR) { if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) { @@ -213,7 +232,7 @@ void Foam::BlockGaussSeidelPrecon::decoupledPreconditionT BlockSweep ( xT, - dDCoeff.asLinear(), + invDiag_.asLinear(), UpperCoeff.asScalar(), LowerCoeff.asScalar(), bT @@ -225,7 +244,7 @@ void Foam::BlockGaussSeidelPrecon::decoupledPreconditionT BlockSweep ( xT, - dDCoeff.asLinear(), + invDiag_.asLinear(), UpperCoeff.asLinear(), LowerCoeff.asLinear(), bT diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/scalarBlockGaussSeidelPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/scalarBlockGaussSeidelPrecon.C index 96b2fcf51..fe6d0070a 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/scalarBlockGaussSeidelPrecon.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/scalarBlockGaussSeidelPrecon.C @@ -44,38 +44,11 @@ namespace Foam { template<> -template<> -void BlockGaussSeidelPrecon::BlockSweep -( - Field& x, - const Field& dD, - const Field& upper, - const Field& b -) const +void BlockGaussSeidelPrecon::calcInvDiag() { - FatalErrorIn - ( - "BlockGaussSeidelPrecon::BlockSweep(...)" - ) << "Function not implemented for Type=scalar. " << endl - << abort(FatalError); -} - -template<> -template<> -void BlockGaussSeidelPrecon::BlockSweep -( - Field& x, - const Field& dD, - const Field& upper, - const Field& lower, - const Field& b -) const -{ - FatalErrorIn - ( - "BlockGaussSeidelPrecon::BlockSweep(...)" - ) << "Function not implemented for Type=scalar. " << endl - << abort(FatalError); + // Direct inversion of diagonal is sufficient, as the diagonal + // is linear. HJ, 20/Aug/2015 + invDiag_ = inv(this->matrix_.diag()); } diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/scalarBlockGaussSeidelPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/scalarBlockGaussSeidelPrecon.H index 2f82080c7..8850dbcf7 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/scalarBlockGaussSeidelPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/scalarBlockGaussSeidelPrecon.H @@ -43,25 +43,7 @@ namespace Foam { template<> -template<> -void BlockGaussSeidelPrecon::BlockSweep -( - scalarField& x, - const scalarField& dD, - const scalarField& upper, - const scalarField& b -) const; - -template<> -template<> -void BlockGaussSeidelPrecon::BlockSweep -( - scalarField& x, - const scalarField& dD, - const scalarField& upper, - const scalarField& lower, - const scalarField& b -) const; +void BlockGaussSeidelPrecon::calcInvDiag(); template<> diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/tensorBlockGaussSeidelPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/tensorBlockGaussSeidelPrecon.C index a15bffd64..7c9d62910 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/tensorBlockGaussSeidelPrecon.C +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/tensorBlockGaussSeidelPrecon.C @@ -47,39 +47,9 @@ namespace Foam { template<> -template<> -void BlockGaussSeidelPrecon::BlockSweep -( - Field& x, - const Field& dD, - const Field& upper, - const Field& b -) const +void BlockGaussSeidelPrecon::calcInvDiag() { - FatalErrorIn - ( - "Foam::BlockGaussSeidelPrecon::BlockSweep(...)" - ) << "Function not implemented for Type=tensor. " << endl - << abort(FatalError); -} - - -template<> -template<> -void BlockGaussSeidelPrecon::BlockSweep -( - Field& x, - const Field& dD, - const Field& upper, - const Field& lower, - const Field& b -) const -{ - FatalErrorIn - ( - "Foam::BlockGaussSeidelPrecon::BlockSweep(...)" - ) << "Function not implemented for Type=tensor. " << endl - << abort(FatalError); + calcDecoupledInvDiag(); } @@ -91,7 +61,7 @@ void BlockGaussSeidelPrecon::precondition ) const { // Decoupled version -// decoupledPrecondition(x, b); + decoupledPrecondition(x, b); } @@ -103,7 +73,7 @@ void BlockGaussSeidelPrecon::preconditionT ) const { // Decoupled version -// decoupledPreconditionT(xT, bT); + decoupledPreconditionT(xT, bT); } diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/tensorBlockGaussSeidelPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/tensorBlockGaussSeidelPrecon.H index 574909110..949b9faa6 100644 --- a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/tensorBlockGaussSeidelPrecon.H +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/tensorBlockGaussSeidelPrecon.H @@ -43,26 +43,7 @@ namespace Foam { template<> -template<> -void BlockGaussSeidelPrecon::BlockSweep -( - Field& x, - const Field& dD, - const Field& upper, - const Field& b -) const; - - -template<> -template<> -void BlockGaussSeidelPrecon::BlockSweep -( - Field& x, - const Field& dD, - const Field& upper, - const Field& lower, - const Field& b -) const; +void BlockGaussSeidelPrecon::calcInvDiag(); template<> From db2d2ff530a3506caef8c6c24552c0aa242cacfc Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 21 Aug 2015 08:16:36 +0100 Subject: [PATCH 7/7] Initial version: block diag Cholesky precon --- src/foam/Make/files | 5 + .../BlockDiagCholeskyPrecon.C | 1006 +++++++++++++++++ .../BlockDiagCholeskyPrecon.H | 235 ++++ .../BlockDiagCholeskyPreconDecoupled.C | 332 ++++++ .../blockDiagCholeskyPrecons.C | 42 + .../blockDiagCholeskyPrecons.H | 66 ++ .../scalarBlockDiagCholeskyPrecon.C | 199 ++++ .../scalarBlockDiagCholeskyPrecon.H | 75 ++ .../tensorBlockDiagCholeskyPrecon.C | 76 ++ .../tensorBlockDiagCholeskyPrecon.H | 75 ++ 10 files changed, 2111 insertions(+) create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPrecon.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPrecon.H create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPreconDecoupled.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/blockDiagCholeskyPrecons.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/blockDiagCholeskyPrecons.H create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/scalarBlockDiagCholeskyPrecon.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/scalarBlockDiagCholeskyPrecon.H create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/tensorBlockDiagCholeskyPrecon.C create mode 100644 src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/tensorBlockDiagCholeskyPrecon.H diff --git a/src/foam/Make/files b/src/foam/Make/files index f5b131619..ec2283a43 100644 --- a/src/foam/Make/files +++ b/src/foam/Make/files @@ -689,6 +689,11 @@ matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/blockGaussSeidelP matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/scalarBlockCholeskyPrecon.C matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/tensorBlockCholeskyPrecon.C matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/blockCholeskyPrecons.C + +matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/scalarBlockDiagCholeskyPrecon.C +matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/tensorBlockDiagCholeskyPrecon.C +matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/blockDiagCholeskyPrecons.C + matrices/blockLduMatrix/BlockLduPrecons/BlockAmgPrecon/blockAmgPrecons.C matrices/blockLduMatrix/BlockLduSmoothers/BlockLduSmoother/blockLduSmoothers.C diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPrecon.C new file mode 100644 index 000000000..f4e539f49 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPrecon.C @@ -0,0 +1,1006 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "error.H" +#include "BlockDiagCholeskyPrecon.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +void Foam::BlockDiagCholeskyPrecon::calcPreconDiag() +{ + // Note: Assuming lower and upper triangle have the same active type + + typedef CoeffField TypeCoeffField; + + if (this->matrix_.symmetric()) + { + const TypeCoeffField& UpperCoeff = this->matrix_.upper(); + + if (preconDiag_.activeType() == blockCoeffBase::SCALAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + diagMultiply + ( + preconDiag_.asScalar(), + UpperCoeff.asScalar() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + diagMultiply + ( + preconDiag_.asLinear(), + UpperCoeff.asLinear() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + // Transpose multiplication + diagMultiplyCoeffT + ( + preconDiag_.asSquare(), + UpperCoeff.asSquare() + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + diagMultiply + ( + preconDiag_.asLinear(), + UpperCoeff.asScalar() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + diagMultiply + ( + preconDiag_.asLinear(), + UpperCoeff.asLinear() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + // Transpose multiplication + diagMultiplyCoeffT + ( + preconDiag_.asSquare(), + UpperCoeff.asSquare() + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::SQUARE) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + diagMultiply + ( + preconDiag_.asSquare(), + UpperCoeff.asScalar() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + diagMultiply + ( + preconDiag_.asSquare(), + UpperCoeff.asLinear() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + // Transpose multiplication + diagMultiplyCoeffT + ( + preconDiag_.asSquare(), + UpperCoeff.asSquare() + ); + } + } + } + else // Asymmetric matrix + { + const TypeCoeffField& LowerCoeff = this->matrix_.lower(); + const TypeCoeffField& UpperCoeff = this->matrix_.upper(); + + if (preconDiag_.activeType() == blockCoeffBase::SCALAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + diagMultiply + ( + preconDiag_.asScalar(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + diagMultiply + ( + preconDiag_.asLinear(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + diagMultiply + ( + preconDiag_.asSquare(), + LowerCoeff.asSquare(), + UpperCoeff.asSquare() + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + diagMultiply + ( + preconDiag_.asLinear(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + diagMultiply + ( + preconDiag_.asLinear(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + diagMultiply + ( + preconDiag_.asSquare(), + LowerCoeff.asSquare(), + UpperCoeff.asSquare() + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::SQUARE) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + diagMultiply + ( + preconDiag_.asSquare(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + diagMultiply + ( + preconDiag_.asSquare(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + diagMultiply + ( + preconDiag_.asSquare(), + LowerCoeff.asSquare(), + UpperCoeff.asSquare() + ); + } + } + } + + // Invert the diagonal +// if (preconDiag_.activeType() == blockCoeffBase::SQUARE) +// { +// typedef typename TypeCoeffField::linearTypeField linearTypeField; +// typedef typename TypeCoeffField::linearType valueType; + +// typedef typename TypeCoeffField::squareTypeField squareTypeField; + +// // Special practice: invert only diagonal of diag coefficient +// Info<< "Special preconDiag" << endl; + +// // Get reference to active diag +// const squareTypeField& activeDiag = preconDiag_.asSquare(); + +// // Get reference to LU: remove diagonal from active diag +// squareTypeField& luDiag = LUDiag_.asSquare(); + +// linearTypeField lf(preconDiag_.size()); + +// // Take out the diagonal from the diag as a linear type +// contractLinear(lf, activeDiag); + +// expandLinear(luDiag, lf); + +// // Keep only off-diagonal in luDiag +// // Note change of sign to avoid multiplication with -1 when moving +// // to the other side. HJ, 20/Aug/2015 +// luDiag -= activeDiag; + +// // Store inverse of diagonal +// preconDiag_.clear(); + +// // Invert the diagonal part into lf +// preconDiag_.asLinear() = +// cmptDivide +// ( +// linearTypeField(lf.size(), pTraits::one), +// lf +// ); +// } +// else + { + preconDiag_ = inv(preconDiag_); + LUDiag_.asScalar() = 0; + } +} + + +template +template +void Foam::BlockDiagCholeskyPrecon::diagMultiply +( + Field& dDiag, + const Field& upper +) +{ + // Precondition the diagonal + + // Get addressing + const unallocLabelList& upperAddr = this->matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = this->matrix_.lduAddr().lowerAddr(); + + // Create multiplication function object + typename BlockCoeff::multiply mult; + + forAll (upper, coeffI) + { + dDiag[upperAddr[coeffI]] -= + mult.tripleProduct + ( + upper[coeffI], + dDiag[lowerAddr[coeffI]], + upper[coeffI] + ); + } +} + + +template +template +void Foam::BlockDiagCholeskyPrecon::diagMultiplyCoeffT +( + Field& dDiag, + const Field& upper +) +{ + // Precondition the diagonal + + // Get addressing + const unallocLabelList& upperAddr = this->matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = this->matrix_.lduAddr().lowerAddr(); + + // Create multiplication function object + typename BlockCoeff::multiply mult; + + forAll (upper, coeffI) + { + dDiag[upperAddr[coeffI]] -= + mult.tripleProduct + ( + upper[coeffI].T(), // Upper coefficient transposed + dDiag[lowerAddr[coeffI]], + upper[coeffI] + ); + } +} + + +template +template +void Foam::BlockDiagCholeskyPrecon::diagMultiply +( + Field& dDiag, + const Field& lower, + const Field& upper +) +{ + // Precondition the diagonal + + // Get addressing + const unallocLabelList& upperAddr = this->matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = this->matrix_.lduAddr().lowerAddr(); + + // Create multiplication function object + typename BlockCoeff::multiply mult; + + forAll (upper, coeffI) + { + dDiag[upperAddr[coeffI]] -= + mult.tripleProduct + ( + lower[coeffI], + dDiag[lowerAddr[coeffI]], + upper[coeffI] + ); + } +} + + +template +template +void Foam::BlockDiagCholeskyPrecon::ILUmultiply +( + Field& x, + const Field& dDiag, + const Field& upper, + const Field& b +) const +{ + // Create multiplication function object + typename BlockCoeff::multiply mult; + + forAll(x, i) + { + x[i] = mult(dDiag[i], b[i]); + } + + const unallocLabelList& upperAddr = this->matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = this->matrix_.lduAddr().lowerAddr(); + + forAll (upper, coeffI) + { + x[upperAddr[coeffI]] -= + mult + ( + dDiag[upperAddr[coeffI]], + mult(upper[coeffI], x[lowerAddr[coeffI]]) + ); + } + + forAllReverse (upper, coeffI) + { + x[lowerAddr[coeffI]] -= + mult + ( + dDiag[lowerAddr[coeffI]], + mult(upper[coeffI], x[upperAddr[coeffI]]) + ); + } +} + + +template +template +void Foam::BlockDiagCholeskyPrecon::ILUmultiplyCoeffT +( + Field& x, + const Field& dDiag, + const Field& upper, + const Field& b +) const +{ + // Create multiplication function object + typename BlockCoeff::multiply mult; + + forAll(x, i) + { + x[i] = mult(dDiag[i], b[i]); + } + + const unallocLabelList& upperAddr = this->matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = this->matrix_.lduAddr().lowerAddr(); + + forAll (upper, coeffI) + { + x[upperAddr[coeffI]] -= + mult + ( + dDiag[upperAddr[coeffI]], + mult(upper[coeffI].T(), x[lowerAddr[coeffI]]) + ); + } + + forAllReverse (upper, coeffI) + { + x[lowerAddr[coeffI]] -= + mult + ( + dDiag[lowerAddr[coeffI]], + mult(upper[coeffI], x[upperAddr[coeffI]]) + ); + } +} + + +template +template +void Foam::BlockDiagCholeskyPrecon::ILUmultiply +( + Field& x, + const Field& dDiag, + const Field& lower, + const Field& upper, + const Field& b +) const +{ + // Create multiplication function object + typename BlockCoeff::multiply mult; + + forAll(x, i) + { + x[i] = mult(dDiag[i], b[i]); + } + + const unallocLabelList& upperAddr = this->matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = this->matrix_.lduAddr().lowerAddr(); + const unallocLabelList& losortAddr = this->matrix_.lduAddr().losortAddr(); + + register label losortCoeff; + + forAll (lower, coeffI) + { + losortCoeff = losortAddr[coeffI]; + + x[upperAddr[losortCoeff]] -= + mult + ( + dDiag[upperAddr[losortCoeff]], + mult(lower[losortCoeff], x[lowerAddr[losortCoeff]]) + ); + } + + forAllReverse (upper, coeffI) + { + x[lowerAddr[coeffI]] -= + mult + ( + dDiag[lowerAddr[coeffI]], + mult(upper[coeffI], x[upperAddr[coeffI]]) + ); + } +} + + +template +template +void Foam::BlockDiagCholeskyPrecon::ILUmultiplyTranspose +( + Field& x, + const Field& dDiag, + const Field& lower, + const Field& upper, + const Field& b +) const +{ + // Create multiplication function object + typename BlockCoeff::multiply mult; + + forAll(x, i) + { + x[i] = mult(dDiag[i], b[i]); + } + + const unallocLabelList& upperAddr = this->matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = this->matrix_.lduAddr().lowerAddr(); + const unallocLabelList& losortAddr = this->matrix_.lduAddr().losortAddr(); + + register label losortCoeff; + + //HJ Not sure if the coefficient itself needs to be transposed. + // HJ, 30/Oct/2007 + forAll (lower, coeffI) + { + x[upperAddr[coeffI]] -= + mult + ( + dDiag[upperAddr[coeffI]], + mult(upper[coeffI], x[lowerAddr[coeffI]]) + ); + } + + forAllReverse (upper, coeffI) + { + losortCoeff = losortAddr[coeffI]; + + x[lowerAddr[losortCoeff]] -= + mult + ( + dDiag[lowerAddr[losortCoeff]], + mult(lower[losortCoeff], x[upperAddr[losortCoeff]]) + ); + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::BlockDiagCholeskyPrecon::BlockDiagCholeskyPrecon +( + const BlockLduMatrix& matrix +) +: + BlockLduPrecon(matrix), + preconDiag_(matrix.diag()), + LUDiag_(matrix.lduAddr().size()), + bPlusLU_() +{ + calcPreconDiag(); +} + + +template +Foam::BlockDiagCholeskyPrecon::BlockDiagCholeskyPrecon +( + const BlockLduMatrix& matrix, + const dictionary& dict +) +: + BlockLduPrecon(matrix), + preconDiag_(matrix.diag()), + LUDiag_(matrix.lduAddr().size()), + bPlusLU_() +{ + calcPreconDiag(); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::BlockDiagCholeskyPrecon::~BlockDiagCholeskyPrecon() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::BlockDiagCholeskyPrecon::precondition +( + Field& x, + const Field& b +) const +{ + typedef CoeffField TypeCoeffField; + + // Note: Assuming lower and upper triangle have the same active type + + if (this->matrix_.symmetric()) + { + const TypeCoeffField& DiagCoeff = this->matrix_.diag(); + const TypeCoeffField& UpperCoeff = this->matrix_.upper(); + + if (DiagCoeff.activeType() == blockCoeffBase::SCALAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + ILUmultiply + ( + x, + preconDiag_.asScalar(), + UpperCoeff.asScalar(), + b + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + ILUmultiply + ( + x, + preconDiag_.asScalar(), + UpperCoeff.asLinear(), + b + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + // Multiplication using transposed upper square coefficient + ILUmultiplyCoeffT + ( + x, + preconDiag_.asScalar(), + UpperCoeff.asSquare(), + b + ); + } + } + else if (DiagCoeff.activeType() == blockCoeffBase::LINEAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + ILUmultiply + ( + x, + preconDiag_.asLinear(), + UpperCoeff.asScalar(), + b + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + ILUmultiply + ( + x, + preconDiag_.asLinear(), + UpperCoeff.asLinear(), + b + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + // Multiplication using transposed upper square coefficient + ILUmultiplyCoeffT + ( + x, + preconDiag_.asLinear(), + UpperCoeff.asSquare(), + b + ); + } + } + else if (DiagCoeff.activeType() == blockCoeffBase::SQUARE) + { + // Add diag coupling to b + if (bPlusLU_.empty()) + { + bPlusLU_.setSize(b.size(), pTraits::zero); + } + + // Multiply overwrites bPlusLU_: no need to initialise + // Change of sign accounted via change of sign of bPlusLU_ + // HJ, 20/Aug/2015 + multiply(bPlusLU_, LUDiag_, x); + bPlusLU_ += b; + + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + // Note linear preconDiag due to decoupling + ILUmultiply + ( + x, + preconDiag_.asLinear(), + UpperCoeff.asScalar(), + bPlusLU_ + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + // Note linear preconDiag due to decoupling + ILUmultiply + ( + x, + preconDiag_.asLinear(), + UpperCoeff.asLinear(), + bPlusLU_ + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + // Multiplication using transposed upper square coefficient + // Note linear preconDiag due to decoupling + ILUmultiplyCoeffT + ( + x, + preconDiag_.asLinear(), + UpperCoeff.asSquare(), + bPlusLU_ + ); + } + } + } + else // Asymmetric matrix + { + const TypeCoeffField& DiagCoeff = this->matrix_.diag(); + const TypeCoeffField& LowerCoeff = this->matrix_.lower(); + const TypeCoeffField& UpperCoeff = this->matrix_.upper(); + + if (DiagCoeff.activeType() == blockCoeffBase::SCALAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + ILUmultiply + ( + x, + preconDiag_.asScalar(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar(), + b + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + ILUmultiply + ( + x, + preconDiag_.asScalar(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear(), + b + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + ILUmultiply + ( + x, + preconDiag_.asScalar(), + LowerCoeff.asSquare(), + UpperCoeff.asSquare(), + b + ); + } + } + else if (DiagCoeff.activeType() == blockCoeffBase::LINEAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + ILUmultiply + ( + x, + preconDiag_.asLinear(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar(), + b + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + ILUmultiply + ( + x, + preconDiag_.asLinear(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear(), + b + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + ILUmultiply + ( + x, + preconDiag_.asLinear(), + LowerCoeff.asSquare(), + UpperCoeff.asSquare(), + b + ); + } + } + else if (DiagCoeff.activeType() == blockCoeffBase::SQUARE) + { + // Add diag coupling to b + if (bPlusLU_.empty()) + { + bPlusLU_.setSize(b.size(), pTraits::zero); + } + + // Multiply overwrites bPlusLU_: no need to initialise + // Change of sign accounted via change of sign of bPlusLU_ + // HJ, 20/Aug/2015 + multiply(bPlusLU_, LUDiag_, x); + bPlusLU_ += b; + + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + // Note linear preconDiag due to decoupling + ILUmultiply + ( + x, + preconDiag_.asLinear(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar(), + bPlusLU_ + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + // Note linear preconDiag due to decoupling + ILUmultiply + ( + x, + preconDiag_.asLinear(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear(), + bPlusLU_ + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + // Note linear preconDiag due to decoupling + ILUmultiply + ( + x, + preconDiag_.asLinear(), + LowerCoeff.asSquare(), + UpperCoeff.asSquare(), + bPlusLU_ + ); + } + } + } +} + + +template +void Foam::BlockDiagCholeskyPrecon::preconditionT +( + Field& xT, + const Field& bT +) const +{ + typedef CoeffField TypeCoeffField; + + // Note: Assuming lower and upper triangle have the same active type + + if (this->matrix_.symmetric()) + { + precondition(xT, bT); + } + else // Asymmetric matrix + { + const TypeCoeffField& DiagCoeff = this->matrix_.diag(); + const TypeCoeffField& LowerCoeff = this->matrix_.lower(); + const TypeCoeffField& UpperCoeff = this->matrix_.upper(); + + if (DiagCoeff.activeType() == blockCoeffBase::SCALAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + ILUmultiplyTranspose + ( + xT, + preconDiag_.asScalar(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar(), + bT + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + ILUmultiplyTranspose + ( + xT, + preconDiag_.asScalar(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear(), + bT + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + ILUmultiplyTranspose + ( + xT, + preconDiag_.asScalar(), + LowerCoeff.asSquare(), + UpperCoeff.asSquare(), + bT + ); + } + } + else if (DiagCoeff.activeType() == blockCoeffBase::LINEAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + ILUmultiplyTranspose + ( + xT, + preconDiag_.asLinear(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar(), + bT + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + ILUmultiplyTranspose + ( + xT, + preconDiag_.asLinear(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear(), + bT + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + ILUmultiplyTranspose + ( + xT, + preconDiag_.asLinear(), + LowerCoeff.asSquare(), + UpperCoeff.asSquare(), + bT + ); + } + } + else if (DiagCoeff.activeType() == blockCoeffBase::SQUARE) + { + // Add diag coupling to b + if (bPlusLU_.empty()) + { + bPlusLU_.setSize(bT.size(), pTraits::zero); + } + + // Multiply overwrites bPlusLU_: no need to initialise + // Change of sign accounted via change of sign of bPlusLU_ + // HJ, 20/Aug/2015 + multiply(bPlusLU_, LUDiag_, xT); + bPlusLU_ += bT; + + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + // Note linear preconDiag due to decoupling + ILUmultiplyTranspose + ( + xT, + preconDiag_.asLinear(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar(), + bPlusLU_ + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + // Note linear preconDiag due to decoupling + ILUmultiplyTranspose + ( + xT, + preconDiag_.asLinear(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear(), + bPlusLU_ + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::SQUARE) + { + // Note linear preconDiag due to decoupling + ILUmultiplyTranspose + ( + xT, + preconDiag_.asLinear(), + LowerCoeff.asSquare(), + UpperCoeff.asSquare(), + bPlusLU_ + ); + } + } + } +} + + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPrecon.H new file mode 100644 index 000000000..0fb05e83c --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPrecon.H @@ -0,0 +1,235 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with foam-extend. If not, see . + +Class + BlockDiagCholeskyPrecon + +Description + Incomplete Cholesky preconditioning with no fill-in, + using the diagonal-of-diagonal for ILU decomposition. + Currently broken + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved. + +SourceFiles + BlockDiagCholeskyPrecon.C + BlockDiagCholeskyPreconDecoupled.C + +\*---------------------------------------------------------------------------*/ + +#ifndef BlockDiagCholeskyPrecon_H +#define BlockDiagCholeskyPrecon_H + +#include "BlockLduPrecon.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class BlockDiagCholeskyPrecon Declaration +\*---------------------------------------------------------------------------*/ + +template +class BlockDiagCholeskyPrecon +: + public BlockLduPrecon +{ + // Private Data + + //- Preconditioned diagonal + mutable CoeffField preconDiag_; + + //- Off-diag part of diagonal + CoeffField LUDiag_; + + //- Temporary space for updated decoupled source + // Initialised with zero size and resized on first use + mutable Field bPlusLU_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + BlockDiagCholeskyPrecon(const BlockDiagCholeskyPrecon&); + + //- Disallow default bitwise assignment + void operator=(const BlockDiagCholeskyPrecon&); + + //- Precondition the diagonal + void calcPreconDiag(); + + // Diagonal multiplication, symmetric matrix + template + void diagMultiply + ( + Field& dDiag, + const Field& upper + ); + + //- Diagonal multiplication with transpose upper square coeff + // for the symmetric matrix + template + void diagMultiplyCoeffT + ( + Field& dDiag, + const Field& upper + ); + + //- Diagonal multiplication, asymmetric matrix + template + void diagMultiply + ( + Field& dDiag, + const Field& lower, + const Field& upper + ); + + //- ILU multiplication, symmetric matrix + template + void ILUmultiply + ( + Field& x, + const Field& dDiag, + const Field& upper, + const Field& b + ) const; + + //- ILU multiplication, with transpose upper square coeff + // for a symmetric matrix + template + void ILUmultiplyCoeffT + ( + Field& x, + const Field& dDiag, + const Field& upper, + const Field& b + ) const; + + //- ILU multiplication, asymmetric matrix + template + void ILUmultiply + ( + Field& x, + const Field& dDiag, + const Field& lower, + const Field& upper, + const Field& b + ) const; + + //- ILU multiplication transposed asymmetric matrix + template + void ILUmultiplyTranspose + ( + Field& x, + const Field& dDiag, + const Field& lower, + const Field& upper, + const Field& b + ) const; + + + // Decoupled operations, used in template specialisation + + //- Precondition the diagonal, decoupled version + void calcDecoupledPreconDiag(); + + //- Execute preconditioning, decoupled version + void decoupledPrecondition + ( + Field& x, + const Field& b + ) const; + + //- Execute preconditioning with matrix transpose, + // decoupled version + void decoupledPreconditionT + ( + Field& xT, + const Field& bT + ) const; + + +public: + + //- Runtime type information + TypeName("Cholesky"); + + + // Constructors + + //- Construct from matrix for smoother use + BlockDiagCholeskyPrecon + ( + const BlockLduMatrix& matrix + ); + + //- Construct from components + BlockDiagCholeskyPrecon + ( + const BlockLduMatrix& matrix, + const dictionary& dict + ); + + + // Destructor + + virtual ~BlockDiagCholeskyPrecon(); + + + // Member Functions + + //- Execute preconditioning + virtual void precondition + ( + Field& x, + const Field& b + ) const; + + //- Execute preconditioning with matrix transpose + virtual void preconditionT + ( + Field& xT, + const Field& bT + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "BlockDiagCholeskyPrecon.C" +# include "BlockDiagCholeskyPreconDecoupled.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPreconDecoupled.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPreconDecoupled.C new file mode 100644 index 000000000..f123d689e --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/BlockDiagCholeskyPreconDecoupled.C @@ -0,0 +1,332 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "error.H" +#include "BlockDiagCholeskyPrecon.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +void Foam::BlockDiagCholeskyPrecon::calcDecoupledPreconDiag() +{ + typedef CoeffField TypeCoeffField; + + // Note: Assuming lower and upper triangle have the same active type + + if (this->matrix_.symmetric()) + { + const TypeCoeffField& UpperCoeff = this->matrix_.upper(); + + if (preconDiag_.activeType() == blockCoeffBase::SCALAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + diagMultiply + ( + preconDiag_.asScalar(), + UpperCoeff.asScalar() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + diagMultiply + ( + preconDiag_.asLinear(), + UpperCoeff.asLinear() + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + diagMultiply + ( + preconDiag_.asLinear(), + UpperCoeff.asScalar() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + diagMultiply + ( + preconDiag_.asLinear(), + UpperCoeff.asLinear() + ); + } + } + } + else // Asymmetric matrix + { + const TypeCoeffField& LowerCoeff = this->matrix_.lower(); + const TypeCoeffField& UpperCoeff = this->matrix_.upper(); + + if (preconDiag_.activeType() == blockCoeffBase::SCALAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + diagMultiply + ( + preconDiag_.asScalar(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + diagMultiply + ( + preconDiag_.asLinear(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear() + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + diagMultiply + ( + preconDiag_.asLinear(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar() + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + diagMultiply + ( + preconDiag_.asLinear(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear() + ); + } + } + } + + // Invert the diagonal + // Note: since square diag type does not exist, simple inversion + // covers all cases + preconDiag_ = inv(preconDiag_); +} + + +template +void Foam::BlockDiagCholeskyPrecon::decoupledPrecondition +( + Field& x, + const Field& b +) const +{ + typedef CoeffField TypeCoeffField; + + // Note: Assuming lower and upper triangle have the same active type + + if (this->matrix_.symmetric()) + { + const TypeCoeffField& UpperCoeff = this->matrix_.upper(); + + if (preconDiag_.activeType() == blockCoeffBase::SCALAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + ILUmultiply + ( + x, + preconDiag_.asScalar(), + UpperCoeff.asScalar(), + b + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + ILUmultiply + ( + x, + preconDiag_.asScalar(), + UpperCoeff.asLinear(), + b + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + ILUmultiply + ( + x, + preconDiag_.asLinear(), + UpperCoeff.asScalar(), + b + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + ILUmultiply + ( + x, + preconDiag_.asLinear(), + UpperCoeff.asLinear(), + b + ); + } + } + } + else // Asymmetric matrix + { + const TypeCoeffField& LowerCoeff = this->matrix_.lower(); + const TypeCoeffField& UpperCoeff = this->matrix_.upper(); + + if (preconDiag_.activeType() == blockCoeffBase::SCALAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + ILUmultiply + ( + x, + preconDiag_.asScalar(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar(), + b + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + ILUmultiply + ( + x, + preconDiag_.asScalar(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear(), + b + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + ILUmultiply + ( + x, + preconDiag_.asLinear(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar(), + b + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + ILUmultiply + ( + x, + preconDiag_.asLinear(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear(), + b + ); + } + } + } +} + + +template +void Foam::BlockDiagCholeskyPrecon::decoupledPreconditionT +( + Field& xT, + const Field& bT +) const +{ + typedef CoeffField TypeCoeffField; + + // Note: Assuming lower and upper triangle have the same active type + + if (this->matrix_.symmetric()) + { + precondition(xT, bT); + } + else // Asymmetric matrix + { + const TypeCoeffField& LowerCoeff = this->matrix_.lower(); + const TypeCoeffField& UpperCoeff = this->matrix_.upper(); + + if (preconDiag_.activeType() == blockCoeffBase::SCALAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + ILUmultiplyTranspose + ( + xT, + preconDiag_.asScalar(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar(), + bT + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + ILUmultiplyTranspose + ( + xT, + preconDiag_.asScalar(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear(), + bT + ); + } + } + else if (preconDiag_.activeType() == blockCoeffBase::LINEAR) + { + if (UpperCoeff.activeType() == blockCoeffBase::SCALAR) + { + ILUmultiplyTranspose + ( + xT, + preconDiag_.asLinear(), + LowerCoeff.asScalar(), + UpperCoeff.asScalar(), + bT + ); + } + else if (UpperCoeff.activeType() == blockCoeffBase::LINEAR) + { + ILUmultiplyTranspose + ( + xT, + preconDiag_.asLinear(), + LowerCoeff.asLinear(), + UpperCoeff.asLinear(), + bT + ); + } + } + } +} + + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/blockDiagCholeskyPrecons.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/blockDiagCholeskyPrecons.C new file mode 100644 index 000000000..4082ef2d6 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/blockDiagCholeskyPrecons.C @@ -0,0 +1,42 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "blockLduMatrices.H" +#include "blockLduPrecons.H" +#include "blockDiagCholeskyPrecons.H" +#include "addToRunTimeSelectionTable.H" + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +makeBlockPrecons(blockDiagCholeskyPrecon); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/blockDiagCholeskyPrecons.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/blockDiagCholeskyPrecons.H new file mode 100644 index 000000000..292ab3771 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/blockDiagCholeskyPrecons.H @@ -0,0 +1,66 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with foam-extend. If not, see . + +Class + BlockDiagCholeskyPrecon + +Description + Typedefs for Incomplete Cholesky preconditioning with no fill-in, + using the diagonal-of-diagonal for ILU decomposition. + Currently broken + + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved. + +SourceFiles + blockDiagCholeskyPrecons.C + +\*---------------------------------------------------------------------------*/ + +#ifndef blockDiagCholeskyPrecons_H +#define blockDiagCholeskyPrecons_H + +#include "scalarBlockDiagCholeskyPrecon.H" +#include "tensorBlockDiagCholeskyPrecon.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +typedef BlockDiagCholeskyPrecon blockDiagCholeskyPreconScalar; +typedef BlockDiagCholeskyPrecon blockDiagCholeskyPreconVector; +typedef BlockDiagCholeskyPrecon blockDiagCholeskyPreconTensor; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/scalarBlockDiagCholeskyPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/scalarBlockDiagCholeskyPrecon.C new file mode 100644 index 000000000..6127ae309 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/scalarBlockDiagCholeskyPrecon.C @@ -0,0 +1,199 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#ifndef scalarBlockDiagCholeskyPrecon_H +#define scalarBlockDiagCholeskyPrecon_H + +#include "BlockDiagCholeskyPrecon.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template<> +void Foam::BlockDiagCholeskyPrecon::calcPreconDiag() +{ + // Precondition the diagonal + if (matrix_.symmetric()) + { + const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); + + // Get off-diagonal matrix coefficients + const scalarField& upper = matrix_.upper(); + + forAll (upper, coeffI) + { + preconDiag_[upperAddr[coeffI]] -= + sqr(upper[coeffI])/preconDiag_[lowerAddr[coeffI]]; + } + } + else if (matrix_.asymmetric()) + { + const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); + + // Get off-diagonal matrix coefficients + const scalarField& upper = matrix_.upper(); + const scalarField& lower = matrix_.lower(); + + forAll (upper, coeffI) + { + preconDiag_[upperAddr[coeffI]] -= + upper[coeffI]*lower[coeffI]/preconDiag_[lowerAddr[coeffI]]; + } + } + + // Invert the diagonal for future use + forAll (preconDiag_, i) + { + preconDiag_[i] = 1.0/preconDiag_[i]; + } +} + + +template<> +void Foam::BlockDiagCholeskyPrecon::precondition +( + scalarField& x, + const scalarField& b +) const +{ + forAll(x, i) + { + x[i] = b[i]*preconDiag_[i]; + } + + if (matrix_.symmetric()) + { + const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); + + // Get off-diagonal matrix coefficients + const scalarField& upper = matrix_.upper(); + + forAll (upper, coeffI) + { + x[upperAddr[coeffI]] -= + preconDiag_[upperAddr[coeffI]]* + upper[coeffI]*x[lowerAddr[coeffI]]; + } + + forAllReverse (upper, coeffI) + { + x[lowerAddr[coeffI]] -= + preconDiag_[lowerAddr[coeffI]]* + upper[coeffI]*x[upperAddr[coeffI]]; + } + } + else if (matrix_.asymmetric()) + { + const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); + const unallocLabelList& losortAddr = matrix_.lduAddr().losortAddr(); + + // Get off-diagonal matrix coefficients + const scalarField& upper = matrix_.upper(); + const scalarField& lower = matrix_.lower(); + + label losortCoeff; + + forAll (lower, coeffI) + { + losortCoeff = losortAddr[coeffI]; + + x[upperAddr[losortCoeff]] -= + preconDiag_[upperAddr[losortCoeff]]* + lower[losortCoeff]*x[lowerAddr[losortCoeff]]; + } + + forAllReverse (upper, coeffI) + { + x[lowerAddr[coeffI]] -= + preconDiag_[lowerAddr[coeffI]]* + upper[coeffI]*x[upperAddr[coeffI]]; + } + } +} + + +template<> +void Foam::BlockDiagCholeskyPrecon::preconditionT +( + scalarField& xT, + const scalarField& bT +) const +{ + if (matrix_.symmetric()) + { + precondition(xT, bT); + } + + forAll(xT, i) + { + xT[i] = bT[i]*preconDiag_[i]; + } + + if (matrix_.asymmetric()) + { + const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr(); + const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr(); + const unallocLabelList& losortAddr = matrix_.lduAddr().losortAddr(); + + // Get off-diagonal matrix coefficients + const scalarField& upper = matrix_.upper(); + const scalarField& lower = matrix_.lower(); + + label losortCoeff; + + forAll (lower, coeffI) + { + xT[upperAddr[coeffI]] -= + preconDiag_[upperAddr[coeffI]]* + upper[coeffI]*xT[lowerAddr[coeffI]]; + } + + forAllReverse (upper, coeffI) + { + losortCoeff = losortAddr[coeffI]; + + xT[lowerAddr[losortCoeff]] -= + preconDiag_[lowerAddr[losortCoeff]]* + lower[losortCoeff]*xT[upperAddr[losortCoeff]]; + } + } +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/scalarBlockDiagCholeskyPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/scalarBlockDiagCholeskyPrecon.H new file mode 100644 index 000000000..8dabdaee4 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/scalarBlockDiagCholeskyPrecon.H @@ -0,0 +1,75 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with foam-extend. If not, see . + +Class + BlockDiagCholeskyPrecon + +Description + Template specialisation for scalar block Cholesky preconditioning + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved. + +SourceFiles + scalarBlockDiagCholeskyPrecon.C + +\*---------------------------------------------------------------------------*/ + +#ifndef scalarBlockDiagCholeskyPrecon_H +#define scalarBlockDiagCholeskyPrecon_H + +#include "BlockDiagCholeskyPrecon.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template<> +void Foam::BlockDiagCholeskyPrecon::calcPreconDiag(); + +template<> +void Foam::BlockDiagCholeskyPrecon::precondition +( + scalarField& x, + const scalarField& b +) const; + + +template<> +void Foam::BlockDiagCholeskyPrecon::preconditionT +( + scalarField& xT, + const scalarField& bT +) const; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/tensorBlockDiagCholeskyPrecon.C b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/tensorBlockDiagCholeskyPrecon.C new file mode 100644 index 000000000..4991d81e2 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/tensorBlockDiagCholeskyPrecon.C @@ -0,0 +1,76 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#ifndef tensorBlockDiagCholeskyPrecon_H +#define tensorBlockDiagCholeskyPrecon_H + +#include "BlockDiagCholeskyPrecon.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template<> +void Foam::BlockDiagCholeskyPrecon::calcPreconDiag() +{ + // Decoupled version + calcDecoupledPreconDiag(); +} + + +template<> +void Foam::BlockDiagCholeskyPrecon::precondition +( + tensorField& x, + const tensorField& b +) const +{ + // Decoupled version + decoupledPrecondition(x, b); +} + + +template<> +void Foam::BlockDiagCholeskyPrecon::preconditionT +( + tensorField& xT, + const tensorField& bT +) const +{ + // Decoupled version + decoupledPreconditionT(xT, bT); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/tensorBlockDiagCholeskyPrecon.H b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/tensorBlockDiagCholeskyPrecon.H new file mode 100644 index 000000000..c55f18096 --- /dev/null +++ b/src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagCholeskyPrecon/tensorBlockDiagCholeskyPrecon.H @@ -0,0 +1,75 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation, either version 3 of the License, or (at your + option) any later version. + + foam-extend is distributed in the hope that it will be useful, but + WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + + You should have received a copy of the GNU General Public License + along with foam-extend. If not, see . + +Class + BlockDiagCholeskyPrecon + +Description + Template specialisation for tensor block Cholesky preconditioning + +Author + Hrvoje Jasak, Wikki Ltd. All rights reserved. + +SourceFiles + tensorBlockDiagCholeskyPrecon.C + +\*---------------------------------------------------------------------------*/ + +#ifndef tensorBlockDiagCholeskyPrecon_H +#define tensorBlockDiagCholeskyPrecon_H + +#include "BlockDiagCholeskyPrecon.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +template<> +void Foam::BlockDiagCholeskyPrecon::calcPreconDiag(); + +template<> +void Foam::BlockDiagCholeskyPrecon::precondition +( + tensorField& x, + const tensorField& b +) const; + + +template<> +void Foam::BlockDiagCholeskyPrecon::preconditionT +( + tensorField& xT, + const tensorField& bT +) const; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //