First version of block matrix selection and interfaces

src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGCycle.C

@@ -99,7 +99,7 @@ void Foam::BlockAMGCycle<Type>::makeCoarseLevels(const label nMaxLevels)
             }
         }
 
-        if (BlockLduMatrix<Type>::debug >= 2)
+        if (blockLduMatrix::debug >= 2)
         {
             Info<< "Created " << nLevels_ << " AMG levels" << endl;
         }

src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGInterfaceFields/BlockAMGInterfaceField/newBlockAMGInterfaceField.C

@@ -45,9 +45,11 @@ Foam::BlockAMGInterfaceField<Type>::New
     {
         FatalErrorIn
         (
-            "BlockAMGInterfaceField::New"
-            "(const AMGInterface& AMGCp, "
-            "const BlockLduInterfaceField<Type>& fineInterface)"
+            "BlockAMGInterfaceField::New\n"
+            "(\n"
+            "    const AMGInterface& AMGCp,\n"
+            "    const BlockLduInterfaceField<Type>& fineInterface\n"
+            ")"
         )   << "Unknown BlockAMGInterfaceField type " << coupleType << ".\n"
             << "Valid BlockAMGInterfaceField types are :"
             << lduInterfaceConstructorTablePtr_->sortedToc()

src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixAgglomeration/BlockMatrixAgglomeration.C

@@ -366,7 +366,7 @@ void Foam::BlockMatrixAgglomeration<Type>::calcAgglomeration()
 
     reduce(coarsen_, andOp<bool>());
 
-    if (BlockLduMatrix<Type>::debug >= 3)
+    if (blockLduMatrix::debug >= 3)
     {
         // Count singleton clusters
         label nSingleClusters = 0;
@@ -893,17 +893,14 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
     // Create coarse-level coupled interfaces
 
     // Create coarse interfaces, addressing and coefficients
-    const label interfaceSize =
-        const_cast<BlockLduMatrix<Type>& >(matrix_).interfaces().size();
-
     const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
         interfaceFields =
         const_cast<BlockLduMatrix<Type>&>(matrix_).interfaces();
 
     // Set the coarse interfaces and coefficients
-    lduInterfacePtrsList coarseInterfaces(interfaceSize);
+    lduInterfacePtrsList coarseInterfaces(interfaceFields.size());
 
-    labelListList coarseInterfaceAddr(interfaceSize);
+    labelListList coarseInterfaceAddr(interfaceFields.size());
 
     // Add the coarse level
 
@@ -921,7 +918,7 @@ Foam::BlockMatrixAgglomeration<Type>::restrictMatrix() const
     );
 
     // Initialise transfer of restrict addressing on the interface
-    // HJ, consider blocking comms.  HJ, 9/Jun/2016
+    // HJ, must use blocking comms.  HJ, 9/Jun/2016
     forAll (interfaceFields, intI)
     {
         if (interfaceFields.set(intI))

src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixClustering/BlockMatrixClustering.C

@@ -495,7 +495,7 @@ void Foam::BlockMatrixClustering<Type>::calcClustering()
 
     reduce(coarsen_, andOp<bool>());
 
-    if (BlockLduMatrix<Type>::debug >= 3)
+    if (blockLduMatrix::debug >= 3)
     {
         // Count singleton clusters
         label nSingleClusters = 0;
@@ -1052,7 +1052,7 @@ Foam::BlockMatrixClustering<Type>::restrictMatrix() const
     );
 
     // Initialise transfer of restrict addressing on the interface
-    // HJ, consider blocking comms.  HJ, 9/Jun/2016
+    // HJ, must use blocking comms.  HJ, 9/Jun/2016
     forAll (interfaceFields, intI)
     {
         if (interfaceFields.set(intI))

src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.C

@@ -60,7 +60,7 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
 //------------------------------------------------------------------------------
 //             MATRIX DATA: ADDRESSING, COEFFICIENTS, COEFF NORMS
 //------------------------------------------------------------------------------
-
+    Info<< "Start equation selection" << endl;
     // Get addressing
     const unallocLabelList& rowStart = matrix_.lduAddr().ownerStartAddr();
     const unallocLabelList& losortAddr = matrix_.lduAddr().losortAddr();
@@ -632,6 +632,8 @@ void Foam::BlockMatrixSelection<Type>::calcCoarsening()
         // Coarsening did not succeed.  Delete Pptr
         deleteDemandDrivenData(Pptr_);
     }
+
+    Info<< "End equation selection" << endl;
 }
 
 
@@ -675,6 +677,7 @@ template<class Type>
 Foam::autoPtr<Foam::BlockAMGLevel<Type> >
 Foam::BlockMatrixSelection<Type>::restrictMatrix() const
 {
+    Info<< "Start matrix restriction" << endl;
     if (!coarsen_)
     {
         FatalErrorIn("autoPtr<amgMatrix> samgPolicy::restrictMatrix() const")
@@ -702,8 +705,11 @@ Foam::BlockMatrixSelection<Type>::restrictMatrix() const
      || crR.nRows() != nCoarseEqns_
     )
     {
-        FatalErrorIn("")
-            << "Incompatible matrices for triple product: "
+        FatalErrorIn
+        (
+            "autoPtr<Foam::BlockAMGLevel<Type> >"
+            "BlockMatrixSelection<Type>::restrictMatrix() const"
+        )   << "Incompatible matrices for triple product: "
             << "R( " << crR.nRows() << " ," << crR.nCols() << ") "
             << "A( " << nEqns << " ," << nEqns << ") "
             << "P( " << crP.nRows() << " ," << crP.nCols() << ") "
@@ -989,6 +995,79 @@ Foam::BlockMatrixSelection<Type>::restrictMatrix() const
         )
     );
 
+    const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
+        interfaceFields =
+        const_cast<BlockLduMatrix<Type>&>(matrix_).interfaces();
+
+    // Set the coarse interfaces and coefficients
+    lduInterfacePtrsList coarseInterfaces(interfaceFields.size());
+
+    labelListList coarseInterfaceAddr(interfaceFields.size());
+
+    // Initialise transfer of restrict addressing on the interface
+    // HJ, reconsider blocking comms.  HJ, 9/Jun/2016
+    forAll (interfaceFields, intI)
+    {
+        if (interfaceFields.set(intI))
+        {
+            interfaceFields[intI].coupledInterface().initInternalFieldTransfer
+            (
+                Pstream::blocking,
+                rowLabel_
+            );
+        }
+    }
+
+    // Store coefficients to avoid tangled communications
+    // HJ, 1/Apr/2009
+    FieldField<Field, label> fineInterfaceAddr(interfaceFields.size());
+
+    forAll (interfaceFields, intI)
+    {
+        if (interfaceFields.set(intI))
+        {
+            const lduInterface& fineInterface =
+                interfaceFields[intI].coupledInterface();
+
+            fineInterfaceAddr.set
+            (
+                intI,
+                new labelField
+                (
+                    fineInterface.internalFieldTransfer
+                    (
+                        Pstream::blocking,
+                        rowLabel_
+                    )
+                )
+            );
+        }
+    }
+
+    // Create AMG interfaces
+    forAll (interfaceFields, intI)
+    {
+        if (interfaceFields.set(intI))
+        {
+            const lduInterface& fineInterface =
+                interfaceFields[intI].coupledInterface();
+
+            coarseInterfaces.set
+            (
+                intI,
+                AMGInterface::New
+                (
+                    coarseAddrPtr(),
+                    coarseInterfaces,
+                    fineInterface,
+                    fineInterface.interfaceInternalField(rowLabel_),
+                    fineInterfaceAddr[intI]
+                ).ptr()
+            );
+        }
+    }
+
+    //HJ: Add interface fields
 //------------------------------------------------------------------------------
 //                            CREATE COARSE MATRIX
 //------------------------------------------------------------------------------
@@ -1274,7 +1353,8 @@ Foam::BlockMatrixSelection<Type>::restrictMatrix() const
         )   << "Matrix diagonal of scalar or linear type not implemented"
             << abort(FatalError);
     }
