Cached block AMG: Work in progress

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

@@ -107,6 +107,22 @@ void Foam::BlockAMGCycle<Type>::makeCoarseLevels(const label nMaxLevels)
 }
 
 
+template<class Type>
+Foam::BlockLduMatrix<Type>& Foam::BlockAMGCycle<Type>::coarseMatrix()
+{
+    if (!coarseLevelPtr_)
+    {
+        FatalErrorIn
+        (
+            "BlockLduMatrix<Type>& BlockAMGCycle<Type>::coarseMatrix()"
+        )   << "Coarse level not available"
+            << abort(FatalError);
+    }
+
+    return coarseLevelPtr_->matrix();
+}
+
+
 template<class Type>
 void Foam::BlockAMGCycle<Type>::fixedCycle
 (
@@ -193,7 +209,21 @@ void Foam::BlockAMGCycle<Type>::fixedCycle
     else
     {
         // Call direct solver
-        levelPtr_->solve(x, b, 1e-9, 0);
+        levelPtr_->solve(x, b, 1e-7, 0);
+    }
+}
+
+
+template<class Type>
+void Foam::BlockAMGCycle<Type>::initMatrix()
+{
+    // Update current level
+    levelPtr_->initLevel(coarseLevelPtr_);
+
+    // If present, update coarse levels recursively
+    if (coarseLevelPtr_)
+    {
+        coarseLevelPtr_->initMatrix();
     }
 }
 

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

@@ -143,6 +143,9 @@ public:
             return nLevels_;
         }
 
+        //- Return coarse matrix
+        BlockLduMatrix<Type>& coarseMatrix();
+
         //- Calculate residual
         void residual
         (
@@ -166,6 +169,9 @@ public:
             const label nPostSweeps,
             const bool scale
         ) const;
+
+        //- Re-initialise preconditioner after matrix coefficient update
+        void initMatrix();
 };
 
 

src/foam/matrices/blockLduMatrix/BlockAMG/BlockAMGLevel.H

@@ -86,6 +86,9 @@ public:
 
     // Member Functions
 
+        //- Return reference to matrix
+        virtual BlockLduMatrix<Type>& matrix() = 0;
+
         //- Return reference to x
         virtual Field<Type>& x() = 0;
 
@@ -144,6 +147,12 @@ public:
 
         //- Create next level from current level
         virtual autoPtr<BlockAMGLevel> makeNextLevel() const = 0;
+
+        //- Re-initialise AMG level after matrix coefficient update
+        virtual void initLevel
+        (
+            autoPtr<Foam::BlockAMGLevel<Type> >& coarseLevelPtr
+        ) = 0;
 };
 
 

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

@@ -192,6 +192,9 @@ public:
             Field<Type>& x,
             const Field<Type>& coarseX
         ) const;
+
+        //- Update coarse matrix using same coefficients
+        virtual void updateMatrix(BlockLduMatrix<Type>& coarseMatrix) const;
 };
 
 

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

@@ -495,6 +495,163 @@ void Foam::BlockMatrixClustering<Type>::calcClustering()
 
     reduce(coarsen_, andOp<bool>());
 
