ENH: generalise masterCoarsest restart (GAMG)

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,11 +48,15 @@ namespace Foam
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
+void Foam::GAMGAgglomeration::compactLevels
+(
+    const label nCreatedLevels,
+    const bool doProcessorAgglomerate
+)
 {
-    nCells_.setSize(nCreatedLevels);
+    nCells_.setSize(nCreatedLevels, 0);
     restrictAddressing_.setSize(nCreatedLevels);
-    nFaces_.setSize(nCreatedLevels);
+    nFaces_.setSize(nCreatedLevels, 0);
     faceRestrictAddressing_.setSize(nCreatedLevels);
     faceFlipMap_.setSize(nCreatedLevels);
     nPatchFaces_.setSize(nCreatedLevels);
@@ -60,8 +64,8 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
     meshLevels_.setSize(nCreatedLevels);
 
     // Have procCommunicator_ always, even if not procAgglomerating
-    procCommunicator_.setSize(nCreatedLevels + 1);
+    procCommunicator_.setSize(nCreatedLevels + 1, -1);
-    if (processorAgglomerate())
+    if (doProcessorAgglomerate && processorAgglomerate())
     {
         procAgglomMap_.setSize(nCreatedLevels);
         agglomProcIDs_.setSize(nCreatedLevels);
@@ -71,156 +75,160 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
         procBoundaryFaceMap_.setSize(nCreatedLevels);
 
         procAgglomeratorPtr_().agglomerate();
+    }
+}
 
 
+void Foam::GAMGAgglomeration::printLevels() const
+{
+    Info<< "GAMGAgglomeration:" << nl
+        << "    local agglomerator     : " << type() << nl;
+    if (processorAgglomerate())
+    {
+        Info<< "    processor agglomerator : "
+            << procAgglomeratorPtr_().type() << nl
+            << nl;
     }
 
-    // Print a bit
+    Info<< setw(36) << "nCells"
-    if (debug)
+        << setw(20) << "nFaces/nCells"
+        << setw(20) << "nInterfaces"
+        << setw(20) << "nIntFaces/nCells"
+        << setw(12) << "profile"
+        << nl
+        << setw(8) << "Level"
+        << setw(8) << "nProcs"
+        << "    "
+        << setw(8) << "avg"
+        << setw(8) << "max"
+        << "    "
+        << setw(8) << "avg"
+        << setw(8) << "max"
+        << "    "
+        << setw(8) << "avg"
+        << setw(8) << "max"
+        << "    "
+        << setw(8) << "avg"
+        << setw(8) << "max"
+        //<< "    "
+        << setw(12) << "avg"
+        << nl
+        << setw(8) << "-----"
+        << setw(8) << "------"
+        << "    "
+        << setw(8) << "---"
+        << setw(8) << "---"
+        << "    "
+        << setw(8) << "---"
+        << setw(8) << "---"
+        << "    "
+        << setw(8) << "---"
+        << setw(8) << "---"
+        << "    "
+        << setw(8) << "---"
+        << setw(8) << "---"
+        //<< "    "
+        << setw(12) << "---"
+        //<< "    "
+        << nl;
+
+    const label maxSize = returnReduce(size(), maxOp<label>());
+
+    for (label levelI = 0; levelI <= maxSize; levelI++)
     {
-        Info<< "GAMGAgglomeration:" << nl
+        label nProcs = 0;
-            << "    local agglomerator     : " << type() << nl;
+        label nCells = 0;
-        if (processorAgglomerate())
+        scalar faceCellRatio = 0;
+        label nInterfaces = 0;
+        label nIntFaces = 0;
+        scalar ratio = 0.0;
+        scalar profile = 0.0;
+
+        if (hasMeshLevel(levelI))
         {
-            Info<< "    processor agglomerator : "
+            nProcs = 1;
-                << procAgglomeratorPtr_().type() << nl
+
-                << nl;
+            const lduMesh& fineMesh = meshLevel(levelI);
+            nCells = fineMesh.lduAddr().size();
+            faceCellRatio =
+                scalar(fineMesh.lduAddr().lowerAddr().size())/nCells;
+
+            const lduInterfacePtrsList interfaces =
+                fineMesh.interfaces();
+            forAll(interfaces, i)
+            {
+                if (interfaces.set(i))
+                {
+                    nInterfaces++;
+                    nIntFaces += interfaces[i].faceCells().size();
+                }
+            }
+            ratio = scalar(nIntFaces)/nCells;
+
+            profile = fineMesh.lduAddr().band().second();
         }
 
-        Info<< setw(36) << "nCells"
+        label totNprocs = returnReduce(nProcs, sumOp<label>());
-            << setw(20) << "nFaces/nCells"
+
-            << setw(20) << "nInterfaces"
+        label maxNCells = returnReduce(nCells, maxOp<label>());
-            << setw(20) << "nIntFaces/nCells"
+        label totNCells = returnReduce(nCells, sumOp<label>());
-            << setw(12) << "profile"
+
-            << nl
+        scalar maxFaceCellRatio =
-            << setw(8) << "Level"
+            returnReduce(faceCellRatio, maxOp<scalar>());
-            << setw(8) << "nProcs"
+        scalar totFaceCellRatio =
+            returnReduce(faceCellRatio, sumOp<scalar>());
+
+        label maxNInt = returnReduce(nInterfaces, maxOp<label>());
+        label totNInt = returnReduce(nInterfaces, sumOp<label>());
+
+        scalar maxRatio = returnReduce(ratio, maxOp<scalar>());
+        scalar totRatio = returnReduce(ratio, sumOp<scalar>());
+
+        scalar totProfile = returnReduce(profile, sumOp<scalar>());
+
+        const int oldPrecision = Info.stream().precision(4);
+
+        Info<< setw(8) << levelI
+            << setw(8) << totNprocs
             << "    "
-            << setw(8) << "avg"
+            << setw(8) << totNCells/totNprocs
-            << setw(8) << "max"
+            << setw(8) << maxNCells
             << "    "
-            << setw(8) << "avg"
+            << setw(8) << totFaceCellRatio/totNprocs
-            << setw(8) << "max"
+            << setw(8) << maxFaceCellRatio
             << "    "
-            << setw(8) << "avg"
+            << setw(8) << scalar(totNInt)/totNprocs
-            << setw(8) << "max"
+            << setw(8) << maxNInt
             << "    "
-            << setw(8) << "avg"
+            << setw(8) << totRatio/totNprocs
-            << setw(8) << "max"
+            << setw(8) << maxRatio
-            //<< "    "
+            << setw(12) << totProfile/totNprocs
-            << setw(12) << "avg"
-            << nl
-            << setw(8) << "-----"
-            << setw(8) << "------"
-            << "    "
-            << setw(8) << "---"
-            << setw(8) << "---"
-            << "    "
-            << setw(8) << "---"
-            << setw(8) << "---"
-            << "    "
-            << setw(8) << "---"
-            << setw(8) << "---"
-            << "    "
-            << setw(8) << "---"
-            << setw(8) << "---"
-            //<< "    "
-            << setw(12) << "---"
-            //<< "    "
             << nl;
 
-        for (label levelI = 0; levelI <= size(); levelI++)
+        Info.stream().precision(oldPrecision);
-        {
-            label nProcs = 0;
-            label nCells = 0;
-            scalar faceCellRatio = 0;
-            label nInterfaces = 0;
-            label nIntFaces = 0;
-            scalar ratio = 0.0;
-            scalar profile = 0.0;
-
-            if (hasMeshLevel(levelI))
-            {
-                nProcs = 1;
-
-                const lduMesh& fineMesh = meshLevel(levelI);
-                nCells = fineMesh.lduAddr().size();
-                faceCellRatio =
-                    scalar(fineMesh.lduAddr().lowerAddr().size())/nCells;
-
-                const lduInterfacePtrsList interfaces =
-                    fineMesh.interfaces();
-                forAll(interfaces, i)
-                {
-                    if (interfaces.set(i))
-                    {
-                        nInterfaces++;
-                        nIntFaces += interfaces[i].faceCells().size();
-                    }
-                }
-                ratio = scalar(nIntFaces)/nCells;
-
-                profile = fineMesh.lduAddr().band().second();
-            }
-
-            label totNprocs = returnReduce(nProcs, sumOp<label>());
-
-            label maxNCells = returnReduce(nCells, maxOp<label>());
-            label totNCells = returnReduce(nCells, sumOp<label>());
-
-            scalar maxFaceCellRatio =
-                returnReduce(faceCellRatio, maxOp<scalar>());
-            scalar totFaceCellRatio =
-                returnReduce(faceCellRatio, sumOp<scalar>());
-
-            label maxNInt = returnReduce(nInterfaces, maxOp<label>());
-            label totNInt = returnReduce(nInterfaces, sumOp<label>());
-
-            scalar maxRatio = returnReduce(ratio, maxOp<scalar>());
-            scalar totRatio = returnReduce(ratio, sumOp<scalar>());
-
-            scalar totProfile = returnReduce(profile, sumOp<scalar>());
-
-            const int oldPrecision = Info.stream().precision(4);
-
-            Info<< setw(8) << levelI
-                << setw(8) << totNprocs
-                << "    "
-                << setw(8) << totNCells/totNprocs
-                << setw(8) << maxNCells
-                << "    "
-                << setw(8) << totFaceCellRatio/totNprocs
-                << setw(8) << maxFaceCellRatio
-                << "    "
-                << setw(8) << scalar(totNInt)/totNprocs
-                << setw(8) << maxNInt
-                << "    "
-                << setw(8) << totRatio/totNprocs
-                << setw(8) << maxRatio
-                << setw(12) << totProfile/totNprocs
-                << nl;
-
-            Info.stream().precision(oldPrecision);
-        }
-        Info<< endl;
     }
+    Info<< endl;
 }
 
 
 bool Foam::GAMGAgglomeration::continueAgglomerating
 (
+    const label nCellsInCoarsestLevel,
     const label nFineCells,
-    const label nCoarseCells
+    const label nCoarseCells,
+    const label comm
 ) const
 {
-    const label nTotalCoarseCells = returnReduce(nCoarseCells, sumOp<label>());
+    const label nTotalCoarseCells =
-    if (nTotalCoarseCells < Pstream::nProcs()*nCellsInCoarsestLevel_)
+        returnReduce(nCoarseCells, sumOp<label>(), UPstream::msgType(), comm);
+    if (nTotalCoarseCells < Pstream::nProcs(comm)*nCellsInCoarsestLevel)
     {
         return false;
     }
     else
     {
-        const label nTotalFineCells = returnReduce(nFineCells, sumOp<label>());
+        const label nTotalFineCells =
+            returnReduce(nFineCells, sumOp<label>(), UPstream::msgType(), comm);
         return nTotalCoarseCells < nTotalFineCells;
     }
 }
