Added master() function in lduInterface: better parallel ILU preconditioning

src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H

@@ -132,9 +132,8 @@ public:
         {}
 
 
-    // Destructor
-
-        virtual ~coupledFaPatch();
+    //- Destructor
+    virtual ~coupledFaPatch();
 
 
     // Member Functions

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
             {

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
 

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
             {

src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H

@@ -136,14 +136,14 @@ public:
         static autoPtr<fvPatch> New(const polyPatch&, const fvBoundaryMesh&);
 
 
-    // Destructor
+    //- Destructor
 
         virtual ~fvPatch();
 
 
     // Member Functions
 
-    // Access
+        // Access
 
             //- Return the polyPatch
             const polyPatch& patch() const

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<labelList>& "
+            "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];
+}
+
+
 // ************************************************************************* //

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<extendedLduAddressing> extendedAddr_;
+
+        // Demand-driven data for ILU precondition with p-order fill in (ILUCp)
+
+            //- Extended addressing with p-order fill-in
+            mutable PtrList<extendedLduAddressing> 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<dynamicLabelList> 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;
 };
 
 

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;
 

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
 

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
 

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
 

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

src/foam/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H

@@ -116,9 +116,8 @@ public:
         );
 
 
-    // Destructor
-
-        virtual ~cyclicPointPatch();
+    //- Destructor
+    virtual ~cyclicPointPatch();
 
 
     // Member Functions

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;
+            }
 };
 
 

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
         {

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

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