-
+    Info<< "End matrix restriction.  Level size: " << nCoarseEqns_
+        << endl;
     // Create and return BlockAMGLevel
     return autoPtr<BlockAMGLevel<Type> >
     (

src/foam/matrices/blockLduMatrix/BlockAMG/BlockMatrixCoarsening/BlockMatrixSelection/BlockMatrixSelection.H

@@ -81,6 +81,7 @@ class BlockMatrixSelection
         crMatrix* Rptr_;
 
         //- Coarsening array - array with labels of coarse/fine equations
+        //  (coarse = coarseID / fine = -1)
         labelList rowLabel_;
 
 
@@ -101,9 +102,9 @@ class BlockMatrixSelection
         // Equation type
         enum equationType
         {
-            UNDECIDED = -1,
-            COARSE = 0,
-            FINE = 1
+            UNDECIDED = -2,
+            FINE = -1,
+            COARSE = 0     // Later replaced with coarse matrix equation index
         };
 
 
@@ -136,8 +137,7 @@ public:
 
     // Member Functions
 
-        //- Return array with labels of equations (coarse = 0 / fine = 1)
-        //- (used only for postprocessing)
+        //- Return array with labels of equations
         const labelList& coarseningLabels() const
         {
             return rowLabel_;

src/foam/matrices/blockLduMatrix/BlockAMG/coarseBlockAMGLevel.C

@@ -235,15 +235,15 @@ void Foam::coarseBlockAMGLevel<Type>::solve
     }
 
     // Switch of debug in top-level direct solve
-    label oldDebug = BlockLduMatrix<Type>::debug();
+    label oldDebug = blockLduMatrix::debug();
 
-    if (BlockLduMatrix<Type>::debug >= 4)
+    if (blockLduMatrix::debug >= 4)
     {
-        BlockLduMatrix<Type>::debug = 1;
+        blockLduMatrix::debug = 1;
     }
     else
     {
-        BlockLduMatrix<Type>::debug = 0;
+        blockLduMatrix::debug = 0;
     }
 
     if (matrixPtr_->symmetric())
@@ -274,7 +274,7 @@ void Foam::coarseBlockAMGLevel<Type>::solve
     }
 
     // Restore debug