@@ -333,7 +341,12 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
                 << exit(FatalError);
         }
 
-        return store(ctorPtr(mesh, controlDict).ptr());
+        auto agglomPtr(ctorPtr(mesh, controlDict));
+        if (debug)
+        {
+            agglomPtr().printLevels();
+        }
+        return store(agglomPtr.ptr());
     }
 }
 
@@ -378,7 +391,12 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
         }
         else
         {
-            return store(ctorPtr(matrix, controlDict).ptr());
+            auto agglomPtr(ctorPtr(matrix, controlDict));
+            if (debug)
+            {
+                agglomPtr().printLevels();
+            }
+            return store(agglomPtr.ptr());
         }
     }
 }
@@ -429,7 +447,7 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
                 << exit(FatalError);
         }
 
-        return store
+        auto agglomPtr
         (
             ctorPtr
             (
@@ -437,8 +455,13 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
                 cellVolumes,
                 faceAreas,
                 controlDict
-            ).ptr()
+            )
         );
+        if (debug)
+        {
+            agglomPtr().printLevels();
+        }
+        return store(agglomPtr.ptr());
     }
 }
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -53,6 +53,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class lduMesh;
 class lduMatrix;
 class mapDistribute;
@@ -151,18 +152,28 @@ protected:
         //- Combine a level with the previous one
         void combineLevels(const label curLevel);
 