+    // If the matrix will be coarsened, create off-diagonal agglomeration
+        // Does the matrix have solo equations
+    bool soloEqns = nSolo_ > 0;
+
+    // Storage for block neighbours and coefficients
+
+    // Guess initial maximum number of neighbours in block
+    label maxNnbrs = 10;
+
+    // Number of neighbours per block
+    labelList blockNnbrs(nCoarseEqns_, 0);
+
+    // Setup initial packed storage for neighbours and coefficients
+    labelList blockNbrsData(maxNnbrs*nCoarseEqns_);
+
+    const label nFineCoeffs = upperAddr.size();
+
+    // Create face-restriction addressing
+    coeffRestrictAddr_.setSize(nFineCoeffs);
+
+    // Initial neighbour array (not in upper-triangle order)
+    labelList initCoarseNeighb(nFineCoeffs);
+
+    // Counter for coarse coeffs
+    label nCoarseCoeffs = 0;
+
+    // Loop through all fine coeffs
+    forAll (upperAddr, fineCoeffi)
+    {
+        label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffi]];
+        label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffi]];
+
+        // If the coefficient touches block zero and solo equations are
+        // present, skip it
+        if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0))
+        {
+            continue;
+        }
+
+        if (rmUpperAddr == rmLowerAddr)
+        {
+            // For each fine coeff inside of a coarse cluster keep the address
+            // of the cluster corresponding to the coeff in the
+            // coeffRestrictAddr_ as a negative index
+            coeffRestrictAddr_[fineCoeffi] = -(rmUpperAddr + 1);
+        }
+        else
+        {
+            // This coeff is a part of a coarse coeff
+
+            label cOwn = rmUpperAddr;
+            label cNei = rmLowerAddr;
+
+            // Get coarse owner and neighbour
+            if (rmUpperAddr > rmLowerAddr)
+            {
+                cOwn = rmLowerAddr;
+                cNei = rmUpperAddr;
+            }
+
+            // Check the neighbour to see if this coeff has already been found
+            bool nbrFound = false;
+            label& ccnCoeffs = blockNnbrs[cOwn];
+
+            for (int i = 0; i < ccnCoeffs; i++)
+            {
+                if (initCoarseNeighb[blockNbrsData[maxNnbrs*cOwn + i]] == cNei)
+                {
+                    nbrFound = true;
+                    coeffRestrictAddr_[fineCoeffi] =
+                        blockNbrsData[maxNnbrs*cOwn + i];
+                    break;
+                }
+            }
+
+            if (!nbrFound)
+            {
+                if (ccnCoeffs >= maxNnbrs)
+                {
+                    // Double the size of list and copy data
+                    label oldMaxNnbrs = maxNnbrs;
+                    maxNnbrs *= 2;
+
+                    // Resize and copy list
+                    const labelList oldBlockNbrsData = blockNbrsData;
+                    blockNbrsData.setSize(maxNnbrs*nCoarseEqns_);
+
+                    forAll (blockNnbrs, i)
+                    {
+                        for (int j = 0; j < blockNnbrs[i]; j++)
+                        {
+                            blockNbrsData[maxNnbrs*i + j] =
+                                oldBlockNbrsData[oldMaxNnbrs*i + j];
+                        }
+                    }
+                }
+
+                blockNbrsData[maxNnbrs*cOwn + ccnCoeffs] = nCoarseCoeffs;
+                initCoarseNeighb[nCoarseCoeffs] = cNei;
+                coeffRestrictAddr_[fineCoeffi] = nCoarseCoeffs;
+                ccnCoeffs++;
+
+                // New coarse coeff created
+                nCoarseCoeffs++;
+            }
+        }
+    } // End for all fine coeffs
+
+
+    // Renumber into upper-triangular order
+
+    // All coarse owner-neighbour storage
+    labelList coarseOwner(nCoarseCoeffs);
+    labelList coarseNeighbour(nCoarseCoeffs);
+    labelList coarseCoeffMap(nCoarseCoeffs);
+
+    label coarseCoeffi = 0;
+
+    forAll (blockNnbrs, cci)
+    {
+        label* cCoeffs = &blockNbrsData[maxNnbrs*cci];
+        label ccnCoeffs = blockNnbrs[cci];
+
+        for (int i = 0; i < ccnCoeffs; i++)
+        {
+            coarseOwner[coarseCoeffi] = cci;
+            coarseNeighbour[coarseCoeffi] = initCoarseNeighb[cCoeffs[i]];
+            coarseCoeffMap[cCoeffs[i]] = coarseCoeffi;
+            coarseCoeffi++;
+        }
+    }
+
+    forAll(coeffRestrictAddr_, fineCoeffi)
+    {
+        label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffi]];
+        label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffi]];
+
+        // If the coefficient touches block zero and solo equations are
+        // present, skip it
+        if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0))
+        {
+            continue;
+        }
+
+        if (coeffRestrictAddr_[fineCoeffi] >= 0)
+        {
+            coeffRestrictAddr_[fineCoeffi] =
+                coarseCoeffMap[coeffRestrictAddr_[fineCoeffi]];
+        }
+    }
+
+    // Clear the temporary storage for the coarse matrix data
+    blockNnbrs.setSize(0);
+    blockNbrsData.setSize(0);
+    initCoarseNeighb.setSize(0);
+    coarseCoeffMap.setSize(0);
+  
     if (blockLduMatrix::debug >= 3)
     {
         // Count singleton clusters
@@ -555,10 +712,8 @@ void Foam::BlockMatrixClustering<Type>::restrictDiag
         squareTypeField& activeCoarseCoeff = coarseCoeff.asSquare();
         const squareTypeField& activeCoeff = Coeff.asSquare();
 
-        forAll (coarseCoeff, i)
+        // Reset coefficients to zero
-        {
+        activeCoarseCoeff = pTraits<squareType>::zero;
-            activeCoarseCoeff[i] = pTraits<squareType>::zero;
-        }
 
         forAll (Coeff, i)
         {
@@ -577,10 +732,8 @@ void Foam::BlockMatrixClustering<Type>::restrictDiag
         linearTypeField& activeCoarseCoeff = coarseCoeff.asLinear();
         const linearTypeField& activeCoeff = Coeff.asLinear();
 
-        forAll (coarseCoeff, i)
+        // Reset coefficients to zero
-        {
+        activeCoarseCoeff = pTraits<linearType>::zero;
-            activeCoarseCoeff[i] = pTraits<linearType>::zero;
-        }
 
         forAll (Coeff, i)
         {
@@ -599,10 +752,8 @@ void Foam::BlockMatrixClustering<Type>::restrictDiag
         scalarTypeField& activeCoarseCoeff = coarseCoeff.asScalar();
         const scalarTypeField& activeCoeff = Coeff.asScalar();
 
-        forAll (coarseCoeff, i)
+        // Reset coefficients to zero
-        {
+        activeCoarseCoeff = pTraits<scalarType>::zero;
-            activeCoarseCoeff[i] = pTraits<scalarType>::zero;
-        }
 
         forAll (Coeff, i)
         {
@@ -624,7 +775,6 @@ template<class Type>
 template<class DiagType, class ULType>
 void Foam::BlockMatrixClustering<Type>::agglomerateCoeffs
 (
-    const labelList& coeffRestrictAddr,
     Field<DiagType>& activeCoarseDiag,
     Field<ULType>& activeCoarseUpper,
     const Field<ULType>& activeFineUpper,
@@ -638,7 +788,10 @@ void Foam::BlockMatrixClustering<Type>::agglomerateCoeffs
     const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
     const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
 
-    forAll(coeffRestrictAddr, fineCoeffI)
+    // Reset coefficients to zero.  Cannot touch the diagonal
+    activeCoarseUpper = pTraits<ULType>::zero;
+    
+    forAll(coeffRestrictAddr_, fineCoeffI)
     {
         label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffI]];
         label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffI]];
@@ -650,7 +803,7 @@ void Foam::BlockMatrixClustering<Type>::agglomerateCoeffs
             continue;
         }
 
-        label cCoeff = coeffRestrictAddr[fineCoeffI];
+        label cCoeff = coeffRestrictAddr_[fineCoeffI];
 
         if (cCoeff >= 0)
         {
@@ -673,7 +826,6 @@ template<class Type>
 template<class DiagType, class ULType>
 void Foam::BlockMatrixClustering<Type>::agglomerateCoeffs
 (
-    const labelList& coeffRestrictAddr,
     Field<DiagType>& activeCoarseDiag,
     Field<ULType>& activeCoarseUpper,
     const Field<ULType>& activeFineUpper,
@@ -688,7 +840,11 @@ void Foam::BlockMatrixClustering<Type>::agglomerateCoeffs
     const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
     const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
 
-    forAll(coeffRestrictAddr, fineCoeffI)
+    // Reset coefficients to zero.  Cannot touch the diagonal
+    activeCoarseUpper = pTraits<ULType>::zero;
+    activeCoarseLower = pTraits<ULType>::zero;
+    
+    forAll(coeffRestrictAddr_, fineCoeffI)
     {
         label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffI]];
         label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffI]];
