ENH: snappyHexMesh: enable late balancing. Fixes #2792

etc/caseDicts/annotated/snappyHexMeshDict

@@ -100,8 +100,8 @@ castellatedMeshControls
     // ~~~~~~~~~~~~~~~~~~~~~
 
     // If local number of cells is >= maxLocalCells on any processor
-    // switches from from refinement followed by balancing
+    // algorithm switches from refinement followed by (weighted) balancing
-    // (current method) to (weighted) balancing before refinement.
+    // to (weighted) balancing followed by refinement.
     maxLocalCells 100000;
 
     // Overall cell limit (approximately). Refinement will stop immediately
@@ -123,9 +123,24 @@ castellatedMeshControls
     // Allow a certain level of imbalance during refining
     // (since balancing is quite expensive)
     // Expressed as fraction of perfect balance (= overall number of cells /
-    // nProcs). 0=balance always.
+    // nProcs). 0=balance always. Can be triggered by using maxCellUnbalance
     maxLoadUnbalance 0.10;
 
+    // Optional absolute unbalance start trigger. Used to avoid balancing
+    // in first refinement iterations since balancing is expensive for
+    // high core counts
+    // The trigger value is used for two checks:
+    // a: if on any processor the number of new added cells <= maxCellUnbalance
+    //    the balancing is skipped
+    // b: if on any processor the value 
+    //    (new cell count of proc - idealCellCount) <= maxCellUnbalance
+    //    the balancing is skipped
+    // Set to -1 to deactivate
+    //maxCellUnbalance 200000;
+
+    // Optional forcing balancing after the mesh is completely refined
+    //balanceAtEnd     false;
+
     // Number of buffer layers between different levels.
     // 1 means normal 2:1 refinement restriction, larger means slower
     // refinement.

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -1271,6 +1271,17 @@ public:
             //- Refine some cells
             autoPtr<mapPolyMesh> refine(const labelList& cellsToRefine);
 
+            //- Balance the mesh
+            autoPtr<mapDistributePolyMesh> balance
+            (
+                const string& msg,
+                decompositionMethod& decomposer,
+                fvMeshDistribute& distributor,
+                labelList& cellsToRefine,
+                const scalar maxLoadUnbalance,
+                const label maxCellUnbalance
+            );
+
             //- Refine some cells and rebalance
             autoPtr<mapDistributePolyMesh> refineAndBalance
             (
@@ -1278,7 +1289,8 @@ public:
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor,
                 const labelList& cellsToRefine,
-                const scalar maxLoadUnbalance
+                const scalar maxLoadUnbalance,
+                const label maxCellUnbalance
             );
 
             //- Balance before refining some cells
@@ -1288,7 +1300,8 @@ public:
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor,
                 const labelList& cellsToRefine,
-                const scalar maxLoadUnbalance
+                const scalar maxLoadUnbalance,
+                const label maxCellUnbalance
             );
 
             //- Calculate list of cells to directionally refine

src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -2598,6 +2598,120 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::refine
 }
 
 