-        //- Shrink the number of levels to that specified
+        //- Shrink the number of levels to that specified. Optionally do
-        void compactLevels(const label nCreatedLevels);
+        //  processor agglomeration
+        void compactLevels
+        (
+            const label nCreatedLevels,
+            const bool doProcessorAgglomerate
+        );
 
         //- Check the need for further agglomeration
         bool continueAgglomerating
         (
+            const label nCellsInCoarsestLevel,
             const label nCells,
-            const label nCoarseCells
+            const label nCoarseCells,
+            const label comm
         ) const;
 
         void clearLevel(const label leveli);
 
+        //- Print level overview
+        void printLevels() const;
+
 
         // Processor agglomeration
 
@@ -468,16 +479,30 @@ public:
             const labelListListList& boundaryFaceMap(const label fineLeveli)
             const;
 
-        //- Given restriction determines if coarse cells are connected.
+
-        //  Return ok is so, otherwise creates new restriction that is
+        // Helpers
-        static bool checkRestriction
+
-        (
+            //- Agglomerate from a starting level. Starting level is usually 0
-            labelList& newRestrict,
+            //- (initial mesh) but sometimes >0 (restarting after processor
-            label& nNewCoarse,
+            //- agglomeration)
-            const lduAddressing& fineAddressing,
+            virtual void agglomerate
-            const labelUList& restriction,
+            (
-            const label nCoarse
+                const label nCellsInCoarsestLevel,
-        );
+                const label startLevel,
+                const scalarField& startFaceWeights,
+                const bool doProcessorAgglomerate = true
+            ) = 0;
+
+            //- Given restriction determines if coarse cells are connected.
+            //  Return ok is so, otherwise creates new restriction that is
+            static bool checkRestriction
+            (
+                labelList& newRestrict,
+                label& nNewCoarse,
+                const lduAddressing& fineAddressing,
+                const labelUList& restriction,
+                const label nCoarse
+            );
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/algebraicPairGAMGAgglomeration/algebraicPairGAMGAgglomeration.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,15 +55,19 @@ Foam::algebraicPairGAMGAgglomeration::algebraicPairGAMGAgglomeration
 :
     pairGAMGAgglomeration(matrix.mesh(), controlDict)
 {
-    const lduMesh& mesh = matrix.mesh();
-
     if (matrix.hasLower())
     {
-        agglomerate(mesh, max(mag(matrix.upper()), mag(matrix.lower())));
+        agglomerate
+        (
+            nCellsInCoarsestLevel_,
+            0,
+            max(mag(matrix.upper()), mag(matrix.lower())),
+            true
+        );
     }
     else
     {
-        agglomerate(mesh, mag(matrix.upper()));
+        agglomerate(nCellsInCoarsestLevel_, 0, mag(matrix.upper()), true);
     }
 }
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/dummyAgglomeration/dummyAgglomeration.C