@@ -700,7 +856,7 @@ void Foam::BlockMatrixClustering<Type>::agglomerateCoeffs
             continue;
         }
 
-        label cCoeff = coeffRestrictAddr[fineCoeffI];
+        label cCoeff = coeffRestrictAddr_[fineCoeffI];
 
         if (cCoeff >= 0)
         {
@@ -739,10 +895,8 @@ void Foam::BlockMatrixClustering<Type>::restrictDiagDecoupled
         linearTypeField& activeCoarseCoeff = coarseCoeff.asLinear();
         const linearTypeField& activeCoeff = Coeff.asLinear();
 
-        forAll (coarseCoeff, i)
+        // Reset coefficients to zero
-        {
+        activeCoarseCoeff = pTraits<linearType>::zero;
-            activeCoarseCoeff[i] = pTraits<linearType>::zero;
-        }
 
         forAll (Coeff, i)
         {
@@ -761,10 +915,8 @@ void Foam::BlockMatrixClustering<Type>::restrictDiagDecoupled
         scalarTypeField& activeCoarseCoeff = coarseCoeff.asScalar();
         const scalarTypeField& activeCoeff = Coeff.asScalar();
 
-        forAll (coarseCoeff, i)
+        // Reset coefficients to zero
-        {
+        activeCoarseCoeff = pTraits<scalarType>::zero;
-            activeCoarseCoeff[i] = pTraits<scalarType>::zero;
-        }
 
         forAll (Coeff, i)
         {
@@ -800,6 +952,7 @@ Foam::BlockMatrixClustering<Type>::BlockMatrixClustering
     maxGroupSize_(readLabel(dict.lookup("maxGroupSize"))),
     normPtr_(BlockCoeffNorm<Type>::New(dict)),
     agglomIndex_(matrix_.lduAddr().size()),
+    coeffRestrictAddr_(),
     groupSize_(groupSize),
     nSolo_(0),
     nCoarseEqns_(0),
@@ -850,8 +1003,6 @@ Foam::BlockMatrixClustering<Type>::restrictMatrix() const
     const unallocLabelList& upperAddr = matrix_.lduAddr().upperAddr();
     const unallocLabelList& lowerAddr = matrix_.lduAddr().lowerAddr();
 
-    const label nFineCoeffs = upperAddr.size();
-
 #   ifdef FULLDEBUG
     if (agglomIndex_.size() != matrix_.lduAddr().size())
     {
@@ -866,161 +1017,6 @@ Foam::BlockMatrixClustering<Type>::restrictMatrix() const
     }
 #   endif
 
-    // Does the matrix have solo equations
-    bool soloEqns = nSolo_ > 0;
-
-    // Storage for block neighbours and coefficients
-
-    // Guess initial maximum number of neighbours in block
-    label maxNnbrs = 10;
-
-    // Number of neighbours per block
-    labelList blockNnbrs(nCoarseEqns_, 0);
-
-    // Setup initial packed storage for neighbours and coefficients
-    labelList blockNbrsData(maxNnbrs*nCoarseEqns_);
-
-    // Create face-restriction addressing
-    labelList coeffRestrictAddr(nFineCoeffs);
-
-    // Initial neighbour array (not in upper-triangle order)
-    labelList initCoarseNeighb(nFineCoeffs);
-
-    // Counter for coarse coeffs
-    label nCoarseCoeffs = 0;
-
-    // Loop through all fine coeffs
-    forAll (upperAddr, fineCoeffi)
-    {
-        label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffi]];
-        label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffi]];
-
-        // If the coefficient touches block zero and solo equations are
-        // present, skip it
-        if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0))
-        {
-            continue;
-        }
-
-        if (rmUpperAddr == rmLowerAddr)
-        {
-            // For each fine coeff inside of a coarse cluster keep the address
-            // of the cluster corresponding to the coeff in the
-            // coeffRestrictAddr as a negative index
-            coeffRestrictAddr[fineCoeffi] = -(rmUpperAddr + 1);
-        }
-        else
-        {
-            // This coeff is a part of a coarse coeff
-
-            label cOwn = rmUpperAddr;
-            label cNei = rmLowerAddr;
-
-            // Get coarse owner and neighbour
-            if (rmUpperAddr > rmLowerAddr)
-            {
-                cOwn = rmLowerAddr;
-                cNei = rmUpperAddr;
-            }
-
-            // Check the neighbour to see if this coeff has already been found
-            bool nbrFound = false;
-            label& ccnCoeffs = blockNnbrs[cOwn];
-
-            for (int i = 0; i < ccnCoeffs; i++)
-            {
-                if (initCoarseNeighb[blockNbrsData[maxNnbrs*cOwn + i]] == cNei)
-                {
-                    nbrFound = true;
-                    coeffRestrictAddr[fineCoeffi] =
-                        blockNbrsData[maxNnbrs*cOwn + i];
-                    break;
-                }
-            }
-
-            if (!nbrFound)
-            {
-                if (ccnCoeffs >= maxNnbrs)
-                {
-                    // Double the size of list and copy data
-                    label oldMaxNnbrs = maxNnbrs;
-                    maxNnbrs *= 2;
-
-                    // Resize and copy list
-                    const labelList oldBlockNbrsData = blockNbrsData;
-                    blockNbrsData.setSize(maxNnbrs*nCoarseEqns_);
-
-                    forAll (blockNnbrs, i)
-                    {
-                        for (int j = 0; j < blockNnbrs[i]; j++)
-                        {
-                            blockNbrsData[maxNnbrs*i + j] =
-                                oldBlockNbrsData[oldMaxNnbrs*i + j];
-                        }
-                    }
-                }
-
-                blockNbrsData[maxNnbrs*cOwn + ccnCoeffs] = nCoarseCoeffs;
-                initCoarseNeighb[nCoarseCoeffs] = cNei;
-                coeffRestrictAddr[fineCoeffi] = nCoarseCoeffs;
-                ccnCoeffs++;
-
-                // New coarse coeff created
-                nCoarseCoeffs++;
-            }
-        }
-    } // End for all fine coeffs
-
-
-    // Renumber into upper-triangular order
-
-    // All coarse owner-neighbour storage
-    labelList coarseOwner(nCoarseCoeffs);
-    labelList coarseNeighbour(nCoarseCoeffs);
-    labelList coarseCoeffMap(nCoarseCoeffs);
-
-    label coarseCoeffi = 0;
-
-    forAll (blockNnbrs, cci)
-    {
-        label* cCoeffs = &blockNbrsData[maxNnbrs*cci];
-        label ccnCoeffs = blockNnbrs[cci];
-
-        for (int i = 0; i < ccnCoeffs; i++)
-        {
-            coarseOwner[coarseCoeffi] = cci;
-            coarseNeighbour[coarseCoeffi] = initCoarseNeighb[cCoeffs[i]];
-            coarseCoeffMap[cCoeffs[i]] = coarseCoeffi;
-            coarseCoeffi++;
-        }
-    }
-
-    forAll(coeffRestrictAddr, fineCoeffi)
-    {
-        label rmUpperAddr = agglomIndex_[upperAddr[fineCoeffi]];
-        label rmLowerAddr = agglomIndex_[lowerAddr[fineCoeffi]];
-
-        // If the coefficient touches block zero and solo equations are
-        // present, skip it
-        if (soloEqns && (rmUpperAddr == 0 || rmLowerAddr == 0))
-        {
-            continue;
-        }
-
-        if (coeffRestrictAddr[fineCoeffi] >= 0)
-        {
-            coeffRestrictAddr[fineCoeffi] =
-                coarseCoeffMap[coeffRestrictAddr[fineCoeffi]];
-        }
-    }
-
-    // Clear the temporary storage for the coarse matrix data
-    blockNnbrs.setSize(0);
-    blockNbrsData.setSize(0);
-    initCoarseNeighb.setSize(0);
-    coarseCoeffMap.setSize(0);
-
-
     // Create coarse-level coupled interfaces
 
     // Create coarse interfaces, addressing and coefficients
