Update SAMG interfaces

src/foam/matrices/lduMatrix/lduAddressing/lduInterfaces/lduInterface/lduInterface.C

@@ -25,7 +25,7 @@ License
 
 #include "lduInterface.H"
 #include "crMatrix.H"
-
+#include "autoPtr.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -43,7 +43,7 @@ Foam::lduInterface::~lduInterface()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::crMatrix> Foam::lduInterface::prolongationTransfer
+Foam::autoPtr<Foam::crMatrix> Foam::lduInterface::prolongationTransfer
 (
     const Pstream::commsTypes commsType,
     const crMatrix& P
@@ -51,16 +51,15 @@ Foam::tmp<Foam::crMatrix> Foam::lduInterface::prolongationTransfer
 {
     notImplemented
     (
-        "tmp<crMatrix> lduInterface::prolongationTransfer\n"
+        "autoPtr<crMatrix> lduInterface::prolongationTransfer\n"
         "(\n"
-        "    const Pstream::commsTypes commsType,\n"
+        "    const Pstream::commsTypes commsType\n"
-        "    const crMatrix& P\n"
         ") const for type " +
         this->type()
     );
 
     // Dummy return to make the compiler happy
-    return P;
+    return autoPtr<crMatrix>(new crMatrix(1, 1, labelList(1, 0)));
 }
 
 

src/foam/matrices/lduMatrix/lduAddressing/lduInterfaces/lduInterface/lduInterface.H

@@ -49,6 +49,9 @@ namespace Foam
 // Forward declaration of classes
 class crMatrix;
 
+template<class T>
+class autoPtr;
+
 /*---------------------------------------------------------------------------*\
                      Class lduInterface Declaration
 \*---------------------------------------------------------------------------*/
@@ -135,16 +138,18 @@ public:
             virtual void initProlongationTransfer
             (
                 const Pstream::commsTypes commsType,
-                const crMatrix& iF
+                const crMatrix& filteredP
             ) const
-            {}
+            {
+                Pout<< "lduInterface::initProlongationTransfer()" << endl;
+            }
 
             //- Transfer and return prolongation matrix adjacent to
             //  the interface
-            virtual tmp<crMatrix> prolongationTransfer
+            virtual autoPtr<crMatrix> prolongationTransfer
             (
                 const Pstream::commsTypes commsType,
-                const crMatrix& iF
+                const crMatrix& filteredP
             ) const;
 };
 

src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/SAMGInterface/SAMGInterface.C

@@ -54,9 +54,10 @@ Foam::tmp<Foam::scalarField> Foam::SAMGInterface::selectCoeffs
     scalarField& coarseCoeffs = tcoarseCoeffs();
 
     // Added weights to account for non-integral matching
-    forAll (fineAddressing_, ffI)
+    forAll (restrictAddressing_, ffi)
     {
-        coarseCoeffs[ffI] = fineWeights_[ffI]*fineCoeffs[fineAddressing_[ffI]];
+        coarseCoeffs[restrictAddressing_[ffi]] +=
+            restrictWeights_[ffi]*fineCoeffs[fineAddressing_[ffi]];
     }
 
     return tcoarseCoeffs;

src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/SAMGInterface/SAMGInterface.H

@@ -43,6 +43,7 @@ SourceFiles
 #include "autoPtr.H"
 #include "lduPrimitiveMesh.H"
 #include "crMatrix.H"
+#include "labelPair.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,15 +74,29 @@ protected:
 
     // Protected data
 
-        //- Local face-cell addressing.  Contains local coarse level addressing
+        //- Face-cell addressing.  Contains coarse level addressing
-        //  Detected as coarse cells on local side
         labelField faceCells_;
 
-        //- Fine addressing.  Contains fine index for each coarse face
+        //- Fine addressing
+        //  On SAMG interfaces, a single fine coefficient may contribute to
+        //  multiple coarse coefficients using different weights.
+        //  To handle this, a fine coefficient may be visited multiple times
+        //  which is recorded in fineAddressing.
+        //  For simple (matching) interfaces, fineAddressing_[i] = i
+        //  HJ, 21/Jun/2016
         labelField fineAddressing_;
 