@@ -74,7 +74,7 @@ Foam::dummyAgglomeration::dummyAgglomeration
     }
 
     // Shrink the storage of the levels to those created
-    compactLevels(nLevels_);
+    compactLevels(nLevels_, true);
 }
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/dummyAgglomeration/dummyAgglomeration.H

@@ -34,8 +34,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef dummyAgglomeration_H
+#ifndef Foam_dummyAgglomeration_H
-#define dummyAgglomeration_H
+#define Foam_dummyAgglomeration_H
 
 #include "GAMGAgglomeration.H"
 
@@ -52,7 +52,7 @@ class dummyAgglomeration
 :
     public GAMGAgglomeration
 {
-    // Private data
+    // Private Data
 
         //- Preset number of levels
         label nLevels_;
@@ -81,6 +81,20 @@ public:
             const lduMesh& mesh,
             const dictionary& controlDict
         );
+
+        //- Agglomerate from a starting level. Starting level is usually 0
+        //- (initial mesh) but sometimes >0 (restarting after processor
+        //- agglomeration)
+        virtual void agglomerate
+        (
+            const label nCellsInCoarsestLevel,
+            const label startLevel,
+            const scalarField& startFaceWeights,
+            const bool doProcessorAgglomerate = true
+        )
+        {
+            NotImplemented;
+        }
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomerate.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,43 +33,79 @@ License
 
 void Foam::pairGAMGAgglomeration::agglomerate
 (
-    const lduMesh& mesh,
+    const label nCellsInCoarsestLevel,
-    const scalarField& faceWeights
+    const label startLevel,
+    const scalarField& startFaceWeights,
+    const bool doProcessorAgglomerate
 )
 {
-    const label nFaces = mesh.lduAddr().lowerAddr().size();
+    if (nCells_.size() < maxLevels_)
-    if (faceWeights.size() != nFaces)
     {
-        FatalErrorInFunction
+        // See compactLevels. Make space if not enough
-            << "Supplied number of face weights " << faceWeights.size()
+        nCells_.resize(maxLevels_);
-            << " does not correspond to the number of faces " << nFaces << endl
+        restrictAddressing_.resize(maxLevels_);
-            << "This may be because of using a geometry-based"
+        nFaces_.resize(maxLevels_);
-            << " agglomeration method instead of a matrix-based one"
+        faceRestrictAddressing_.resize(maxLevels_);
-            << exit(FatalError);
+        faceFlipMap_.resize(maxLevels_);
+        nPatchFaces_.resize(maxLevels_);
+        patchFaceRestrictAddressing_.resize(maxLevels_);
+        meshLevels_.resize(maxLevels_);
+        // Have procCommunicator_ always, even if not procAgglomerating.
+        // Use value -1 to indicate nothing is proc-agglomerated
+        procCommunicator_.resize(maxLevels_ + 1, -1);
+        if (processorAgglomerate())
+        {
+            procAgglomMap_.resize(maxLevels_);
+            agglomProcIDs_.resize(maxLevels_);
+            procCommunicator_.resize(maxLevels_);
+            procCellOffsets_.resize(maxLevels_);
+            procFaceMap_.resize(maxLevels_);
+            procBoundaryMap_.resize(maxLevels_);
+            procBoundaryFaceMap_.resize(maxLevels_);
+        }
     }
 
 
     // Start geometric agglomeration from the given faceWeights
-    scalarField* faceWeightsPtr = const_cast<scalarField*>(&faceWeights);
+    scalarField faceWeights = startFaceWeights;
 
     // Agglomerate until the required number of cells in the coarsest level
     // is reached
 
     label nPairLevels = 0;
-    label nCreatedLevels = 0;
+    label nCreatedLevels = startLevel;
 
     while (nCreatedLevels < maxLevels_ - 1)
     {
+        if (!hasMeshLevel(nCreatedLevels))
+        {
+            FatalErrorInFunction<< "No mesh at nCreatedLevels:"
+                << nCreatedLevels
+                << exit(FatalError);
+        }
+
+        const auto& fineMesh = meshLevel(nCreatedLevels);
+
+
         label nCoarseCells = -1;
 
         tmp<labelField> finalAgglomPtr = agglomerate
         (
             nCoarseCells,
-            meshLevel(nCreatedLevels).lduAddr(),
+            fineMesh.lduAddr(),
-            *faceWeightsPtr
+            faceWeights
         );
 
-        if (continueAgglomerating(finalAgglomPtr().size(), nCoarseCells))
+        if
+        (
+            continueAgglomerating
+            (
+                nCellsInCoarsestLevel,
+                finalAgglomPtr().size(),
+                nCoarseCells,
+                fineMesh.comm()
+            )
+        )
         {
             nCells_[nCreatedLevels] = nCoarseCells;
             restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
@@ -78,32 +115,25 @@ void Foam::pairGAMGAgglomeration::agglomerate
             break;
         }
 
+        // Create coarse mesh
         agglomerateLduAddressing(nCreatedLevels);
 
         // Agglomerate the faceWeights field for the next level
         {
-            scalarField* aggFaceWeightsPtr
+            scalarField aggFaceWeights
             (
-                new scalarField
+                meshLevels_[nCreatedLevels].upperAddr().size(),
-                (
+                0.0
-                    meshLevels_[nCreatedLevels].upperAddr().size(),
-                    0.0
-                )
             );
 
             restrictFaceField
             (
-                *aggFaceWeightsPtr,
+                aggFaceWeights,
-                *faceWeightsPtr,
+                faceWeights,
                 nCreatedLevels
             );
 
-            if (nCreatedLevels)
+            faceWeights = std::move(aggFaceWeights);
-            {
-                delete faceWeightsPtr;
-            }
-
-            faceWeightsPtr = aggFaceWeightsPtr;
         }
 
         if (nPairLevels % mergeLevels_)
@@ -119,13 +149,7 @@ void Foam::pairGAMGAgglomeration::agglomerate
     }
 
     // Shrink the storage of the levels to those created
-    compactLevels(nCreatedLevels);
+    compactLevels(nCreatedLevels, doProcessorAgglomerate);
-
-    // Delete temporary geometry storage
-    if (nCreatedLevels)
-    {
-        delete faceWeightsPtr;
-    }
 }
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -66,13 +67,6 @@ protected:
 
     // Protected Member Functions
 
