From 46e1b00c34217a19fdf4574263c41263be8fea3d Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Fri, 15 Mar 2024 19:42:39 +0100
Subject: [PATCH] ENH: handle sub-mesh connectivity by subsetting of adjacency
 matrix

- in renumberMesh replace calculation of a subMesh connectivity
  with calculation of the full mesh connectivity followed by subsetting
  of the full adjacency matrix. This should reduce the overall number of
  operations. (MR !669)
---
 .../manipulation/renumberMesh/renumberMesh.C  | 245 ++++++++++++------
 .../renumberMesh/subsetAdjacency.H            | 168 ++++++++++++
 2 files changed, 327 insertions(+), 86 deletions(-)
 create mode 100644 applications/utilities/mesh/manipulation/renumberMesh/subsetAdjacency.H

diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
index 64011258c6..3702cfe8ce 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
@@ -145,6 +145,7 @@ Usage
 #include "hexRef8Data.H"
 #include "regionProperties.H"
 #include "polyMeshTools.H"
+#include "subsetAdjacency.H"
 
 using namespace Foam;
 
@@ -194,21 +195,56 @@ tmp<volScalarField> createScalarField
 }
 
 
-// Calculate band of matrix
-label getBand(const labelList& owner, const labelList& neighbour)
+// Calculate band of mesh
+// label getBand(const labelUList& owner, const labelUList& neighbour)
+// {
+//     label bandwidth = 0;
+//
+//     forAll(neighbour, facei)
+//     {
+//         const label width = neighbour[facei] - owner[facei];
+//
+//         if (bandwidth < width)
+//         {
+//             bandwidth = width;
+//         }
+//     }
+//     return bandwidth;
+// }
+
+
+// Calculate band and profile of matrix. Profile is scalar to avoid overflow
+Tuple2<label, scalar> getBand
+(
+    const CompactListList<label>& mat
+)
 {
-    label band = 0;
+    Tuple2<label, scalar> metrics(0, 0);
 
-    forAll(neighbour, facei)
+    auto& bandwidth = metrics.first();
+    auto& profile = metrics.second();
+
+    forAll(mat, celli)
     {
-        label diff = neighbour[facei] - owner[facei];
+        const auto& neighbours = mat[celli];
 
-        if (diff > band)
+        const label nNbr = neighbours.size();
+
+        if (nNbr)
         {
-            band = diff;
+            // Max distance
+            const label width = (neighbours[nNbr-1] - celli);
+
+            if (bandwidth < width)
+            {
+                bandwidth = width;
+            }
+
+            profile += scalar(width);
         }
     }
-    return band;
+
+    return metrics;
 }
 
 
@@ -217,27 +253,35 @@ void getBand
 (
     const bool calculateIntersect,
     const label nCells,
-    const labelList& owner,
-    const labelList& neighbour,
+    const labelUList& owner,
+    const labelUList& neighbour,
     label& bandwidth,
     scalar& profile,            // scalar to avoid overflow
     scalar& sumSqrIntersect     // scalar to avoid overflow
 )
 {
     labelList cellBandwidth(nCells, Foam::zero{});
-    scalarField nIntersect(nCells, Foam::zero{});
+
+    bandwidth = 0;
 
     forAll(neighbour, facei)
     {
-        label own = owner[facei];
-        label nei = neighbour[facei];
+        const label own = owner[facei];
+        const label nei = neighbour[facei];
 
         // Note: mag not necessary for correct (upper-triangular) ordering.
-        label diff = nei-own;
-        cellBandwidth[nei] = max(cellBandwidth[nei], diff);
-    }
+        const label width = nei - own;
 
-    bandwidth = max(cellBandwidth);
+        if (cellBandwidth[nei] < width)
+        {
+            cellBandwidth[nei] = width;
+
+            if (bandwidth < width)
+            {
+                bandwidth = width;
+            }
+        }
+    }
 
     // Do not use field algebra because of conversion label to scalar
     profile = 0;
@@ -246,14 +290,16 @@ void getBand
         profile += scalar(width);
     }
 
