From 7794e764174f0f91553cbee36da07c308e78fced Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Fri, 7 Jul 2017 13:29:27 +0100 Subject: [PATCH] Added master() function in lduInterface: better parallel ILU preconditioning --- .../faPatches/basic/coupled/coupledFaPatch.H | 5 +- .../constraint/cyclic/cyclicFaPatch.H | 13 +- .../constraint/processor/processorFaPatch.H | 8 +- .../fvPatches/basic/coupled/coupledFvPatch.H | 11 +- .../fvMesh/fvPatches/fvPatch/fvPatch.H | 4 +- .../lduMatrix/lduAddressing/lduAddressing.C | 180 +++++++++++++++++- .../lduMatrix/lduAddressing/lduAddressing.H | 52 ++++- .../lduInterfaces/lduInterface/lduInterface.H | 3 + .../cyclicAMGInterface/cyclicAMGInterface.H | 9 + .../processorAMGInterface.H | 6 + .../processorSAMGInterface.H | 6 + .../basic/coupled/coupledPointPatch.H | 5 +- .../constraint/cyclic/cyclicPointPatch.H | 5 +- .../derived/coupled/coupledFacePointPatch.H | 14 +- .../derived/global/globalPointPatch.H | 14 +- .../constraint/global/globalTetPolyPatch.H | 14 +- .../processor/processorTetPolyPatch.H | 11 +- 17 files changed, 324 insertions(+), 36 deletions(-) diff --git a/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H b/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H index 693dba339..f46e51672 100644 --- a/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H +++ b/src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H @@ -132,9 +132,8 @@ public: {} - // Destructor - - virtual ~coupledFaPatch(); + //- Destructor + virtual ~coupledFaPatch(); // Member Functions diff --git a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H index 2ba7e8cbc..d81481e94 100644 --- a/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H +++ b/src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H @@ -98,16 +98,21 @@ public: {} - // Destructor - - virtual ~cyclicFaPatch() - {} + //- Destructor + virtual ~cyclicFaPatch() + {} // Member functions // Access + //- Is this the master side? Yes: it contains both sets of faces + virtual bool master() const + { + return true; + } + //- Return face transformation tensor virtual const tensorField& forwardT() const { diff --git a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H index cb3314fb5..4fe7a59cb 100644 --- a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H +++ b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H @@ -176,7 +176,7 @@ public: return size(); } - //- Return processor number + //- Return processor number int myProcNo() const { return myProcNo_; @@ -201,6 +201,12 @@ public: } } + //- Is this the master side? + virtual bool master() const + { + return (myProcNo_ < neighbProcNo_); + } + // Communications support diff --git a/src/finiteVolume/fvMesh/fvPatches/basic/coupled/coupledFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/basic/coupled/coupledFvPatch.H index 1592aaa1f..bf96c2f24 100644 --- a/src/finiteVolume/fvMesh/fvPatches/basic/coupled/coupledFvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/basic/coupled/coupledFvPatch.H @@ -92,9 +92,8 @@ public: {} - // Destructor - - virtual ~coupledFvPatch(); + //- Destructor + virtual ~coupledFvPatch(); // Member Functions @@ -107,6 +106,12 @@ public: return coupledPolyPatch_.coupled(); } + //- Return true if patch is coupled + virtual bool master() const + { + return coupledPolyPatch_.master(); + } + //- Return face transformation tensor const tensorField& forwardT() const { diff --git a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H index fd0df4d41..12a7722af 100644 --- a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H @@ -136,14 +136,14 @@ public: static autoPtr New(const polyPatch&, const fvBoundaryMesh&); - // Destructor + //- Destructor virtual ~fvPatch(); // Member Functions - // Access + // Access //- Return the polyPatch const polyPatch& patch() const diff --git a/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.C b/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.C index 32098657d..ca87c406c 100644 --- a/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.C +++ b/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.C @@ -26,6 +26,7 @@ License #include "lduAddressing.H" #include "extendedLduAddressing.H" #include "demandDrivenData.H" +#include "dynamicLabelList.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -170,6 +171,120 @@ void Foam::lduAddressing::calcLosortStart() const } +void Foam::lduAddressing::calcInternalBoundaryEqnCoeffs +( + const lduInterfaceFieldPtrsList& lduInterfaces +) const +{ + if + ( + internalEqnCoeffsPtr_ + || flippedInternalEqnCoeffsPtr_ + || boundaryEqnCoeffs_.size() + ) + { + FatalErrorIn("lduAddressing::calcInternalBoundaryCoeffs() const") + << "Internal/boundary equation coefficients already calculated." + << abort(FatalError); + } + + // Get number of internal coefficients (number of faces) + const label nInternalCoeffs = upperAddr().size(); + + // Allocate insertion-friendly storage with enough memory + dynamicLabelList internalCoeffs(nInternalCoeffs); + dynamicLabelList flippedInternalCoeffs(nInternalCoeffs); + + // Initialise boundary equation coefficients for all coupled patches with + // enough storage + boundaryEqnCoeffs_.setSize(lduInterfaces.size()); + forAll (lduInterfaces, intI) + { + if (lduInterfaces.set(intI)) + { + boundaryEqnCoeffs_.set + ( + intI, + new dynamicLabelList + ( + lduInterfaces[intI].coupledInterface().faceCells().size() + ) + ); + } + } + + // First, we need to mark boundary equations with associated interface + // index. Note: this does not take into account corner cells that have + // faces on more than one coupled interface. Don't care about them during + // preconditioning at the moment. VV, 23/Jun/2017. + labelList boundaryEqnInterfaceIndex(size_, -1); + + // Loop through interfaces + forAll (lduInterfaces, intI) + { + // Check whether the interface is set + if (lduInterfaces.set(intI)) + { + // Get boundary equations/rows (face cells) + const unallocLabelList& boundaryEqns = + lduInterfaces[intI].coupledInterface().faceCells(); + + // Loop through boundary equations and mark them + forAll (boundaryEqns, beI) + { + boundaryEqnInterfaceIndex[boundaryEqns[beI]] = intI; + } + } + } + + // Get lower/upper (owner/neighbour) addressing + const unallocLabelList& own = lowerAddr(); + const unallocLabelList& nei = upperAddr(); + + // Loop through upper triangle and filter coefficients (faces) + forAll (own, coeffI) + { + // Get owner/neighbour (row/column) for this coefficient + const label& ownI = own[coeffI]; + const label& neiI = nei[coeffI]; + + if + ( + boundaryEqnInterfaceIndex[ownI] != -1 + && boundaryEqnInterfaceIndex[neiI] != -1 + ) + { + // Both owner and neigbour are at the boundary, append to boundary + // coeffs list. Note: it is possible that owner/neighbour do not + // have the same interface index since we ignored corner cells. Take + // the owner index + boundaryEqnCoeffs_[boundaryEqnInterfaceIndex[ownI]].append(coeffI); + } + else + { + // If owner is a boundary cell and neighbour is not a boundary cell, + // we need to mark this face for flipping + if + ( + boundaryEqnInterfaceIndex[ownI] != -1 + && boundaryEqnInterfaceIndex[neiI] == -1 + ) + { + flippedInternalCoeffs.append(coeffI); + } + else + { + internalCoeffs.append(coeffI); + } + } + } + + // Reuse dynamic lists to initialise data members + internalEqnCoeffsPtr_ = new labelList(internalCoeffs.xfer()); + flippedInternalEqnCoeffsPtr_ = new labelList(flippedInternalCoeffs.xfer()); +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::lduAddressing::lduAddressing(const label nEqns) @@ -178,7 +293,10 @@ Foam::lduAddressing::lduAddressing(const label nEqns) losortPtr_(NULL), ownerStartPtr_(NULL), losortStartPtr_(NULL), - extendedAddr_(5) + extendedAddr_(5), + internalEqnCoeffsPtr_(NULL), + flippedInternalEqnCoeffsPtr_(NULL), + boundaryEqnCoeffs_() {} @@ -189,6 +307,8 @@ Foam::lduAddressing::~lduAddressing() deleteDemandDrivenData(losortPtr_); deleteDemandDrivenData(ownerStartPtr_); deleteDemandDrivenData(losortStartPtr_); + deleteDemandDrivenData(internalEqnCoeffsPtr_); + deleteDemandDrivenData(flippedInternalEqnCoeffsPtr_); } @@ -284,4 +404,62 @@ Foam::lduAddressing::extendedAddr(const label p) const } +const Foam::unallocLabelList& Foam::lduAddressing::internalEqnCoeffs +( + const lduInterfaceFieldPtrsList& lduInterfaces +) const +{ + if (!internalEqnCoeffsPtr_) + { + calcInternalBoundaryEqnCoeffs(lduInterfaces); + } + + return *internalEqnCoeffsPtr_; +} + + +const Foam::unallocLabelList& Foam::lduAddressing::flippedInternalEqnCoeffs +( + const lduInterfaceFieldPtrsList& lduInterfaces +) const +{ + if (!flippedInternalEqnCoeffsPtr_) + { + calcInternalBoundaryEqnCoeffs(lduInterfaces); + } + + return *flippedInternalEqnCoeffsPtr_; +} + + +const Foam::dynamicLabelList& Foam::lduAddressing::boundaryEqnCoeffs +( + const lduInterfaceFieldPtrsList& lduInterfaces, + const label intI +) const +{ + if (intI > lduInterfaces.size() - 1 || intI < 0) + { + FatalErrorIn + ( + "const Foam::PtrList& " + "Foam::lduAddressing::boundaryEqnCoeffs" + "\n(" + "\n const lduInterfaceFieldPtrsList& lduInterfaces," + "\n const label p" + "\n) const" + ) << "Invalid interface index specified: " << intI << nl + << "Number of coupled interfaces: " << lduInterfaces.size() + << abort(FatalError); + } + + if (!boundaryEqnCoeffs_.set(intI)) + { + calcInternalBoundaryEqnCoeffs(lduInterfaces); + } + + return boundaryEqnCoeffs_[intI]; +} + + // ************************************************************************* // diff --git a/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.H b/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.H index 6597ac879..c54fb7e56 100644 --- a/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.H +++ b/src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.H @@ -98,6 +98,8 @@ SourceFiles #include "labelList.H" #include "lduSchedule.H" +#include "lduInterfaceFieldPtrsList.H" +#include "boolList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -130,8 +132,27 @@ class lduAddressing //- Losort start addressing mutable labelList* losortStartPtr_; - //- Extended addressing with p-order fill-in - mutable PtrList extendedAddr_; + + // Demand-driven data for ILU precondition with p-order fill in (ILUCp) + + //- Extended addressing with p-order fill-in + mutable PtrList extendedAddr_; + + // Demand-driven data for parallelisation of ILU0 and Cholesky + // preconditioners where the order of factorisation is important + + //- List of coefficients for internal equations (e.g. faces where + // at least one cell is internal cell) + mutable labelList* internalEqnCoeffsPtr_; + + //- List of coefficients where the owner cell is at the boundary and + // the neighbour cell is internal. Needed for correct ordering + mutable labelList* flippedInternalEqnCoeffsPtr_; + + //- PtrList of coefficients for boundary equations for each coupled + // interface (e.g. faces where both cells are boundary cells on a + // given processor patch) + mutable PtrList boundaryEqnCoeffs_; // Private Member Functions @@ -151,6 +172,13 @@ class lduAddressing //- Calculate losort start void calcLosortStart() const; + //- Calculate internal/boundary equation coefficients given a list of + // interfaces + void calcInternalBoundaryEqnCoeffs + ( + const lduInterfaceFieldPtrsList& lduInterfaces + ) const; + public: @@ -202,6 +230,26 @@ public: //- Return extended addressing given p index const extendedLduAddressing& extendedAddr(const label p) const; + + //- Return list of coefficients for internal equations + const unallocLabelList& internalEqnCoeffs + ( + const lduInterfaceFieldPtrsList& lduInterfaces + ) const; + + //- Return internal cell/boundary cell face flips + const unallocLabelList& flippedInternalEqnCoeffs + ( + const lduInterfaceFieldPtrsList& lduInterfaces + ) const; + + //- Return list of coefficients for boundary equations given list of + // ldu interface fields and coupled interface index + const dynamicLabelList& boundaryEqnCoeffs + ( + const lduInterfaceFieldPtrsList& lduInterfaces, + const label intI + ) const; }; diff --git a/src/foam/matrices/lduMatrix/lduAddressing/lduInterfaces/lduInterface/lduInterface.H b/src/foam/matrices/lduMatrix/lduAddressing/lduInterfaces/lduInterface/lduInterface.H index edcb900fb..cecffef7b 100644 --- a/src/foam/matrices/lduMatrix/lduAddressing/lduInterfaces/lduInterface/lduInterface.H +++ b/src/foam/matrices/lduMatrix/lduAddressing/lduInterfaces/lduInterface/lduInterface.H @@ -91,6 +91,9 @@ public: //- Return true if interface is coupled virtual bool coupled() const = 0; + //- Does this side own the patch ? + virtual bool master() const = 0; + //- Return faceCell addressing virtual const unallocLabelList& faceCells() const = 0; diff --git a/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/AMGInterfaces/cyclicAMGInterface/cyclicAMGInterface.H b/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/AMGInterfaces/cyclicAMGInterface/cyclicAMGInterface.H index 90d260e35..1ac72dec9 100644 --- a/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/AMGInterfaces/cyclicAMGInterface/cyclicAMGInterface.H +++ b/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/AMGInterfaces/cyclicAMGInterface/cyclicAMGInterface.H @@ -102,6 +102,15 @@ public: return true; } + //- Does this side own the patch ? + // HACKED: the patch contains both master and slave + // This influences parallel Cholesky and ILU preconditioning + // Please use cyclicGgi instead. HJ, 22/Jun/2017 + virtual bool master() const + { + return true; + } + // Interface transfer functions diff --git a/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/AMGInterfaces/processorAMGInterface/processorAMGInterface.H b/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/AMGInterfaces/processorAMGInterface/processorAMGInterface.H index 716b76d3a..2c8a0b18c 100644 --- a/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/AMGInterfaces/processorAMGInterface/processorAMGInterface.H +++ b/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/AMGInterfaces/processorAMGInterface/processorAMGInterface.H @@ -114,6 +114,12 @@ public: return true; } + //- Does this side own the patch ? + virtual bool master() const + { + return myProcNo() < neighbProcNo(); + } + // Communications support diff --git a/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/processorSAMGInterface/processorSAMGInterface.H b/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/processorSAMGInterface/processorSAMGInterface.H index 1d924aa27..1fb93244b 100644 --- a/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/processorSAMGInterface/processorSAMGInterface.H +++ b/src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/processorSAMGInterface/processorSAMGInterface.H @@ -114,6 +114,12 @@ public: return true; } + //- Does this side own the patch ? + virtual bool master() const + { + return myProcNo() < neighbProcNo(); + } + // Communications support diff --git a/src/foam/meshes/pointMesh/pointPatches/basic/coupled/coupledPointPatch.H b/src/foam/meshes/pointMesh/pointPatches/basic/coupled/coupledPointPatch.H index f5ec5bf14..77c17a0a1 100644 --- a/src/foam/meshes/pointMesh/pointPatches/basic/coupled/coupledPointPatch.H +++ b/src/foam/meshes/pointMesh/pointPatches/basic/coupled/coupledPointPatch.H @@ -97,9 +97,8 @@ public: coupledPointPatch(const pointBoundaryMesh& bm); - // Destructor - - virtual ~coupledPointPatch(); + //- Destructor + virtual ~coupledPointPatch(); // Member Functions diff --git a/src/foam/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H b/src/foam/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H index 678cc5b7d..166e037a2 100644 --- a/src/foam/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H +++ b/src/foam/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H @@ -116,9 +116,8 @@ public: ); - // Destructor - - virtual ~cyclicPointPatch(); + //- Destructor + virtual ~cyclicPointPatch(); // Member Functions diff --git a/src/foam/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H b/src/foam/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H index 390ace2cc..1c74c2fa2 100644 --- a/src/foam/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H +++ b/src/foam/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H @@ -86,9 +86,8 @@ public: ); - // Destructor - - virtual ~coupledFacePointPatch(); + //- Destructor + virtual ~coupledFacePointPatch(); // Member Functions @@ -100,6 +99,15 @@ public: { return true; } + + //- Does this side own the patch ? + // HACKED: the patch contains both master and slave + // This influences parallel Cholesky and ILU preconditioning + // Please use cyclicGgi instead. HJ, 22/Jun/2017 + virtual bool master() const + { + return true; + } }; diff --git a/src/foam/meshes/pointMesh/pointPatches/derived/global/globalPointPatch.H b/src/foam/meshes/pointMesh/pointPatches/derived/global/globalPointPatch.H index 25afaf5ee..adec71633 100644 --- a/src/foam/meshes/pointMesh/pointPatches/derived/global/globalPointPatch.H +++ b/src/foam/meshes/pointMesh/pointPatches/derived/global/globalPointPatch.H @@ -119,9 +119,8 @@ public: ); - // Destructor - - virtual ~globalPointPatch(); + //- Destructor + virtual ~globalPointPatch(); // Member functions @@ -153,6 +152,15 @@ public: } } + //- Does this side own the patch ? + // HACKED: the patch contains both master and slave + // This influences parallel Cholesky and ILU preconditioning + // HJ, 22/Jun/2017 + virtual bool master() const + { + return true; + } + //- Return number of faces virtual label nFaces() const { diff --git a/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/constraint/global/globalTetPolyPatch.H b/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/constraint/global/globalTetPolyPatch.H index 2a7cfe0ae..6d3f1eb16 100644 --- a/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/constraint/global/globalTetPolyPatch.H +++ b/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/constraint/global/globalTetPolyPatch.H @@ -197,9 +197,9 @@ public: ); - // Destructor + //- Destructor + virtual ~globalTetPolyPatch(); - virtual ~globalTetPolyPatch(); // Member functions @@ -223,6 +223,15 @@ public: return 0; } + //- Does this side own the patch ? + // HACKED: the patch contains both master and slave + // This influences parallel Cholesky and ILU preconditioning + // HJ, 22/Jun/2017 + virtual bool master() const + { + return true; + } + //- Return total number of shared points virtual label globalPointSize() const { @@ -292,6 +301,7 @@ public: return meshCutEdgeMask_; } + // Access functions for demand driven data //- Return list of edge indices for edges local to the patch diff --git a/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/constraint/processor/processorTetPolyPatch.H b/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/constraint/processor/processorTetPolyPatch.H index 3a787778d..c03a09049 100644 --- a/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/constraint/processor/processorTetPolyPatch.H +++ b/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/constraint/processor/processorTetPolyPatch.H @@ -172,9 +172,8 @@ public: ); - // Destructor - - virtual ~processorTetPolyPatch(); + //- Destructor + virtual ~processorTetPolyPatch(); // Member functions @@ -192,15 +191,15 @@ public: } //- Is this a master patch - bool isMaster() const + bool master() const { return myProcNo() < neighbProcNo(); } //- Is this a slave patch - bool isSlave() const + bool slave() const { - return !isMaster(); + return !master(); } //- Return the underlying processorPolyPatch