@@ -1028,8 +1024,7 @@ Foam::BlockMatrixClustering<Type>::restrictMatrix() const
         const_cast<BlockLduMatrix<Type>& >(matrix_).interfaces().size();
 
     const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
-        interfaceFields =
+        interfaceFields = matrix_.interfaces();
-        const_cast<BlockLduMatrix<Type>&>(matrix_).interfaces();
 
     // Set the coarse interfaces and coefficients
     lduInterfacePtrsList coarseInterfaces(interfaceSize);
@@ -1046,7 +1041,7 @@ Foam::BlockMatrixClustering<Type>::restrictMatrix() const
             nCoarseEqns_,
             coarseOwner,
             coarseNeighbour,
-            Pstream::worldComm, //HJ, AMG Comm fineMesh.comm(),
+            Pstream::worldComm,
             true
         )
     );
@@ -1140,7 +1135,39 @@ Foam::BlockMatrixClustering<Type>::restrictMatrix() const
     );
     BlockLduMatrix<Type>& coarseMatrix = coarseMatrixPtr();
 
+    // Build matrix coefficients
+    this->updateMatrix(coarseMatrix);
+
+    // Create and return BlockAMGLevel
+    return autoPtr<BlockAMGLevel<Type> >
+    (
+        new coarseBlockAMGLevel<Type>
+        (
+            coarseAddrPtr,
+            coarseMatrixPtr,
+            this->dict(),
+            this->type(),
+            this->groupSize(),
+            this->minCoarseEqns()
+        )
+    );
+}
+
+
+template<class Type>
+void Foam::BlockMatrixClustering<Type>::updateMatrix
+(
+    BlockLduMatrix<Type>& coarseMatrix
+) const
+{
+    // Get interfaces from fine matrix
+    const typename BlockLduInterfaceFieldPtrsList<Type>::Type&
+        interfaceFields = matrix_.interfaces();
+
     // Get interfaces from coarse matrix
+    lduInterfacePtrsList coarseInterfaces = coarseMatrix.mesh().interfaces();
+    
+    // Get interfaces fields from coarse matrix
     typename BlockLduInterfaceFieldPtrsList<Type>::Type&
         coarseInterfaceFieldsTransfer = coarseMatrix.interfaces();
 
@@ -1242,7 +1269,6 @@ Foam::BlockMatrixClustering<Type>::restrictMatrix() const
 
             agglomerateCoeffs
             (
-                coeffRestrictAddr,
                 activeCoarseDiag,
                 activeCoarseUpper,
                 activeFineUpper,
@@ -1264,7 +1290,6 @@ Foam::BlockMatrixClustering<Type>::restrictMatrix() const
 
             agglomerateCoeffs
             (
-                coeffRestrictAddr,
                 activeCoarseDiag,
                 activeCoarseUpper,
                 activeFineUpper,
@@ -1286,7 +1311,6 @@ Foam::BlockMatrixClustering<Type>::restrictMatrix() const
 
             agglomerateCoeffs
             (
-                coeffRestrictAddr,
                 activeCoarseDiag,
                 activeCoarseUpper,
                 activeFineUpper,
@@ -1326,7 +1350,6 @@ Foam::BlockMatrixClustering<Type>::restrictMatrix() const
 
             agglomerateCoeffs
             (
-                coeffRestrictAddr,
                 activeCoarseDiag,
                 activeCoarseUpper,
                 activeFineUpper,
@@ -1352,7 +1375,6 @@ Foam::BlockMatrixClustering<Type>::restrictMatrix() const
 
             agglomerateCoeffs
             (
-                coeffRestrictAddr,
                 activeCoarseDiag,
                 activeCoarseUpper,
                 activeFineUpper,
@@ -1378,7 +1400,6 @@ Foam::BlockMatrixClustering<Type>::restrictMatrix() const
 
             agglomerateCoeffs
             (
-                coeffRestrictAddr,
                 activeCoarseDiag,
                 activeCoarseUpper,
                 activeFineUpper,
@@ -1396,20 +1417,6 @@ Foam::BlockMatrixClustering<Type>::restrictMatrix() const
                 << abort(FatalError);
         }
     }
-
-    // Create and return BlockAMGLevel
-    return autoPtr<BlockAMGLevel<Type> >
-    (
-        new coarseBlockAMGLevel<Type>
-        (
-            coarseAddrPtr,
-            coarseMatrixPtr,
-            this->dict(),
-            this->type(),
-            this->groupSize(),
-            this->minCoarseEqns()
-        )
-    );
 }
 
 

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

@@ -76,7 +76,10 @@ class BlockMatrixClustering
         autoPtr<BlockCoeffNorm<Type> > normPtr_;
 
         //- Child array: for each fine equation give a clustering index
-        labelField agglomIndex_;
+        labelList agglomIndex_;
+
+        //- Face-restriction addressing
+        labelList coeffRestrictAddr_;
 
         //- Group size
         label groupSize_;
@@ -116,7 +119,6 @@ class BlockMatrixClustering
         template<class DiagType, class ULType>
         void agglomerateCoeffs
         (
-            const labelList& coeffRestrictAddr,
             Field<DiagType>& activeCoarseDiag,
             Field<ULType>& activeCoarseUpper,
             const Field<ULType>& activeFineUpper,
@@ -127,7 +129,6 @@ class BlockMatrixClustering
         template<class DiagType, class ULType>
         void agglomerateCoeffs
         (
-            const labelList& coeffRestrictAddr,
             Field<DiagType>& activeCoarseDiag,
             Field<ULType>& activeCoarseUpper,
             const Field<ULType>& activeFineUpper,
@@ -197,6 +198,9 @@ public:
             Field<Type>& x,
             const Field<Type>& coarseX
         ) const;
+
+        //- Update coarse matrix using same coefficients
+        virtual void updateMatrix(BlockLduMatrix<Type>& coarseMatrix) const;
 };
 
 

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

@@ -172,8 +172,7 @@ public:
         virtual bool coarsen() const = 0;
 
         //- Restrict matrix
-        virtual autoPtr<BlockAMGLevel<Type> > restrictMatrix
+        virtual autoPtr<BlockAMGLevel<Type> > restrictMatrix() const = 0;
-        () const = 0;
 
         //- Restrict residual
         virtual void restrictResidual
@@ -189,6 +188,8 @@ public:
             const Field<Type>& coarseX
         ) const = 0;
 
+        //- Update coarse matrix using same coefficients
+        virtual void updateMatrix(BlockLduMatrix<Type>& coarseMatrix) const = 0;
 };
 
 
@@ -203,15 +204,15 @@ public:
 #endif
 
 #define makeBlockMatrixCoarsening(BlockMatrixCoarseningType, typeBlockMatrixCoarseningType) \
-                                                                             \
+                                                                               \
-defineNamedTemplateTypeNameAndDebug(typeBlockMatrixCoarseningType, 0);              \
+defineNamedTemplateTypeNameAndDebug(typeBlockMatrixCoarseningType, 0);         \
-                                                                             \
+                                                                               \
 addToRunTimeSelectionTable(BlockMatrixCoarseningType, typeBlockMatrixCoarseningType, matrix);
 
-#define makeBlockMatrixCoarsenings(blockAmgCoarseningType)                               \
+#define makeBlockMatrixCoarsenings(blockAmgCoarseningType)                     \
-                                                                             \
+                                                                               \
-makeBlockMatrixCoarsening(blockScalarMatrixCoarsening, blockAmgCoarseningType##Scalar);         \
+makeBlockMatrixCoarsening(blockScalarMatrixCoarsening, blockAmgCoarseningType##Scalar); \
-makeBlockMatrixCoarsening(blockVectorMatrixCoarsening, blockAmgCoarseningType##Vector);         \
+makeBlockMatrixCoarsening(blockVectorMatrixCoarsening, blockAmgCoarseningType##Vector); \
 makeBlockMatrixCoarsening(blockTensorMatrixCoarsening, blockAmgCoarseningType##Tensor);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -989,8 +989,7 @@ Foam::BlockMatrixSelection<Type>::restrictMatrix() const
             nCoarseEqns_,
             coarseOwner,
             coarseNeighbour,
-            // true
+            false           // DO NOT REUSE STORAGE!
-            false           // DO NOT REUSE STORAGE!  HJ, HERE
         )
     );
 

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

@@ -104,6 +104,13 @@ Foam::coarseBlockAMGLevel<Type>::~coarseBlockAMGLevel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+Foam::BlockLduMatrix<Type>& Foam::coarseBlockAMGLevel<Type>::matrix()
+{
+    return matrixPtr_();
+}
+
+
 template<class Type>
 Foam::Field<Type>& Foam::coarseBlockAMGLevel<Type>::x()
 {
@@ -334,4 +341,23 @@ Foam::coarseBlockAMGLevel<Type>::makeNextLevel() const
 }
 
 
+template<class Type>
+void Foam::coarseBlockAMGLevel<Type>::initLevel
+(
+    autoPtr<Foam::BlockAMGLevel<Type> >& coarseLevelPtr
+)
+{
+    // Update matrix by repeating coarsening
+    
+    // Update smoother for new matrix
+    smootherPtr_->initMatrix();
+
+    // Update coarse matrix if it exists
+    if (coarseLevelPtr.valid())
+    {
+        coarseningPtr_->updateMatrix(coarseLevelPtr->matrix());
+    }
+}
+
+
 // ************************************************************************* //

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

@@ -125,6 +125,9 @@ public:
             return dict_;
         }
 
+        //- Return reference to matrix
+        virtual BlockLduMatrix<Type>& matrix();
+
         //- Return reference to x
         virtual Field<Type>& x();
 
@@ -184,6 +187,13 @@ public:
 
         //- Create next level from current level
         virtual autoPtr<BlockAMGLevel<Type> > makeNextLevel() const;
+
+        //- Re-initialise AMG level after matrix coefficient update
+        //  and update matrix coefficients on coarse level using same coarsening
+        virtual void initLevel
+        (
+            autoPtr<Foam::BlockAMGLevel<Type> >& coarseLevelPtr
+        );
 };
 
 

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

@@ -81,11 +81,23 @@ Foam::fineBlockAMGLevel<Type>::fineBlockAMGLevel
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+Foam::BlockLduMatrix<Type>& Foam::fineBlockAMGLevel<Type>::matrix()
+{
+    FatalErrorIn("Field<Type>& Foam::fineBlockAMGLevel<Type>::matrix()")
+        << "Access to matrix is not available on fine level."
+        << abort(FatalError);
+
+    // Dummy return
+    return const_cast<BlockLduMatrix<Type>&>(matrix_);
+}
+
+
 template<class Type>
 Foam::Field<Type>& Foam::fineBlockAMGLevel<Type>::x()
 {
     FatalErrorIn("Field<Type>& Foam::fineBlockAMGLevel<Type>::x()")
-        << "x is not available."
+        << "x is not available on fine level."
         << abort(FatalError);
 
     // Dummy return
@@ -97,7 +109,7 @@ template<class Type>
 Foam::Field<Type>& Foam::fineBlockAMGLevel<Type>::b()
 {
     FatalErrorIn("Field<Type>& Foam::fineBlockAMGLevel<Type>::b()")
-        << "b is not available."
+        << "b is not available on fine level."
         << abort(FatalError);
 
     // Dummy return
@@ -279,4 +291,23 @@ Foam::fineBlockAMGLevel<Type>::makeNextLevel() const
 }
 
 
+template<class Type>
+void Foam::fineBlockAMGLevel<Type>::initLevel
+(
+    autoPtr<Foam::BlockAMGLevel<Type> >& coarseLevelPtr
+)
+{
+    // Fine matrix has been updated externally
+
+    // Update smoother for new matrix
+    smootherPtr_->initMatrix();
+
+    // Update coarse matrix if it exists
+    if (coarseLevelPtr.valid())
+    {
+        coarseningPtr_->updateMatrix(coarseLevelPtr->matrix());
+    }
+}
+
+
 // ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockAMG/fineBlockAMGLevel.H

@@ -125,6 +125,9 @@ public:
             return dict_;
         }
 
+        //- Return reference to matrix
+        virtual BlockLduMatrix<Type>& matrix();
+
         //- Return reference to x
         virtual Field<Type>& x();
 
@@ -183,6 +186,13 @@ public:
 
         //- Create next level from current level
         virtual autoPtr<BlockAMGLevel<Type> > makeNextLevel() const;
+
+        //- Re-initialise AMG level after matrix coefficient update
+        //  and update matrix coefficients on coarse level using same coarsening
+        virtual void initLevel
+        (
+            autoPtr<Foam::BlockAMGLevel<Type> >& coarseLevelPtr
+        );
 };
 
 

src/foam/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrix.C

@@ -67,7 +67,6 @@ Foam::BlockLduMatrix<Type>::BlockLduMatrix(const lduMesh& ldu)
         coupleUpper_.set(i, new CoeffField<Type>(addr.patchAddr(i).size()));
         coupleLower_.set(i, new CoeffField<Type>(addr.patchAddr(i).size()));
     }
-
 }
 
 

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockAMGPrecon/BlockAMGPrecon.C