-        //- Agglomerate all levels starting from the given face weights
-        void agglomerate
-        (
-            const lduMesh& mesh,
-            const scalarField& faceWeights
-        );
-
         //- No copy construct
         pairGAMGAgglomeration(const pairGAMGAgglomeration&) = delete;
 
@@ -102,6 +96,17 @@ public:
             const lduAddressing& fineMatrixAddressing,
             const scalarField& faceWeights
         );
+
+        //- Agglomerate from a starting level. Starting level is usually 0
+        //- (initial mesh) but sometimes >0 (restarting after processor
+        //- agglomeration)
+        virtual void agglomerate
+        (
+            const label nCellsInCoarsestLevel,
+            const label startLevel,
+            const scalarField& startFaceWeights,
+            const bool doProcessorAgglomerate = true
+        );
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/GAMGProcAgglomeration/GAMGProcAgglomeration.C

@@ -51,7 +51,8 @@ void Foam::GAMGProcAgglomeration::printStats
     {
         if (agglom.hasMeshLevel(levelI))
         {
-            os  << agglom.meshLevel(levelI).info() << endl;
+            os  << "Level " << levelI << " mesh:"
+                << agglom.meshLevel(levelI).info() << endl;
         }
         else
         {

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/masterCoarsestGAMGProcAgglomeration/masterCoarsestGAMGProcAgglomeration.C

@@ -62,6 +62,10 @@ Foam::masterCoarsestGAMGProcAgglomeration::masterCoarsestGAMGProcAgglomeration
             0,
             keyType::LITERAL
         )
+    ),
+    nMasterCoarsestCells_
+    (
+        controlDict.getOrDefault<label>("nCellsInMasterLevel", -1)
     )
 {
     const auto* ePtr = controlDict.findEntry("nMasters", keyType::LITERAL);
@@ -208,9 +212,36 @@ bool Foam::masterCoarsestGAMGProcAgglomeration::agglomerate()
                         agglomProcIDs,
                         comms_.back()
                     );