+// Load balancing
+Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
+(
+    const string& msg,
+    decompositionMethod& decomposer,
+    fvMeshDistribute& distributor,
+    labelList& cellsToRefine,
+    const scalar maxLoadUnbalance,
+    const label maxCellUnbalance
+)
+{
+    autoPtr<mapDistributePolyMesh> distMap;
+
+    if (Pstream::nProcs() > 1)
+    {
+        // First check if we need to balance at all. Precalculate number of
+        // cells after refinement and see what maximum difference is.
+        const scalar nNewCells =
+            scalar(mesh_.nCells() + 7*cellsToRefine.size());
+        const scalar nNewCellsAll = returnReduce(nNewCells, sumOp<scalar>());
+        const scalar nIdealNewCells = nNewCellsAll / Pstream::nProcs();
+        const scalar unbalance = returnReduce
+        (
+            mag(1.0-nNewCells/nIdealNewCells),
+            maxOp<scalar>()
+        );
+
+        // Trigger the balancing to avoid too early balancing for better
+        // scaling performance.
+        const scalar nNewCellsOnly = scalar(7*cellsToRefine.size());
+
+        const label maxNewCells =
+            label(returnReduce(nNewCellsOnly, maxOp<scalar>()));
+
+        const label maxDeltaCells =
+            label(mag(returnReduce(nNewCells, maxOp<scalar>())-nIdealNewCells));
+
+        // New trigger to avoid too early balancing
+        // 1. Check if globally one proc exceeds the maxCellUnbalance value
+        //    related to the new added cells at the refinement loop
+        // 2. Check if globally one proc exceeds the maxCellUnbalance based on
+        //    the average cell count a proc should have
+        if
+        (
+            (maxNewCells <= maxCellUnbalance)
+         && (maxDeltaCells <= maxCellUnbalance)
+        )
+        {
+            Info<< "Skipping balancing since trigger value not reached:" << "\n"
+                << "    Trigger cell count: " << maxCellUnbalance << nl
+                << "    Max new cell count in proc: " << maxNewCells << nl
+                << "    Max difference between new cells and balanced: "
+                << maxDeltaCells << nl
+                << "    Max load unbalance " << maxLoadUnbalance
+                << nl <<endl;
+        }
+        else
+        {
+            if (unbalance <= maxLoadUnbalance)
+            {
+                Info<< "Skipping balancing since max unbalance " << unbalance
+                    << " is less than allowable " << maxLoadUnbalance
+                    << endl;
+            }
+            else
+            {
+                scalarField cellWeights(mesh_.nCells(), 1);
+                forAll(cellsToRefine, i)
+                {
+                    cellWeights[cellsToRefine[i]] += 7;
+                }
+
+                distMap = balance
+                (
+                    false,  //keepZoneFaces
+                    false,  //keepBaffles
+                    cellWeights,
+                    decomposer,
+                    distributor
+                );
+
+                // Update cells to refine
+                distMap().distributeCellIndices(cellsToRefine);
+
+                Info<< "Balanced mesh in = "
+                    << mesh_.time().cpuTimeIncrement() << " s" << endl;
+
+                printMeshInfo(debug, "After balancing " + msg);
+
+
+                if (debug&meshRefinement::MESH)
+                {
+                    Pout<< "Writing balanced " << msg
+                        << " mesh to time " << timeName() << endl;
+                    write
+                    (
+                        debugType(debug),
+                        writeType(writeLevel() | WRITEMESH),
+                        mesh_.time().path()/timeName()
+                    );
+                    Pout<< "Dumped debug data in = "
+                        << mesh_.time().cpuTimeIncrement() << " s" << endl;
+
+                    // test all is still synced across proc patches
+                    checkData();
+                }
+            }
+        }
+    }
+
+    return distMap;
+}
+
+
 // Do refinement of consistent set of cells followed by truncation and
 // load balancing.
 Foam::autoPtr<Foam::mapDistributePolyMesh>
@@ -2607,10 +2721,13 @@ Foam::meshRefinement::refineAndBalance
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor,
     const labelList& cellsToRefine,
-    const scalar maxLoadUnbalance
+    const scalar maxLoadUnbalance,
+    const label maxCellUnbalance
 )
 {
-    // Do all refinement
+    // Refinement
+    // ~~~~~~~~~~
+
     refine(cellsToRefine);
 
     if (debug&meshRefinement::MESH)
@@ -2638,63 +2755,17 @@ Foam::meshRefinement::refineAndBalance
     // Load balancing
     // ~~~~~~~~~~~~~~
 
-    autoPtr<mapDistributePolyMesh> distMap;
+    labelList noCellsToRefine;
 
-    if (Pstream::nProcs() > 1)
+    auto distMap = balance
-    {
+    (
-        scalar nIdealCells =
+        msg,
-            mesh_.globalData().nTotalCells()
+        decomposer,
-          / Pstream::nProcs();
+        distributor,
-
+        noCellsToRefine,    // mesh is already refined; no need to predict
-        scalar unbalance = returnReduce
+        maxLoadUnbalance,
-        (
+        maxCellUnbalance
-            mag(1.0-mesh_.nCells()/nIdealCells),
+    );
-            maxOp<scalar>()
-        );
-
-        if (unbalance <= maxLoadUnbalance)
-        {
-            Info<< "Skipping balancing since max unbalance " << unbalance
-                << " is less than allowable " << maxLoadUnbalance
-                << endl;
-        }
-        else
-        {
-            scalarField cellWeights(mesh_.nCells(), 1);
-
-            distMap = balance
-            (
-                false,  //keepZoneFaces
-                false,  //keepBaffles
-                cellWeights,
-                decomposer,
-                distributor
-            );
-
-            Info<< "Balanced mesh in = "
-                << mesh_.time().cpuTimeIncrement() << " s" << endl;
-
-            printMeshInfo(debug, "After balancing " + msg);
-
-
-            if (debug&meshRefinement::MESH)
-            {
-                Pout<< "Writing balanced " << msg
-                    << " mesh to time " << timeName() << endl;
-                write
-                (
-                    debugType(debug),
-                    writeType(writeLevel() | WRITEMESH),
-                    mesh_.time().path()/timeName()
-                );
-                Pout<< "Dumped debug data in = "
-                    << mesh_.time().cpuTimeIncrement() << " s" << endl;
-
-                // test all is still synced across proc patches
-                checkData();
-            }
-        }
-    }
 
     return distMap;
 }