@@ -141,4 +141,11 @@ void Foam::BlockAMGPrecon<Type>::precondition
 }
 
 
+template<class Type>
+void Foam::BlockAMGPrecon<Type>::initMatrix()
+{
+    amgPtr_->initMatrix();
+}
+
+
 // ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockAMGPrecon/BlockAMGPrecon.H

@@ -25,7 +25,7 @@ Class
     BlockAMGPrecon
 
 Description
-    AMG preconditioning
+    Algebraic Multigrid preconditioning for block matrices
 
 Author
     Klas Jareteg, 2012-12-13
@@ -55,7 +55,7 @@ template<class Type>
 class BlockAMGCycle;
 
 /*---------------------------------------------------------------------------*\
-                       Class BlockAMGPrecon Declaration
+                         Class BlockAMGPrecon Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Type>
@@ -141,6 +141,9 @@ public:
             Field<Type>& x,
             const Field<Type>& b
         ) const;
+
+        //- Re-initialise preconditioner after matrix coefficient update
+        virtual void initMatrix();
 };
 
 

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.C

@@ -914,7 +914,7 @@ Foam::BlockCholeskyPrecon<Type>::BlockCholeskyPrecon
     BlockLduPrecon<Type>(matrix),
     preconDiag_(matrix.diag())
 {
-    calcPreconDiag();
+    this->calcPreconDiag();
 }
 
 
@@ -1301,4 +1301,12 @@ void Foam::BlockCholeskyPrecon<Type>::preconditionT
 }
 
 
+template<class Type>
+void Foam::BlockCholeskyPrecon<Type>::initMatrix()
+{
+    preconDiag_.clear();
+    this->calcPreconDiag();
+}
+
+
 // ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockCholeskyPrecon/BlockCholeskyPrecon.H

@@ -25,7 +25,7 @@ Class
     BlockCholeskyPrecon
 
 Description
-    Incomplete Cholesky preconditioning with no fill-in.
+    Incomplete Cholesky preconditioning with no fill-in for block matrices.
 
 Author
     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
@@ -215,6 +215,9 @@ public:
             Field<Type>& xT,
             const Field<Type>& bT
         ) const;
+
+        //- Re-initialise preconditioner after matrix coefficient update
+        virtual void initMatrix();
 };
 
 

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockDiagonalPrecon/BlockDiagonalPrecon.H

@@ -107,6 +107,10 @@ public:
         {
             return precondition(xT, bT);
         }
+
+        //- Re-initialise preconditioner after matrix coefficient update
+        virtual void initMatrix()
+        {}
 };
 
 

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.C

@@ -946,4 +946,14 @@ void Foam::BlockGaussSeidelPrecon<Type>::preconditionT
 }
 
 
+template<class Type>
+void Foam::BlockGaussSeidelPrecon<Type>::initMatrix()
+{
+    invDiag_.clear();
+    LUDiag_.clear();
+
+    this->calcInvDiag();
+}
+
+
 // ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockGaussSeidelPrecon/BlockGaussSeidelPrecon.H

@@ -172,6 +172,9 @@ public:
             Field<Type>& xT,
             const Field<Type>& bT
         ) const;
+
+        //- Re-initialise preconditioner after matrix coefficient update
+        virtual void initMatrix();
 };
 
 

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.C

@@ -840,4 +840,17 @@ void Foam::BlockILUCpPrecon<Type>::preconditionT
 }
 
 
+template<class Type>
+void Foam::BlockILUCpPrecon<Type>::initMatrix()
+{
+    // Clear extended matrix
+    extBlockMatrix_.extendedLower().clear();
+    extBlockMatrix_.extendedLower().clear();
+
+    preconDiag_.clear();
+
+    this->calcFactorization();
+}
+
+
 // ************************************************************************* //

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockILUCpPrecon/BlockILUCpPrecon.H

@@ -168,6 +168,9 @@ public:
             Field<Type>& xT,
             const Field<Type>& bT
         ) const;
+
+        //- Re-initialise preconditioner after matrix coefficient update
+        virtual void initMatrix();
 };
 
 

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockLduPrecon/BlockLduPrecon.H

@@ -140,6 +140,9 @@ public:
                 "(Field<Type>& xT, const Field<Type>& bT) const"
             );
         }
+
+        //- Re-initialise preconditioner after matrix coefficient update
+        virtual void initMatrix() = 0;
 };
 
 

src/foam/matrices/blockLduMatrix/BlockLduPrecons/BlockNoPrecon/BlockNoPrecon.H

@@ -109,6 +109,10 @@ public:
         {
             precondition(xT, bT);
         }
+
+        //- Re-initialise preconditioner after matrix coefficient update
+        virtual void initMatrix()
+        {}
 };
 
 

src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockGaussSeidelSmoother/BlockGaussSeidelSmoother.H

@@ -111,6 +111,12 @@ public:
                 gs_.precondition(x, b);
             }
         }
+
+        //- Re-initialise preconditioner after matrix coefficient update
+        virtual void initMatrix()
+        {
+            gs_.initMatrix();
+        }
 };
 
 

src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUCpSmoother/BlockILUCpSmoother.H

@@ -127,6 +127,12 @@ public:
                 x += xCorr_;
             }
         }
+
+        //- Re-initialise preconditioner after matrix coefficient update
+        virtual void initMatrix()
+        {
+            precon_.initMatrix();
+        }
 };
 
 

src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockILUSmoother/BlockILUSmoother.H

@@ -133,6 +133,12 @@ public:
                 x += xCorr_;
             }
         }
+
+        //- Re-initialise preconditioner after matrix coefficient update
+        virtual void initMatrix()
+        {
+            precon_.initMatrix();
+        }
 };
 
 

src/foam/matrices/blockLduMatrix/BlockLduSmoothers/BlockLduSmoother/BlockLduSmoother.H

@@ -125,6 +125,9 @@ public:
             const Field<Type>& b,
             const label nSweeps
         ) const = 0;
+
+        //- Re-initialise preconditioner after matrix coefficient update
+        virtual void initMatrix() = 0;
 };
 
 

src/foam/matrices/blockLduMatrix/BlockLduSolvers/BlockAMGSolver/BlockAMGSolver.H

@@ -51,7 +51,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class BlockAMGSolver Declaration
+                      Class BlockAMGSolver Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Type>
@@ -93,10 +93,9 @@ public:
         );
 
 
-    // Destructor
+    //- Destructor
-
+    virtual ~BlockAMGSolver()
-        virtual ~BlockAMGSolver()
+    {}
-        {}
 
 
     // Member Functions