+
+                    if (nMasterCoarsestCells_ > 0)
+                    {
+                        const label levelI = agglom_.size();
+                        if (agglom_.hasMeshLevel(levelI))
+                        {
+                            const lduMesh& fineMesh = agglom_.meshLevel(levelI);
+                            const auto& addr = fineMesh.lduAddr();
+                            const scalarField weights
+                            (
+                                addr.lowerAddr().size(),
+                                1.0
+                            );
+                            agglom_.agglomerate
+                            (
+                                nMasterCoarsestCells_,
+                                levelI,
+                                weights,
+                                false
+                            );
+                        }
+                    }
                 }
             }
         }
+
+
+        // Note that at this point for nMasterCoarsestCells_ the non-master
+        // processors will have less levels. This does/should not matter since
+        // they are not involved in those levels
     }
 
     // Print a bit

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/masterCoarsestGAMGProcAgglomeration/masterCoarsestGAMGProcAgglomeration.H

@@ -50,6 +50,11 @@ Description
 
             // Alternative : specify number of (equi-distributed) masters
             nMasters                2;
+
+            // Restart local agglomeration after processor agglomeration. Used
+            // as nCellsInCoarsestLevel. Only applicable for pair-wise local
+            // agglomeration methods.
+            nMasterCoarsestCells    1;
         }
 
 SourceFiles
@@ -67,6 +72,9 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
+class GAMGAgglomeration;
+
 /*---------------------------------------------------------------------------*\
               Class masterCoarsestGAMGProcAgglomeration Declaration
 \*---------------------------------------------------------------------------*/
@@ -79,6 +87,22 @@ class masterCoarsestGAMGProcAgglomeration
 
         label nProcessorsPerMaster_;
 
+        const label nMasterCoarsestCells_;
+
+        DynamicList<label> comms_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        masterCoarsestGAMGProcAgglomeration
+        (
+            const masterCoarsestGAMGProcAgglomeration&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const masterCoarsestGAMGProcAgglomeration&) = delete;
+
 
 public:
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -298,7 +299,7 @@ void Foam::GAMGSolver::gatherMatrices
     List<List<label>>& otherRanks
 ) const
 {
-    if (debug)
+    if (debug & 2)
     {
         Pout<< "GAMGSolver::gatherMatrices :"
             << " collecting matrices from procs:" << procIDs

src/finiteVolume/fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -67,7 +68,8 @@ Foam::faceAreaPairGAMGAgglomeration::faceAreaPairGAMGAgglomeration
     //agglomerate(mesh, sqrt(fvmesh.magSf().primitiveField()));
     agglomerate
     (
-        mesh,
+        nCellsInCoarsestLevel_,
+        0,          //mesh,
         mag
         (
             cmptMultiply
@@ -77,7 +79,8 @@ Foam::faceAreaPairGAMGAgglomeration::faceAreaPairGAMGAgglomeration
                 vector(1, 1.01, 1.02)
                 //vector::one
             )
-        )
+        ),
+        true
     );
 }
 