@@ -2708,7 +2779,8 @@ Foam::meshRefinement::balanceAndRefine
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor,
     const labelList& initCellsToRefine,
-    const scalar maxLoadUnbalance
+    const scalar maxLoadUnbalance,
+    const label maxCellUnbalance
 )
 {
     labelList cellsToRefine(initCellsToRefine);
@@ -2737,90 +2809,24 @@ Foam::meshRefinement::balanceAndRefine
     //}
 
 
+
     // Load balancing
     // ~~~~~~~~~~~~~~
 
-    autoPtr<mapDistributePolyMesh> distMap;
+    auto distMap = balance
-
+    (
-    if (Pstream::nProcs() > 1)
+        msg,
-    {
+        decomposer,
-        // First check if we need to balance at all. Precalculate number of
+        distributor,
-        // cells after refinement and see what maximum difference is.
+        cellsToRefine,
-        scalar nNewCells = scalar(mesh_.nCells() + 7*cellsToRefine.size());
+        maxLoadUnbalance,
-        scalar nIdealNewCells =
+        maxCellUnbalance
-            returnReduce(nNewCells, sumOp<scalar>())
+    );
-          / Pstream::nProcs();
-        scalar unbalance = returnReduce
-        (
-            mag(1.0-nNewCells/nIdealNewCells),
-            maxOp<scalar>()
-        );
-
-        if (unbalance <= maxLoadUnbalance)
-        {
-            Info<< "Skipping balancing since max unbalance " << unbalance
-                << " is less than allowable " << maxLoadUnbalance
-                << endl;
-        }
-        else
-        {
-            scalarField cellWeights(mesh_.nCells(), 1);
-            forAll(cellsToRefine, i)
-            {
-                cellWeights[cellsToRefine[i]] += 7;
-            }
-
-            distMap = balance
-            (
-                false,  //keepZoneFaces
-                false,  //keepBaffles
-                cellWeights,
-                decomposer,
-                distributor
-            );
-
-            // Update cells to refine
-            distMap().distributeCellIndices(cellsToRefine);
-
-            Info<< "Balanced mesh in = "
-                << mesh_.time().cpuTimeIncrement() << " s" << endl;
-        }
-
-        //{
-        //    globalIndex globalCells(mesh_.nCells());
-        //
-        //    Info<< "** Distribution after balancing:" << endl;
-        //    for (const int procI : Pstream::allProcs())
-        //    {
-        //        Info<< "    " << procI << '\t'
-        //            << globalCells.localSize(procI) << endl;
-        //    }
-        //    Info<< endl;
-        //}
-
-        printMeshInfo(debug, "After balancing " + msg);
-
-        if (debug&meshRefinement::MESH)
-        {
-            Pout<< "Writing balanced " << msg
-                << " mesh to time " << timeName() << endl;
-            write
-            (
-                debugType(debug),
-                writeType(writeLevel() | WRITEMESH),
-                mesh_.time().path()/timeName()
-            );
-            Pout<< "Dumped debug data in = "
-                << mesh_.time().cpuTimeIncrement() << " s" << endl;
-
-            // test all is still synced across proc patches
-            checkData();
-        }
-    }
 
 
     // Refinement
     // ~~~~~~~~~~
+    // Note: uses updated cellsToRefine
 
     refine(cellsToRefine);
 

src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2015-2020,2022 OpenCFD Ltd.
+    Copyright (C) 2015-2020,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -85,6 +85,7 @@ Foam::refinementParameters::refinementParameters
         dict.getOrDefault("useTopologicalSnapDetection", true)
     ),
     maxLoadUnbalance_(dict.getOrDefault<scalar>("maxLoadUnbalance", 0)),