-    BlockLduMatrix<Type>::debug = oldDebug;
+    blockLduMatrix::debug = oldDebug;
 
     // Escape cases of top-level solver divergence
     if
@@ -294,7 +294,7 @@ void Foam::coarseBlockAMGLevel<Type>::solve
         coarseSolverPerf.print();
     }
 
-    if (BlockLduMatrix<Type>::debug >= 3)
+    if (blockLduMatrix::debug >= 3)
     {
         coarseSolverPerf.print();
     }

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockIterativeSolver/BlockIterativeSolver.C

@@ -79,7 +79,7 @@ Foam::scalar Foam::BlockIterativeSolver<Type>::normFactor
 
     scalar normFactor = gSum(mag(wA - pA) + mag(b - pA)) + this->small_;
 
-    if (BlockLduMatrix<Type>::debug >= 2)
+    if (blockLduMatrix::debug >= 2)
     {
         Info<< "Iterative solver normalisation factor = "
             << normFactor << endl;

src/foam/matrices/blockLduMatrix/BlockLduSolvers/Segregated/SegregatedSolver.C

@@ -109,7 +109,7 @@ Foam::BlockSolverPerformance<Type> Foam::SegregatedSolver<Type>::solve
 
     if (blockMatrix.componentCoupled())
     {
-        if (BlockLduMatrix<Type>::debug >= 2)
+        if (blockLduMatrix::debug >= 2)
         {
             Info << " Component coupled segregation" << endl;
         }

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

@@ -72,6 +72,7 @@ public:
     //- Runtime type information
     TypeName("ggiLduInterface");
 
+
     // Constructors
 
         //- Construct null

src/foam/matrices/lduMatrix/solvers/AMG/interfaceFields/AMGInterfaceFields/AMGInterfaceField/newAMGInterfaceField.C

@@ -42,9 +42,11 @@ Foam::autoPtr<Foam::AMGInterfaceField> Foam::AMGInterfaceField::New
     {
         FatalErrorIn
         (
-            "AMGInterfaceField::New"
-            "(const AMGInterface& AMGCp, "
-            "const lduInterfaceField& fineInterface)"
+            "AMGInterfaceField::New\n"
+            "(\n"
+            "    const AMGInterface& AMGCp,\n"
+            "    const lduInterfaceField& fineInterface\n"
+            ")"
         )   << "Unknown AMGInterfaceField type " << coupleType << ".\n"
             << "Valid AMGInterfaceField types are :"
             << lduInterfaceConstructorTablePtr_->sortedToc()

src/foam/matrices/lduMatrix/solvers/AMG/interfaces/AMGInterfaces/AMGInterface/AMGInterface.H

@@ -162,10 +162,7 @@ public:
 
         //- Construct from fine-level interface,
         //  local and neighbour restrict addressing
-        AMGInterface
-        (
-            const lduPrimitiveMesh& lduMesh
-        )
+        AMGInterface(const lduPrimitiveMesh& lduMesh)
         :
             lduMesh_(lduMesh)
         {}