@@ -95,7 +98,8 @@ Foam::faceAreaPairGAMGAgglomeration::faceAreaPairGAMGAgglomeration
     //agglomerate(mesh, sqrt(mag(faceAreas)));
     agglomerate
     (
-        mesh,
+        nCellsInCoarsestLevel_,
+        0,          //mesh,
         mag
         (
             cmptMultiply
@@ -105,7 +109,8 @@ Foam::faceAreaPairGAMGAgglomeration::faceAreaPairGAMGAgglomeration
                 vector(1, 1.01, 1.02)
                 //vector::one
             )
-        )
+        ),
+        true
     );
 }
 

src/finiteVolume/lduPrimitiveMeshAssembly/assemblyFaceAreaPairGAMGAgglomeration/assemblyFaceAreaPairGAMGAgglomeration.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -181,7 +181,8 @@ assemblyFaceAreaPairGAMGAgglomeration
 
         agglomerate
         (
-            mesh,
+            nCellsInCoarsestLevel_,
+            0,      //mesh,
             mag
             (
                 cmptMultiply
@@ -189,7 +190,8 @@ assemblyFaceAreaPairGAMGAgglomeration
                     faceAreas/sqrt(mag(faceAreas)),
                     vector(1, 1.01, 1.02)
                 )
-            )
+            ),
+            true
         );
     }
     else
@@ -199,7 +201,8 @@ assemblyFaceAreaPairGAMGAgglomeration
 
         agglomerate
         (
-            matrix.mesh(),
+            nCellsInCoarsestLevel_,
+            0,      //matrix.mesh(),
             mag
             (
                 cmptMultiply
@@ -209,7 +212,8 @@ assemblyFaceAreaPairGAMGAgglomeration
                     vector(1, 1.01, 1.02)
                     //vector::one
                 )
-            )
+            ),
+            true
         );
     }
 }
@@ -228,7 +232,8 @@ assemblyFaceAreaPairGAMGAgglomeration
 {
     agglomerate
     (
-        matrix.mesh(),
+        nCellsInCoarsestLevel_,
+        0,      //matrix.mesh(),
         mag
         (
             cmptMultiply
@@ -236,7 +241,8 @@ assemblyFaceAreaPairGAMGAgglomeration
                 faceAreas/sqrt(mag(faceAreas)),
                 vector(1, 1.01, 1.02)
             )
-        )
+        ),
+        true
     );
 }
 

src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomeration.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -290,7 +290,16 @@ Foam::MGridGenGAMGAgglomeration::MGridGenGAMGAgglomeration
             );
         }
 
-        if (continueAgglomerating(finalAgglomPtr().size(), nCoarseCells))
+        if
+        (
+            continueAgglomerating
+            (
+                nCellsInCoarsestLevel_,
+                finalAgglomPtr().size(),
+                nCoarseCells,
+                meshLevel(nCreatedLevels).comm()
+            )
+        )
         {
             nCells_[nCreatedLevels] = nCoarseCells;
             restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
@@ -359,7 +368,7 @@ Foam::MGridGenGAMGAgglomeration::MGridGenGAMGAgglomeration
     }
 
     // Shrink the storage of the levels to those created
-    compactLevels(nCreatedLevels);
+    compactLevels(nCreatedLevels, true);
 
     // Delete temporary geometry storage
     if (nCreatedLevels)

tutorials/compressible/rhoSimpleFoam/squareBendLiq/Allrun

@@ -20,7 +20,11 @@ runParallel $(getApplication)
 foamDictionary -entry solvers.p.processorAgglomerator -set 'masterCoarsest' system/fvSolution
 runParallel -s masterCoarsest $(getApplication)
 
-#- Run with processorAgglomerator - multiple masters
+#- Run with processorAgglomerator+nCellsInMasterLevel
+foamDictionary -entry solvers.p.nCellsInMasterLevel -set '1' system/fvSolution
+runParallel -s nCellsInMasterLevel $(getApplication)
+
+#- Run with processorAgglomerator - multiple masters 
 foamDictionary -entry solvers.p.nMasters -set '2' system/fvSolution
 runParallel -s masterCoarsest2 $(getApplication)