+    maxCellUnbalance_(dict.getOrDefault<label>("maxCellUnbalance", -1)),
     handleSnapProblems_
     (
         dict.getOrDefault<Switch>("handleSnapProblems", true)
@@ -97,6 +98,10 @@ Foam::refinementParameters::refinementParameters
     nFilterIter_(dict.getOrDefault<label>("nFilterIter", 2)),
     minCellFraction_(dict.getOrDefault<scalar>("minCellFraction", 0)),
     nMinCells_(dict.getOrDefault<label>("nMinCells", 0)),
+    balanceAtEnd_
+    (
+        dict.getOrDefault("balanceAtEnd", false)
+    ),
     dryRun_(dryRun)
 {
     point locationInMesh;

src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2014 OpenFOAM Foundation
-    Copyright (C) 2015-2020,2022 OpenCFD Ltd.
+    Copyright (C) 2015-2020,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -108,6 +108,9 @@ class refinementParameters
         //- Allowed load unbalance
         const scalar maxLoadUnbalance_;
 
+        //- Trigger cell count to start balancing
+        const label maxCellUnbalance_;
+
         const Switch handleSnapProblems_;
 
         const Switch interfaceRefine_;
@@ -120,6 +123,9 @@ class refinementParameters
 
         const label nMinCells_;
 
+        //- Force final balancing after castellation
+        const bool balanceAtEnd_;
+
         const bool dryRun_;
 
 
@@ -228,6 +234,12 @@ public:
                 return maxLoadUnbalance_;
             }
 
+            //- Trigger cell count to start balancing
+            label maxCellUnbalance() const
+            {
+                return maxCellUnbalance_;
+            }
+
             bool handleSnapProblems() const
             {
                 return handleSnapProblems_;
@@ -266,6 +278,12 @@ public:
                 return nMinCells_;
             }
 
+            //- Force final balancing after castallation
+            bool balanceAtEnd() const
+            {
+                return balanceAtEnd_;
+            }
+
 
         // Other
 

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C

@@ -179,7 +179,8 @@ Foam::label Foam::snappyRefineDriver::featureEdgeRefine
                     decomposer_,
                     distributor_,
                     cellsToRefine,
-                    refineParams.maxLoadUnbalance()
+                    refineParams.maxLoadUnbalance(),
+                    refineParams.maxCellUnbalance()
                 );
             }
             else
@@ -190,7 +191,8 @@ Foam::label Foam::snappyRefineDriver::featureEdgeRefine
                     decomposer_,
                     distributor_,
                     cellsToRefine,
-                    refineParams.maxLoadUnbalance()
+                    refineParams.maxLoadUnbalance(),
+                    refineParams.maxCellUnbalance()
                 );
             }
         }
@@ -313,7 +315,8 @@ Foam::label Foam::snappyRefineDriver::smallFeatureRefine
                 decomposer_,
                 distributor_,
                 cellsToRefine,
-                refineParams.maxLoadUnbalance()
+                refineParams.maxLoadUnbalance(),
+                refineParams.maxCellUnbalance()
             );
         }
         else
@@ -324,7 +327,8 @@ Foam::label Foam::snappyRefineDriver::smallFeatureRefine
                 decomposer_,
                 distributor_,
                 cellsToRefine,
-                refineParams.maxLoadUnbalance()
+                refineParams.maxLoadUnbalance(),
+                refineParams.maxCellUnbalance()
             );
         }
     }
@@ -490,7 +494,8 @@ Foam::label Foam::snappyRefineDriver::surfaceOnlyRefine
                 decomposer_,
                 distributor_,
                 cellsToRefine,
-                refineParams.maxLoadUnbalance()
+                refineParams.maxLoadUnbalance(),
+                refineParams.maxCellUnbalance()
             );
         }
         else
@@ -501,7 +506,8 @@ Foam::label Foam::snappyRefineDriver::surfaceOnlyRefine
                 decomposer_,
                 distributor_,
                 cellsToRefine,
-                refineParams.maxLoadUnbalance()
+                refineParams.maxLoadUnbalance(),
+                refineParams.maxCellUnbalance()
             );
         }
     }
@@ -742,7 +748,8 @@ Foam::label Foam::snappyRefineDriver::gapOnlyRefine
                 decomposer_,
                 distributor_,
                 cellsToRefine,
-                refineParams.maxLoadUnbalance()
+                refineParams.maxLoadUnbalance(),
+                refineParams.maxCellUnbalance()
             );
         }
         else