-    sumSqrIntersect = 0.0;
+    sumSqrIntersect = 0;
     if (calculateIntersect)
     {
+        scalarField nIntersect(nCells, Foam::zero{});
+
         forAll(nIntersect, celli)
         {
             for (label colI = celli-cellBandwidth[celli]; colI <= celli; colI++)
             {
-                nIntersect[colI] += 1.0;
+                nIntersect[colI] += scalar(1);
             }
         }
 
@@ -675,9 +721,8 @@ CompactListList<label> regionRenumber
 
         forAll(regionCellOrder, regioni)
         {
-            // Info<< "    region " << regioni
-            //     << " starts at " << regionCellOrder.localStart(regioni)
-            //     << nl;
+            // Info<< "    region " << regioni << " starts at "
+            //     << regionCellOrder.localStart(regioni) << nl;
 
             // No parallel communication
             const bool oldParRun = UPstream::parRun(false);
@@ -703,27 +748,42 @@ CompactListList<label> regionRenumber
     {
         timer.resetTimeIncrement();
 
+        // Create adjacency matrix of the full mesh and subset subsequently.
+        // This is more efficient than creating adjacency matrices of
+        // sub-meshes.
+
+        // No parallel communication
+        const bool oldParRun = UPstream::parRun(false);
+
+        // The local connectivity of the full (non-subsetted) mesh
+        CompactListList<label> meshCellCells;
+        globalMeshData::calcCellCells(mesh, meshCellCells);
+        UPstream::parRun(oldParRun);  // Restore parallel state
+
+        timings[TimingType::CELL_CELLS] += timer.timeIncrement();
+
+        // For the respective subMesh selections
+        bitSet subsetCells(mesh.nCells());
+
         forAll(regionCellOrder, regioni)
         {
-            // Info<< "    region " << regioni
-            //     << " starts at " << regionCellOrder.localStart(regioni)
-            //     << nl;
+            // Info<< "    region " << regioni << " starts at "
+            //     << regionCellOrder.localStart(regioni) << nl;
+
+            subsetCells = false;
+            subsetCells.set(regionCellOrder[regioni]);
+
+            // Connectivity of local sub-mesh
+            labelList cellMap;
+            CompactListList<label> subCellCells =
+                subsetAdjacency(subsetCells, meshCellCells, cellMap);
+
+            timings[TimingType::CELL_CELLS] += timer.timeIncrement();
 
             // No parallel communication
             const bool oldParRun = UPstream::parRun(false);
 
-            // Connectivity of local sub-mesh
-            CompactListList<label> cellCells;
-            labelList cellMap = globalMeshData::calcCellCells
-            (
-                mesh,
-                regionCellOrder[regioni],
-                cellCells
-            );
-
-            timings[TimingType::CELL_CELLS] += timer.timeIncrement();
-
-            labelList subCellOrder = method.renumber(cellCells);
+            labelList subCellOrder = method.renumber(subCellCells);
 
             UPstream::parRun(oldParRun);  // Restore parallel state
 
@@ -835,10 +895,10 @@ int main(int argc, char *argv[])
     const bool dryrun = args.dryRun();
 
     const bool readDict = args.found("dict");
-    const bool doFrontWidth = args.found("frontWidth");
+    const bool doDecompose = args.found("decompose");
     const bool overwrite = args.found("overwrite");
     const bool doFields = !args.found("no-fields");
-    const bool doDecompose = args.found("decompose");
+    const bool doFrontWidth = args.found("frontWidth") && !doDecompose;
 
     word renumberMethodName;
     args.readIfPresent("renumber-method", renumberMethodName);
@@ -846,8 +906,7 @@ int main(int argc, char *argv[])
     if (doDecompose && UPstream::parRun())
     {
         FatalErrorIn(args.executable())
-            << "Cannot use -decompose option in parallel"
-            << " ... giving up" << nl
+            << "Cannot use -decompose option in parallel ... giving up" << nl
             << exit(FatalError);
     }
 
@@ -908,21 +967,21 @@ int main(int argc, char *argv[])
 
         reduce(band, maxOp<label>());
         reduce(profile, sumOp<scalar>());
-        reduce(sumSqrIntersect, sumOp<scalar>());
-
-        scalar rmsFrontwidth = Foam::sqrt
-        (
-            sumSqrIntersect/mesh.globalData().nTotalCells()
-        );
 
         Info<< "Mesh " << mesh.name()
             << " size: " << mesh.globalData().nTotalCells() << nl
-            << "Before renumbering :" << nl
+            << "Before renumbering" << nl
             << "    band           : " << band << nl
             << "    profile        : " << profile << nl;
 
         if (doFrontWidth)
         {
+            reduce(sumSqrIntersect, sumOp<scalar>());
+            scalar rmsFrontwidth = Foam::sqrt
+            (
+                sumSqrIntersect/mesh.globalData().nTotalCells()
+            );
+
             Info<< "    rms frontwidth : " << rmsFrontwidth << nl;
         }
 
@@ -1091,10 +1150,7 @@ int main(int argc, char *argv[])
         );
 
 
-        // List of objects read from time directory
-        // List of stored objects to clear from mesh
-
-        IOobjectList objects;
+        // List of stored objects to clear from mesh (after reading)
         DynamicList<regIOobject*> storedObjects;
 
         if (!dryrun && doFields)
@@ -1103,55 +1159,60 @@ int main(int argc, char *argv[])
 
             timer.resetTimeIncrement();
 
-            objects = IOobjectList(mesh, runTime.timeName());
+            IOobjectList objects(mesh, runTime.timeName());
             storedObjects.reserve(objects.size());
 
             const predicates::always nameMatcher;
 
             // Read GeometricFields
 
-            #undef  ReadFields
-            #define ReadFields(FieldType)                                     \
+            #undef  doLocalCode
+            #define doLocalCode(FieldType)                                    \
             readFields<FieldType>(mesh, objects, nameMatcher, storedObjects);
 
             // Read volume fields
-            ReadFields(volScalarField);
-            ReadFields(volVectorField);
-            ReadFields(volSphericalTensorField);
-            ReadFields(volSymmTensorField);
-            ReadFields(volTensorField);
+            doLocalCode(volScalarField);
+            doLocalCode(volVectorField);
+            doLocalCode(volSphericalTensorField);
+            doLocalCode(volSymmTensorField);
+            doLocalCode(volTensorField);
 
             // Read internal fields
-            ReadFields(volScalarField::Internal);
-            ReadFields(volVectorField::Internal);
-            ReadFields(volSphericalTensorField::Internal);
-            ReadFields(volSymmTensorField::Internal);
-            ReadFields(volTensorField::Internal);
+            doLocalCode(volScalarField::Internal);
+            doLocalCode(volVectorField::Internal);
+            doLocalCode(volSphericalTensorField::Internal);
+            doLocalCode(volSymmTensorField::Internal);
+            doLocalCode(volTensorField::Internal);
 
             // Read surface fields
-            ReadFields(surfaceScalarField);
-            ReadFields(surfaceVectorField);
-            ReadFields(surfaceSphericalTensorField);
-            ReadFields(surfaceSymmTensorField);
-            ReadFields(surfaceTensorField);
+            doLocalCode(surfaceScalarField);
+            doLocalCode(surfaceVectorField);
+            doLocalCode(surfaceSphericalTensorField);
+            doLocalCode(surfaceSymmTensorField);
+            doLocalCode(surfaceTensorField);
 
 
             // Read point fields
             const pointMesh& pMesh = pointMesh::New(mesh);
-            #undef  ReadPointFields
-            #define ReadPointFields(FieldType)                                \
+            #undef  doLocalCode
+            #define doLocalCode(FieldType)                                    \
             readFields<FieldType>(pMesh, objects, nameMatcher, storedObjects);
 
-            ReadPointFields(pointScalarField);
-            ReadPointFields(pointVectorField);
-            ReadPointFields(pointSphericalTensorField);
-            ReadPointFields(pointSymmTensorField);
-            ReadPointFields(pointTensorField);
+            doLocalCode(pointScalarField);
+            doLocalCode(pointVectorField);
+            doLocalCode(pointSphericalTensorField);
+            doLocalCode(pointSymmTensorField);
+            doLocalCode(pointTensorField);
 
-            #undef ReadFields
-            #undef ReadPointFields
+            #undef doLocalCode
 
             timings[TimingType::READ_FIELDS] += timer.timeIncrement();
+
+            // Write loaded fields when mesh.write() is called
+            for (auto* fldptr : storedObjects)
+            {
+                fldptr->writeOpt(IOobject::AUTO_WRITE);
+            }
         }
 
 
@@ -1224,7 +1285,13 @@ int main(int argc, char *argv[])
 
 
             CompactListList<label> regionCellOrder =
-                regionRenumber(renumberPtr(), mesh, cellToRegion);
+                regionRenumber
+                (
+                    renumberPtr(),
+                    mesh,
+                    cellToRegion,
+                    decomposePtr().nDomains()
+                );
 
             cellOrder = regionCellOrder.values();
 
@@ -1586,19 +1653,25 @@ int main(int argc, char *argv[])
             );
             reduce(band, maxOp<label>());
             reduce(profile, sumOp<scalar>());
-            reduce(sumSqrIntersect, sumOp<scalar>());
 
-            scalar rmsFrontwidth = Foam::sqrt
-            (
-                sumSqrIntersect/mesh.globalData().nTotalCells()
-            );
+            Info<< "After renumbering";
+            if (doDecompose)
+            {
+                Info<< " [values are misleading with -decompose option]";
+            }
 
-            Info<< "After renumbering :" << nl
+            Info<< nl
                 << "    band           : " << band << nl
                 << "    profile        : " << profile << nl;
 
             if (doFrontWidth)
             {
+                reduce(sumSqrIntersect, sumOp<scalar>());
+                scalar rmsFrontwidth = Foam::sqrt
+                (
+                    sumSqrIntersect/mesh.globalData().nTotalCells()
+                );
+
                 Info<< "    rms frontwidth : " << rmsFrontwidth << nl;
             }
 
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/subsetAdjacency.H b/applications/utilities/mesh/manipulation/renumberMesh/subsetAdjacency.H
new file mode 100644
index 0000000000..0c098616cf
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/renumberMesh/subsetAdjacency.H
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+   This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+    Subsetting of an adjacency matrix (as CompactListList).
+    Can be relocated elsewhere.
+
+\*---------------------------------------------------------------------------*/
+
+#include "CompactListList.H"
+#include "bitSet.H"
+#include "ListOps.H"
+#include "Map.H"
+
+namespace Foam
+{
+
+// Perform a subset of the adjacency matrix
+CompactListList<label> subsetAdjacency
+(
+    const bitSet& select,   // could also be labelHashSet
+    const CompactListList<label>& input,
+    labelList& subMap
+)
+{
+    // Corresponds to cellMap etc (the original selection)
+    subMap = select.sortedToc();
+
+    // Ensure that the subMap corresponds to a valid subset
+    {
+        label validSize = 0;
+
+        const label nTotal = input.size();
+
+        forAllReverse(subMap, i)
+        {
+            if (subMap[i] < nTotal)
+            {
+                validSize = i + 1;
+                break;
+            }
+        }
+
+        subMap.resize(validSize);
+    }
+
+
+    // Assumed to be sparse - use Map for reverse lookup
+    const Map<label> reverseMap(invertToMap(subMap));
+
+
+    // Pass 1: determine the selected sub-sizes
+    labelList sizes(subMap.size(), Foam::zero{});
+
+    forAll(subMap, idx)
+    {
+        for (const label nbr : input[subMap[idx]])
+        {
+            if
+            (
+                select.test(nbr)
+             && reverseMap.contains(nbr)  // extra consistency (paranoid)
+            )
+            {
+                ++sizes[idx];
+            }
+        }
+    }
+
+
+    CompactListList<label> output(sizes);
+
+    // Reuse sizes as output offset into output.values()
+    sizes = labelList::subList(output.offsets(), output.size());
+    labelList& values = output.values();
+
+
+    // Pass 2: extract sub-adjacent matrix
+
+    label newNbr = -1;
+
+    forAll(subMap, idx)
+    {
+        for (const label nbr : input[subMap[idx]])
+        {
+            if
+            (
+                select.test(nbr)
+             && (newNbr = reverseMap.lookup(nbr, -1)) >= 0
+            )
+            {
+                values[sizes[idx]++] = newNbr;
+            }
+        }
+    }
+
+    return output;
+}
+
+
+// Perform a subset of the adjacency matrix
+CompactListList<label> subsetAdjacency
+(
+    const labelRange& slice,
+    const CompactListList<label>& input,
+    labelList& subMap
+)
+{
+    // Ensure that the selection corresponds to a valid subset
+    const labelRange select = slice.subset0(input.size());
+
+    // Corresponds to cellMap etc (the original selection)
+    subMap = Foam::identity(select);
+
+
+    // Pass 1: determine the selected sub-sizes
+    labelList sizes(subMap.size(), Foam::zero{});
+
+    forAll(subMap, idx)
+    {
+        for (const label nbr : input[subMap[idx]])
+        {
+            if (select.contains(nbr))
+            {
+                ++sizes[idx];
+            }
+        }
+    }
+
+
+    CompactListList<label> output(sizes);
+
+    // Reuse sizes as output offset into output.values()
+    sizes = labelList::subList(output.offsets(), output.size());
+    labelList& values = output.values();
+
+
+    // Pass 2: extract sub-adjacent matrix
+
+    const label localOffset = select.start();
+
+    forAll(subMap, idx)
+    {
+        for (const label nbr : input[subMap[idx]])
+        {
+            if (select.contains(nbr))
+            {
+                values[sizes[idx]++] = nbr - localOffset;
+            }
+        }
+    }
+
+    return output;
+}
+
+} // End namespace Foam
+
+
+// ************************************************************************* //