pairGAMGAgglomeration: Support reverse cell-scanning without reversing the...

pairGAMGAgglomeration: Support reverse cell-scanning without reversing the order of the generated agglomeration
Gives better cache coherency

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

@@ -26,7 +26,96 @@ License
 #include "pairGAMGAgglomeration.H"
 #include "lduAddressing.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+void Foam::pairGAMGAgglomeration::agglomerate
+(
+    const lduMesh& mesh,
+    const scalarField& faceWeights
+)
+{
+    // Start geometric agglomeration from the given faceWeights
+    scalarField* faceWeightsPtr = const_cast<scalarField*>(&faceWeights);
+
+    // Agglomerate until the required number of cells in the coarsest level
+    // is reached
+
+    label nPairLevels = 0;
+    label nCreatedLevels = 0;
+
+    while (nCreatedLevels < maxLevels_ - 1)
+    {
+        label nCoarseCells = -1;
+
+        tmp<labelField> finalAgglomPtr = agglomerate
+        (
+            nCoarseCells,
+            meshLevel(nCreatedLevels).lduAddr(),
+            *faceWeightsPtr
+        );
+
+        if (continueAgglomerating(nCoarseCells))
+        {
+            nCells_[nCreatedLevels] = nCoarseCells;
+            restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
+        }
+        else
+        {
+            break;
+        }
+
+        agglomerateLduAddressing(nCreatedLevels);
+
+        // Agglomerate the faceWeights field for the next level
+        {
+            scalarField* aggFaceWeightsPtr
+            (
+                new scalarField
+                (
+                    meshLevels_[nCreatedLevels].upperAddr().size(),
+                    0.0
+                )
+            );
+
+            restrictFaceField
+            (
+                *aggFaceWeightsPtr,
+                *faceWeightsPtr,
+                nCreatedLevels
+            );
+
+            if (nCreatedLevels)
+            {
+                delete faceWeightsPtr;
+            }
+
+            faceWeightsPtr = aggFaceWeightsPtr;
+        }
+
+        if (nPairLevels % mergeLevels_)
+        {
+            combineLevels(nCreatedLevels);
+        }
+        else
+        {
+            nCreatedLevels++;
+        }
+
+        nPairLevels++;
+    }
+
+    // Shrink the storage of the levels to those created
+    compactLevels(nCreatedLevels);
+
+    // Delete temporary geometry storage
+    if (nCreatedLevels)
+    {
+        delete faceWeightsPtr;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::labelField> Foam::pairGAMGAgglomeration::agglomerate
 (
@@ -95,9 +184,13 @@ Foam::tmp<Foam::labelField> Foam::pairGAMGAgglomeration::agglomerate
     labelField& coarseCellMap = tcoarseCellMap();
 
     nCoarseCells = 0;
+    label celli;
 
-    for (label celli=0; celli<nFineCells; celli++)
+    for (label cellfi=0; cellfi<nFineCells; cellfi++)
     {
+        // Change cell ordering depending on direction for this level
+        celli = forward_ ? cellfi : nFineCells - cellfi - 1;
+
         if (coarseCellMap[celli] < 0)
         {
             label matchFaceNo = -1;
@@ -171,11 +264,13 @@ Foam::tmp<Foam::labelField> Foam::pairGAMGAgglomeration::agglomerate
         }
     }
 
-
     // Check that all cells are part of clusters,
     // if not create single-cell "clusters" for each
-    for (label celli=0; celli<nFineCells; celli++)
+    for (label cellfi=0; cellfi<nFineCells; cellfi++)
     {
+        // Change cell ordering depending on direction for this level
+        celli = forward_ ? cellfi : nFineCells - cellfi - 1;
+
         if (coarseCellMap[celli] < 0)
         {
             coarseCellMap[celli] = nCoarseCells;
@@ -183,105 +278,23 @@ Foam::tmp<Foam::labelField> Foam::pairGAMGAgglomeration::agglomerate
         }
     }
 
-    // Reverse the map ordering to improve the next level of agglomeration
-    // (doesn't always help and is sometimes detrimental)
-    nCoarseCells--;
-
-    forAll(coarseCellMap, celli)
-    {
-        coarseCellMap[celli] = nCoarseCells - coarseCellMap[celli];
-    }
-
-    nCoarseCells++;
-
-    return tcoarseCellMap;
-}
-
-
-void Foam::pairGAMGAgglomeration::agglomerate
-(
-    const lduMesh& mesh,
-    const scalarField& faceWeights
-)
-{
-    // Start geometric agglomeration from the given faceWeights
-    scalarField* faceWeightsPtr = const_cast<scalarField*>(&faceWeights);
-
-    // Agglomerate until the required number of cells in the coarsest level
-    // is reached
-
-    label nPairLevels = 0;
-    label nCreatedLevels = 0;
-
-    while (nCreatedLevels < maxLevels_ - 1)
+    if (!forward_)
     {
-        label nCoarseCells = -1;
-
-        tmp<labelField> finalAgglomPtr = agglomerate
-        (
-            nCoarseCells,
-            meshLevel(nCreatedLevels).lduAddr(),
-            *faceWeightsPtr
-        );
-
-        if (continueAgglomerating(nCoarseCells))
-        {
-            nCells_[nCreatedLevels] = nCoarseCells;
-            restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
-        }
-        else
-        {
-            break;
-        }
-
-        agglomerateLduAddressing(nCreatedLevels);
-
-        // Agglomerate the faceWeights field for the next level
-        {
-            scalarField* aggFaceWeightsPtr
-            (
-                new scalarField
-                (
-                    meshLevels_[nCreatedLevels].upperAddr().size(),
-                    0.0
-                )
-            );
-
-            restrictFaceField
-            (
-                *aggFaceWeightsPtr,
-                *faceWeightsPtr,
-                nCreatedLevels
-            );
-
-            if (nCreatedLevels)
-            {
-                delete faceWeightsPtr;
-            }
+        nCoarseCells--;
 
-            faceWeightsPtr = aggFaceWeightsPtr;
-        }
-
-        if (nPairLevels % mergeLevels_)
-        {
-            combineLevels(nCreatedLevels);
-        }
-        else
+        forAll(coarseCellMap, celli)
         {
-            nCreatedLevels++;
+            coarseCellMap[celli] = nCoarseCells - coarseCellMap[celli];
         }
 
-        nPairLevels++;
+        nCoarseCells++;
     }
 
-    // Shrink the storage of the levels to those created
-    compactLevels(nCreatedLevels);
+    // Reverse the map ordering for the next level
+    // to improve the next level of agglomeration
+    forward_ = !forward_;
 
-    // Delete temporary geometry storage
-    if (nCreatedLevels)
-    {
-        delete faceWeightsPtr;
-    }
+    return tcoarseCellMap;
 }
 
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,7 +29,8 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(pairGAMGAgglomeration, 0);
+    defineTypeNameAndDebug(pairGAMGAgglomeration, 0);
+    bool pairGAMGAgglomeration::forward_(true);
 }
 
 

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

@@ -56,6 +56,9 @@ class pairGAMGAgglomeration
         //- Number of levels to merge, 1 = don't merge, 2 = merge pairs etc.
         label mergeLevels_;
 
+        //- Direction of cell loop for the current level
+        static bool forward_;
+
 
 protected:
 
@@ -97,7 +100,6 @@ public:
             const lduAddressing& fineMatrixAddressing,
             const scalarField& faceWeights
         );
-
 };