@@ -753,7 +760,8 @@ Foam::label Foam::snappyRefineDriver::gapOnlyRefine
                 decomposer_,
                 distributor_,
                 cellsToRefine,
-                refineParams.maxLoadUnbalance()
+                refineParams.maxLoadUnbalance(),
+                refineParams.maxCellUnbalance()
             );
         }
     }
@@ -967,7 +975,8 @@ Foam::label Foam::snappyRefineDriver::bigGapOnlyRefine
                 decomposer_,
                 distributor_,
                 cellsToRefine,
-                refineParams.maxLoadUnbalance()
+                refineParams.maxLoadUnbalance(),
+                refineParams.maxCellUnbalance()
             );
         }
         else
@@ -978,7 +987,8 @@ Foam::label Foam::snappyRefineDriver::bigGapOnlyRefine
                 decomposer_,
                 distributor_,
                 cellsToRefine,
-                refineParams.maxLoadUnbalance()
+                refineParams.maxLoadUnbalance(),
+                refineParams.maxCellUnbalance()
             );
         }
     }
@@ -1090,6 +1100,22 @@ Foam::label Foam::snappyRefineDriver::danglingCellRefine
         {
             Info<< "Stopping refining since too few cells selected."
                 << nl << endl;
+
+            if (refineParams.balanceAtEnd())
+            {
+                Info<< "Final mesh balancing" << endl;
+
+                meshRefiner_.balance
+                (
+                    "",
+                    decomposer_,
+                    distributor_,
+                    cellsToRefine,
+                    0,
+                    -1
+                );
+            }
+
             break;
         }
 
@@ -1108,7 +1134,8 @@ Foam::label Foam::snappyRefineDriver::danglingCellRefine
                 decomposer_,
                 distributor_,
                 cellsToRefine,
-                refineParams.maxLoadUnbalance()
+                refineParams.maxLoadUnbalance(),
+                refineParams.maxCellUnbalance()
             );
         }
         else
@@ -1119,7 +1146,8 @@ Foam::label Foam::snappyRefineDriver::danglingCellRefine
                 decomposer_,
                 distributor_,
                 cellsToRefine,
-                refineParams.maxLoadUnbalance()
+                refineParams.maxLoadUnbalance(),
+                refineParams.maxCellUnbalance()
             );
         }
     }
@@ -1453,7 +1481,8 @@ Foam::label Foam::snappyRefineDriver::refinementInterfaceRefine
                     decomposer_,
                     distributor_,
                     cellsToRefine,
-                    refineParams.maxLoadUnbalance()
+                    refineParams.maxLoadUnbalance(),
+                    refineParams.maxCellUnbalance()
                 );
             }
             else
@@ -1464,7 +1493,8 @@ Foam::label Foam::snappyRefineDriver::refinementInterfaceRefine
                     decomposer_,
                     distributor_,
                     cellsToRefine,
-                    refineParams.maxLoadUnbalance()
+                    refineParams.maxLoadUnbalance(),
+                    refineParams.maxCellUnbalance()
                 );
             }
         }
@@ -1701,7 +1731,8 @@ Foam::label Foam::snappyRefineDriver::boundaryRefinementInterfaceRefine
                     decomposer_,
                     distributor_,
                     cellsToRefine,
-                    refineParams.maxLoadUnbalance()
+                    refineParams.maxLoadUnbalance(),
+                    refineParams.maxCellUnbalance()
                 );
             }
             else
@@ -1712,7 +1743,8 @@ Foam::label Foam::snappyRefineDriver::boundaryRefinementInterfaceRefine
                     decomposer_,
                     distributor_,
                     cellsToRefine,
-                    refineParams.maxLoadUnbalance()
+                    refineParams.maxLoadUnbalance(),
+                    refineParams.maxCellUnbalance()
                 );
             }
         }
@@ -1961,7 +1993,8 @@ Foam::label Foam::snappyRefineDriver::shellRefine
                 decomposer_,
                 distributor_,
                 cellsToRefine,
-                refineParams.maxLoadUnbalance()
+                refineParams.maxLoadUnbalance(),
+                refineParams.maxCellUnbalance()
             );
         }
         else
@@ -1972,7 +2005,8 @@ Foam::label Foam::snappyRefineDriver::shellRefine
                 decomposer_,
                 distributor_,
                 cellsToRefine,
-                refineParams.maxLoadUnbalance()
+                refineParams.maxLoadUnbalance(),
+                refineParams.maxCellUnbalance()
             );
         }
     }