-        //- Fine weights
+        //- Restrict addressing
-        scalarField fineWeights_;
+        //  For each fine coefficient, list coarse cluster index it will be
+        //  agglomerated into
+        //  For cases where the fineAddressing is used, restrict addressing
+        //  and weights are expanded to match multiple hits for a single
+        //  fine coefficient, as dictated by fineAddressing
+        //  HJ, 21/Jun/2016
+        labelField restrictAddressing_;
+
+        //- Fine level agglomeration weights
+        scalarField restrictWeights_;
 
 
 private:
@@ -171,12 +186,19 @@ public:
                 return lduMesh_;
             }
 
-            //- Return reference to prolongation matrix
+            //- Return reference to filtered prolongation matrix
             const crMatrix& prolongation() const
             {
                 return prolongation_;
             }
 
+            //- Return reference to neighbour's filtered prolongation matrix
+            const crMatrix& nbrInterfaceProlongation() const
+            {
+                return nbrInterfaceProlongation_;
+            }
+
+
             //- Return local size
             virtual label size() const
             {
@@ -201,10 +223,16 @@ public:
                 return fineAddressing_;
             }
 
-            //- Return fine weights
+            //- Return restrict addressing
-            const scalarField& fineWeights() const
+            const labelField& restrictAddressing() const
             {
-                return fineWeights_;
+                return restrictAddressing_;
+            }
+
+            //- Return fine level agglomeration weights
+            const scalarField& restrictWeights() const
+            {
+                return restrictWeights_;
             }
 
             //- Return the interface internal field of the given field
@@ -225,15 +253,15 @@ public:
             virtual void initProlongationTransfer
             (
                 const Pstream::commsTypes commsType,
-                const crMatrix& iF
+                const crMatrix& filteredP
             ) const = 0;
 
             //- Transfer and return prolongation matrix adjacent to
             //  the interface
-            virtual tmp<crMatrix> prolongationTransfer
+            virtual autoPtr<crMatrix> prolongationTransfer
             (
                 const Pstream::commsTypes commsType,
-                const crMatrix& iF
+                const crMatrix& filteredP
             ) const = 0;
 
 

src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/processorSAMGInterface/processorSAMGInterface.C

@@ -57,10 +57,246 @@ Foam::processorSAMGInterface::processorSAMGInterface
     tag_(fineProcInterface_.tag())
 {
     /*
+    // On master:
+    //  - receive slave's filtered prolongation
+    //  - create a transpose of your local prolongation
+    //  - do the triple product: store the coefficients and the addresses
+    //    row = row of restriction, column = column of prolongation
+    //  - the coefficients should be stored row by row, but with randomized
+    //    columns
+    //  - send the addressing to slave to sort its coefficients to match
+    //    the ordering on master (is it possible to sort at the same time?)
+
+    // Needed:
+    //  - labelList masterOwner,
+    //  - labelList masterNeighbour,
+    //  - label nCoarseCoeffs
+
+
+    // Question: Should we put the masterOwner and masterNeighbour into a linked
+    // list which can be searched quickly? In this case we can create the
+    // addressing on master, send and at the same time, calculate the coarse
+    // coefficients on both sides and store into corresponding order
+
+    // Question: How can I save the addressing of the boundary coeffs for
+    // the next level?
+
+    // Question: With such boundary coefficients, my vector-matrix product is a
+    // matrix-matrix product - change multiplication rule in SAMGInterfaceFields
+    // or account for it with addressing?
+//------------------------------------------------------------------------------
+
+    // Count coarse coefficients and create addressing on master
+    if (myProcNo() < neighbProcNo())
+    {
+//------------------------------------------------------------------------------
+//                              COUNT COARSE COEFFS
+//------------------------------------------------------------------------------
+        // Filtered prolongation matrix from my side
+        const labelList& localOwner = prolongation.crAddr().column();
+
+        // Filtered prolongation matrix from my side
+        tmp<crMatrix> tProlongationT(prolongation.T());
+        const labelList& rRowStart = tProlongationT->crAddr().rowStart();
+        const labelList& rColumn = tProlongationT->crAddr().column();
+        const label rNRows = tProlongationT->crAddr().nRows();
+
+        // Filtered prolongation matrix from neighbour - to obtain restriction,
+        // make a transpose
+        const labelList& pRowStart =
+            nbrInterfaceProlongation.crAddr().rowStart();
+        const labelList& pColumn = nbrInterfaceProlongation.crAddr().column();
+        const label pNCols = nbrInterfaceProlongation.crAddr().nCols();
+
+        labelList coeffMark(pNCols, -1);
+        label nCoarseCoeffs = 0;
+
+        // Row of R
+        for (label ir = 0; ir < rNRows; ir++)
+        {
+            //Restart coeffMark for each restriction row
+            coeffMark = -1;
+
+            // Row start addressing R
+            for
+            (
+                label indexR = rRowStart[ir];
+                indexR < rRowStart[ir + 1];
+                indexR++
+            )
+            {
+                // Column of coefficient in R
+                const label jr = rColumn[indexR];
+
+                for
+                (
+                    label indexP = pRowStart[jr];
+                    indexP < pRowStart[jr + 1];
+                    indexP++
+                )
+                {
+                    label jp = pColumn[indexP];
+
+                    if (coeffMark[jp] == -1)
+                    {
+                        // Found a new coarse coeff!
+                        coeffMark[jp] = nCoarseCoeffs;
+                        nCoarseCoeffs++;
+                    }
+                }
+            }
+        }
+
+//------------------------------------------------------------------------------
+//                      CREATE HASH TABLE FOR ADDRESSING
+//------------------------------------------------------------------------------
+
+        HashTable<label, labelPair, Hash<labelPair> > coarseLevel(nCoarseCoeffs);
+        faceCells_.setSize(nCoarseCoeffs);
+
+        // Reset for creating addressing
+        nCoarseCoeffs = 0;
+        coeffMark = -1;
+
+        // Row of R
+        for (label ir = 0; ir < rNRows; ir++)
+        {
+            // Restart coeffMark for each restriction row
+            coeffMark = -1;
+
+            // Row start addressing R
+            for
+            (
+                label indexR = rRowStart[ir];
+                indexR < rRowStart[ir + 1];
+                indexR++
+            )
+            {
+                // Column of coefficient in R
+                const label jr = rColumn[indexR];
+
+                // Grab column of the original prolongation, this is the
+                // faceCell_!
+                const label ja = localOwner[ir];
+
+                for
+                (
+                    label indexP = pRowStart[jr];
+                    indexP < pRowStart[jr + 1];
+                    indexP++
+                )
+                {
+                    label jp = pColumn[indexP];
+
+                    label identify = coeffMark[jp];
+
+                    if (identify == -1)
+                    {
+                        // Found a new coarse coeff!
+                        identify = nCoarseCoeffs;
+                        coeffMark[jp] = nCoarseCoeffs;
+
+                        faceCells_[identify] = ja;
+
+                        // Save address and coeff into HashTable
+                        coarseLevel.insert(labelPair(ir, jp), nCoarseCoeffs);
+                        nCoarseCoeffs++;
+                    }
+                }
+            }
+        }
+
+        // Send to slave
+        OPstream toNbr(Pstream::blocking, neighbProcNo());
+
+        toNbr<< coarseLevel;
+
+    }
+    // Slave
+    else
+    {
+        // Slave must receive the hash table sent by the master to know how to
+        // sort its boundary coefficients
+
+        IPstream fromNbr(Pstream::blocking, neighbProcNo());
+
+        masterCoarseLevel_ =
+            HashTable<label, labelPair, Hash<labelPair> >(fromNbr);
+
+        // Filtered prolongation matrix from my side
+        const labelList& localOwner = prolongation.crAddr().column();
+
+        // Filtered prolongation matrix from my side
+        tmp<crMatrix> tProlongationT(prolongation.T());
+        const labelList& rRowStart = tProlongationT->crAddr().rowStart();
+        const labelList& rColumn = tProlongationT->crAddr().column();
+        const label rNRows = tProlongationT->crAddr().nRows();
+
+        // Filtered prolongation matrix from neighbour - to obtain restriction,
+        // make
+        // a transpose
+        const labelList& pRowStart = nbrInterfaceProlongation.crAddr().rowStart();
+        const labelList& pColumn = nbrInterfaceProlongation.crAddr().column();
+        const label pNCols = nbrInterfaceProlongation.crAddr().nCols();
+
+        faceCells_.setSize(masterCoarseLevel_.size());
+
+        labelList coeffMark(pNCols, -1);
+
+        // Row of R
+        for (label ir = 0; ir < rNRows; ir++)
+        {
+            // Restart coeffMark for each restriction row
+            coeffMark = -1;
+
+            // Row start addressing R
+            for
+            (
+                label indexR = rRowStart[ir];
+                indexR < rRowStart[ir + 1];
+                indexR++
+            )
+            {
+                // Column of coefficient in R
+                const label jr = rColumn[indexR];
+
+                // Grab column of the original prolongation -  this is the
+                // faceCell_!!!
+                label ja = localOwner[ir];
+
+                for
+                (
+                    label indexP = pRowStart[jr];
+                    indexP < pRowStart[jr + 1];
+                    indexP++
+                )
+                {
+                    label jp = pColumn[indexP];
+
+                    label identify = coeffMark[jp];
+
+                    if (identify == -1)
+                    {
+                        label address = masterCoarseLevel_[labelPair(jp,ja)];
+
+                        // Found a new coarse coeff!
+                        identify = address;
+
+                        coeffMark[jp] = address;
+                        faceCells_[identify] = ja;
+                    }
+                }
+            }
+        }
+    }
+
+
     // Analyse the local and neighbour row label:
     //  local coarse, remote coarse = regular coarse face
     //  local coarse, remote fine = local expanded face: receive prolonged data
     //  local fine, remote coarse = neighbour expanded face: send prolonged data
+    //  local fine, remote fine = local and neighbour expanded face: send and
+    //                            receive
 
     // Algorithm
     // Go through local row labels and examine cases
@@ -89,6 +325,19 @@ Foam::processorSAMGInterface::processorSAMGInterface
         Foam::max(128, fineProcInterface_.interfaceSize()/4)
     );
 
+    // For FINE to FINE communication - local boundary coefficient multiplyed by
+    // local weight - send to the other side for triple product
+    HashTable<labelList, label, Hash<label> > neighboursFromLocalFineToFine
+    (
+        Foam::max(128, fineProcInterface_.interfaceSize()/4)
+    );
+
+    HashTable<scalarField, label, Hash<label> > weightsFromLocalFineToFine
+    (
+        Foam::max(128, fineProcInterface_.interfaceSize()/4)
+    );
+
+
     // Get access to the prolongation addressing and coefficients
     const labelList& pRowStart = prolongation.crAddr().rowStart();
     const labelList& pColumn = prolongation.crAddr().column();
@@ -98,6 +347,7 @@ Foam::processorSAMGInterface::processorSAMGInterface
     const labelList& fineFaceCells = fineInterface.faceCells();
 
     // Collect local fine to neighbour coarse connections for communication
+    // TU, 19 May 2017: and fine to fine (communication through triple product)
     forAll (localRowLabel, faceI)
     {
         if
@@ -127,11 +377,45 @@ Foam::processorSAMGInterface::processorSAMGInterface
             neighboursFromLocalFine.insert(neighbourRowLabel[faceI], nbrs);
             weightsFromLocalFine.insert(neighbourRowLabel[faceI], weights);
         }
+        // FINE to FINE communication - TU, May 2017
+        if
+        (
+            localRowLabel[faceI] < 0
+         && neighbourRowLabel[faceI] < 0
+        )
+        {
+            // Found local FINE to neighbour FINE interface
+
+            // First, multiply all local interface boundary coefficients with
+            // local prolongation weights - you will send this to the other side
+            // in localBoundaryProlongation
+
+            // Collect local coarse neighbours and weights from the
+            // prolongation matrix to send to other processor
+            const label curStart = pRowStart[fineFaceCells[faceI]];
+            const label curEnd = pRowStart[fineFaceCells[faceI] + 1];
+            const label nCoeffs = curEnd - curStart;
+
+            labelList nbrs(nCoeffs);
+            scalarField weights(nCoeffs);
+
+            forAll (nbrs, i)
+            {
+                nbrs[i] = pColumn[curStart + i];
+                weights[i] = pCoeffs[curStart + i];
+            }
+
+            // Insert neighbours under remote coarse index
+            neighboursFromLocalFineToFine.insert(neighbourRowLabel[faceI], nbrs);
+            weightsFromLocalFineToFine.insert(neighbourRowLabel[faceI], weights);
+        }
     }
 
     // Receive remote prolongation data
     HashTable<labelList, label, Hash<label> > neighboursFromRemoteFine;
     HashTable<scalarField, label, Hash<label> > weightsFromRemoteFine;
+    HashTable<labelList, label, Hash<label> > neighboursFromRemoteFineToFine;
+    HashTable<scalarField, label, Hash<label> > weightsFromRemoteFineToFine;
 
     // Send and receive the addressing from the other side
     {
@@ -147,6 +431,13 @@ Foam::processorSAMGInterface::processorSAMGInterface
 
         weightsFromRemoteFine =
             HashTable<scalarField, label, Hash<label> >(fromNbr);
+
+        neighboursFromRemoteFineToFine =
+            HashTable<labelList, label, Hash<label> >(fromNbr);
+
+        weightsFromRemoteFineToFine =
+            HashTable<scalarField, label, Hash<label> >(fromNbr);
+
     }
 
     // Assemble connectivity
@@ -232,7 +523,36 @@ Foam::processorSAMGInterface::processorSAMGInterface
         }
         else
         {
-            // Fine to fine.  Ignore
+            // Fine to fine.
+            // Establish connection as in triple product: multiply restriction
+            // matrix on my side, send to the other side and multiply with
+            // prolongation there
+
+            // Get fine faceCells
+            const scalarField& fineLocalCoeffs =
+            // Pick up local prolongation coarse entries and for all
+            // add a new face with appropriate weight
+            // Note: faceCells changes, but how???
+
+            const labelList& curLocalCoarseNbrs =
+                neighboursFromLocalFineToFine[neighbourRowLabel[faceI]];
+
+            const scalarField& curLocalCoarseWeights =
+                weightsFromLocalFineToFine[neighbourRowLabel[faceI]];
+
+            forAll (curLocalCoarseNbrs, curNbrI)
+            {
+                // Create new coarse face
+                faceCells_[nCoarseFaces] = curLocalCoarseNbrs[curNbrI];
+                fineAddressing_[nCoarseFaces] = faceI;
+                fineWeights_[nCoarseFaces] =
+                    curLocalCoarseWeights[curNbrI]*fineLocalCoeffs[faceI];
+                nCoarseFaces++;
+
+                // Store fineFaceCells, recording the fine index of prolonged
+                // cell locally and all coarse neighbours on the other side
+                // HJ, HERE!!!
+            }
         }
     }
 
@@ -305,28 +625,27 @@ Foam::tmp<Foam::labelField> Foam::processorSAMGInterface::internalFieldTransfer
 void Foam::processorSAMGInterface::initProlongationTransfer
 (
     const Pstream::commsTypes commsType,
-    const crMatrix& P
+    const crMatrix& filteredP
 ) const
 {
-    // Select the part of the prolongation matrix to send
+    // Send prolongation matrix, using IOstream operators
-
+    OPstream toNbr(Pstream::blocking, neighbProcNo());
-    // Send prolongation matrix
+    toNbr<< filteredP;
-    Pout<< "HJ, in send" << endl;
 }
 
-Foam::tmp<Foam::crMatrix> Foam::processorSAMGInterface::prolongationTransfer
+
+Foam::autoPtr<Foam::crMatrix>
+Foam::processorSAMGInterface::prolongationTransfer
 (
     const Pstream::commsTypes commsType,
-    const crMatrix& P
+    const crMatrix& filteredP
 ) const
 {
-    // Receive and return prolongation matrix
+    IPstream fromNbr(Pstream::blocking, neighbProcNo());
-    Pout<< "HJ, in receive" << endl;
 
-    // Dummy
+    autoPtr<crMatrix> tnbrP(new crMatrix(fromNbr));
-    tmp<crMatrix> tcr(new crMatrix(5, 5, labelList(5)));
 
-    return tcr;
+    return tnbrP;
 }
 
 

src/foam/matrices/lduMatrix/solvers/AMG/interfaces/SAMGInterfaces/processorSAMGInterface/processorSAMGInterface.H

@@ -165,15 +165,15 @@ public:
             virtual void initProlongationTransfer
             (
                 const Pstream::commsTypes commsType,
-                const crMatrix& iF
+                const crMatrix& P
             ) const;
 
             //- Transfer and return prolongation matrix adjacent to
             //  the interface
-            virtual tmp<crMatrix> prolongationTransfer
+            virtual autoPtr<crMatrix> prolongationTransfer
             (
                 const Pstream::commsTypes commsType,
-                const crMatrix& iF
+                const crMatrix& P
             ) const;