From 61aaacd088fc9fcf3afb40526f1c943580fdcbf7 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Mon, 4 Mar 2024 16:48:59 +0100
Subject: [PATCH] ENH: adjust renumbering methods, extend renumberMesh options

- renumberMesh now has -dry-run, -write-maps, -no-fields,
  -renumber-method, -renumber-coeffs options.

  * Use -dry-run with -write-maps to visualize the before/after
    effects of renumbering (creates a VTK file).

  * -no-fields to renumber the mesh only.
    This is useful and faster when the input fields are uniform
    and the -overwrite option is specified.

  * -renumber-method allows a quick means of specifying a different
    default renumber method (instead of Cuthill-McKee).

    The -renumber-coeffs option allows passing of dictionary content
    for the method.

    Examples,

       // Different ways to specify reverse Cuthill-McKee

       *  -renumber-method RCM
       *  -renumber-coeffs 'reverse true;'
       *  -renumber-method CuthillMcKee
       *  -renumber-coeffs 'reverse true;'
       *  -renumber-coeffs 'method CuthillMcKee; reverse true;'

       // Other (without dictionary coefficients)
       *  renumberMesh -renumber-method random

       // Other (with dictionary coefficients)
       renumberMesh \
           -renumber-method spring \
           -renumber-coeffs 'maxCo 0.1; maxIter 1000; freezeFraction 0.99;'

       // Other (with additional libraries)
       renumberMesh -renumber-method zoltan -lib zoltanRenumber

COMP: build zoltan renumbering to MPI-specific location

- zoltan and Sloan renumbering are now longer automatically linked to
  the renumberMesh utility but must be separately loaded by a
  command-line option or through a dictionary "libs" entry.

ENH: add output cellID for decomposePar -dry-run -cellDist
---
 .../cavity/system/renumberMeshDict-random     |  112 --
 .../mesh/manipulation/renumberMesh/Allwmake   |   26 -
 .../manipulation/renumberMesh/Make/options    |   11 +-
 .../manipulation/renumberMesh/renumberMesh.C  | 1078 +++++++++++------
 .../domainDecompositionDryRunWrite.C          |    1 +
 bin/tools/CleanFunctions                      |    2 +-
 etc/caseDicts/annotated/renumberMeshDict      |   45 +-
 etc/config.sh/zoltan                          |   34 +
 src/Allwmake                                  |    2 +-
 src/renumber/Allwmake                         |    1 +
 src/renumber/SloanRenumber/SloanRenumber.C    |   12 +-
 src/renumber/SloanRenumber/SloanRenumber.H    |   22 +-
 .../CuthillMcKeeRenumber.C                    |   68 +-
 .../CuthillMcKeeRenumber.H                    |   72 +-
 src/renumber/renumberMethods/Make/files       |    5 +-
 .../manualRenumber/manualRenumber.C           |   16 +-
 .../manualRenumber/manualRenumber.H           |   20 +-
 .../renumberMethods/noRenumber/noRenumber.C   |  102 ++
 .../renumberMethods/noRenumber/noRenumber.H   |  124 ++
 .../randomRenumber/randomRenumber.C           |   65 +-
 .../randomRenumber/randomRenumber.H           |   42 +-
 .../renumberMethod/renumberMethod.C           |   27 +-
 .../renumberMethod/renumberMethod.H           |   91 +-
 .../springRenumber/springRenumber.C           |   48 +-
 .../springRenumber/springRenumber.H           |    3 +
 .../structuredRenumber/structuredRenumber.C   |   12 +-
 .../structuredRenumber/structuredRenumber.H   |   15 +-
 src/renumber/zoltanRenumber/Make/files        |    3 +-
 src/renumber/zoltanRenumber/Make/options      |    6 +-
 src/renumber/zoltanRenumber/zoltanRenumber.C  |  344 ++++--
 src/renumber/zoltanRenumber/zoltanRenumber.H  |    6 +-
 tutorials/mesh/parallel/cavity/Allrun         |    3 +-
 .../cavity/system/renumberMeshDict-random     |  109 --
 .../geometric/Allclean                        |    1 -
 wmake/scripts/have_zoltan                     |   31 +-
 35 files changed, 1626 insertions(+), 933 deletions(-)
 delete mode 100644 applications/test/mapDistributePolyMesh/cavity/system/renumberMeshDict-random
 delete mode 100755 applications/utilities/mesh/manipulation/renumberMesh/Allwmake
 create mode 100644 etc/config.sh/zoltan
 create mode 100644 src/renumber/renumberMethods/noRenumber/noRenumber.C
 create mode 100644 src/renumber/renumberMethods/noRenumber/noRenumber.H
 delete mode 100644 tutorials/mesh/parallel/cavity/system/renumberMeshDict-random

diff --git a/applications/test/mapDistributePolyMesh/cavity/system/renumberMeshDict-random b/applications/test/mapDistributePolyMesh/cavity/system/renumberMeshDict-random
deleted file mode 100644
index d3d55f4fd1..0000000000
--- a/applications/test/mapDistributePolyMesh/cavity/system/renumberMeshDict-random
+++ /dev/null
@@ -1,112 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2312                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    note        "mesh renumbering dictionary";
-    object      renumberMeshDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Write maps from renumbered back to original mesh
-writeMaps true;
-
-// Optional entry: sort cells on coupled boundaries to last for use with
-// e.g. nonBlockingGaussSeidel.
-sortCoupledFaceCells false;
-
-// Optional entry: renumber on a block-by-block basis. It uses a
-// blockCoeffs dictionary to construct a decompositionMethod to do
-// a block subdivision) and then applies the renumberMethod to each
-// block in turn. This can be used in large cases to keep the blocks
-// fitting in cache with all the cache misses bunched at the end.
-// This number is the approximate size of the blocks - this gets converted
-// to a number of blocks that is the input to the decomposition method.
-//blockSize 1000;
-
-// Optional entry: sort points into internal and boundary points
-//orderPoints false;
-
-// Optional: suppress renumbering cellSets,faceSets,pointSets
-//renumberSets false;
-
-
-//method          CuthillMcKee;
-//method          Sloan;
-//method          manual;
-method          random;
-//method          structured;
-//method          spring;
-//method          zoltan;             // only if compiled with zoltan support
-
-//CuthillMcKeeCoeffs
-//{
-//    // Reverse CuthillMcKee (RCM) or plain
-//    reverse true;
-//}
-
-manualCoeffs
-{
-    // In system directory: new-to-original (i.e. order) labelIOList
-    dataFile "cellMap";
-}
-
-
-// For extruded (i.e. structured in one direction) meshes
-structuredCoeffs
-{
-    // Patches that mesh was extruded from. These determine the starting
-    // layer of cells
-    patches (movingWall);
-    // Method to renumber the starting layer of cells
-    method  random;
-
-    // Renumber in columns (depthFirst) or in layers
-    depthFirst true;
-
-    // Reverse ordering
-    reverse false;
-}
-
-
-springCoeffs
-{
-    // Maximum jump of cell indices. Is fraction of number of cells
-    maxCo 0.01;
-
-    // Limit the amount of movement; the fraction maxCo gets decreased
-    // with every iteration
-    freezeFraction 0.999;
-
-    // Maximum number of iterations
-    maxIter 1000;
-}
-
-
-blockCoeffs
-{
-    method          scotch;
-    //method hierarchical;
-    //hierarchicalCoeffs
-    //{
-    //    n           (1 2 1);
-    //    delta       0.001;
-    //    order       xyz;
-    //}
-}
-
-
-zoltanCoeffs
-{
-    ORDER_METHOD    LOCAL_HSFC;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/Allwmake b/applications/utilities/mesh/manipulation/renumberMesh/Allwmake
deleted file mode 100755
index 0f0d32738e..0000000000
--- a/applications/utilities/mesh/manipulation/renumberMesh/Allwmake
+++ /dev/null
@@ -1,26 +0,0 @@
-#!/bin/sh
-cd "${0%/*}" || exit                                # Run from this directory
-. ${WM_PROJECT_DIR:?}/wmake/scripts/AllwmakeParseArguments # (error catching)
-. ${WM_PROJECT_DIR:?}/wmake/scripts/sysFunctions    # General system functions
-. ${WM_PROJECT_DIR:?}/wmake/scripts/have_zoltan
-
-#------------------------------------------------------------------------------
-
-unset COMP_FLAGS LINK_FLAGS
-
-if findLibrary "$FOAM_LIBBIN/libSloanRenumber" > /dev/null
-then
-    echo "    found libSloanRenumber -- enabling sloan renumbering support."
-    export LINK_FLAGS="$LINK_FLAGS -lSloanRenumber"
-fi
-
-if findLibrary "$FOAM_LIBBIN/libzoltanRenumber" > /dev/null && have_zoltan
-then
-    echo "    found libzoltanRenumber -- enabling zoltan renumbering support."
-    export COMP_FLAGS="$COMP_FLAGS -DHAVE_ZOLTAN"
-    export LINK_FLAGS="$LINK_FLAGS -lzoltanRenumber -L$ZOLTAN_LIB_DIR -lzoltan"
-fi
-
-wmake $targetType
-
-#------------------------------------------------------------------------------
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/Make/options b/applications/utilities/mesh/manipulation/renumberMesh/Make/options
index 004f82049f..8aca81f2dc 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/Make/options
+++ b/applications/utilities/mesh/manipulation/renumberMesh/Make/options
@@ -1,21 +1,22 @@
 EXE_INC = \
-    $(COMP_FLAGS) \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/renumber/renumberMethods/lnInclude \
-    -I$(LIB_SRC)/renumber/zoltanRenumber/lnInclude \
+    -I$(LIB_SRC)/parallel/decompose/decompose/lnInclude \
     -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/parallel/reconstruct/reconstruct/lnInclude
 
 EXE_LIBS = \
-    -lfiniteVolume \
+    -lfileFormats \
     -lmeshTools \
     -ldynamicMesh \
+    -lfiniteVolume \
     -lgenericPatchFields \
     -lrenumberMethods \
+    -ldecompose \
     -lreconstruct \
-    $(LINK_FLAGS) \
     -ldecompositionMethods \
     -L$(FOAM_LIBBIN)/dummy \
     -lkahipDecomp -lmetisDecomp -lscotchDecomp
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
index b35202a4ad..7c08c679a5 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,12 +34,96 @@ Description
     Renumbers the cell list in order to reduce the bandwidth, reading and
     renumbering all fields from all the time directories.
 
-    By default uses bandCompression (CuthillMcKee) but will
-    read system/renumberMeshDict if -dict option is present
+    By default uses bandCompression (Cuthill-McKee) or the method
+    specified by the -renumber-method option, but will read
+    system/renumberMeshDict if -dict option is present
+
+Usage
+    \b renumberMesh [OPTIONS]
+
+    Options:
+      - \par -allRegions
+        Use all regions in regionProperties
+
+      - \par -case \<dir\>
+        Specify case directory to use (instead of the cwd).
+
+      - \par -constant
+        Include the 'constant/' dir in the times list.
+
+      - \par -decompose
+        Aggregate initially with a decomposition method (serial only)
+
+      - \par -decomposeParDict \<file\>
+        Use specified file for decomposePar dictionary.
+
+      - \par -dict \<file\>
+        Use specified file for renumberMeshDict dictionary.
+
+      - \par -dry-run
+        Test only
+
+      - \par -frontWidth
+        Calculate the rms of the front-width
+
+      - \par -latestTime
+        Select the latest time.
+
+      - \par -lib \<name\>
+        Additional library or library list to load (can be used multiple times).
+
+      - \par -no-fields
+        Suppress renumber of fields
+
+      - \par -noZero
+        Exclude the \a 0 dir from the times list.
+
+      - \par -overwrite
+        Overwrite existing mesh/results files
+
+      - \par -parallel
+        Run in parallel
+
+      - \par -region \<regionName\>
+        Renumber named region.
+
+      - \par -regions \<wordRes\>
+        Renumber named regions.
+
+      - \par -renumber-coeffs \<string-content\>
+        String to create renumber dictionary contents.
+
+      - \par -renumber-method \<name\>
+        Specify renumber method (default: CuthillMcKee) without dictionary
+
+      - \par -time \<value\>
+        Specify time to select
+
+      - \par -verbose
+        Additional verbosity.
+
+      - \par -doc
+        Display documentation in browser.
+
+      - \par -doc-source
+        Display source code in browser.
+
+      - \par -help
+        Display short help and exit.
+
+      - \par -help-man
+        Display full help (manpage format) and exit.
+
+      - \par -help-notes
+        Display help notes (description) and exit.
+
+      - \par -help-full
+        Display full help and exit.
 
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
+#include "timeSelector.H"
 #include "IOobjectList.H"
 #include "fvMesh.H"
 #include "polyTopoChange.H"
@@ -48,8 +132,9 @@ Description
 #include "surfaceFields.H"
 #include "SortableList.H"
 #include "decompositionMethod.H"
+#include "decompositionModel.H"
 #include "renumberMethod.H"
-#include "zeroGradientFvPatchFields.H"
+#include "foamVtkInternalMeshWriter.H"
 #include "CuthillMcKeeRenumber.H"
 #include "fvMeshSubset.H"
 #include "cellSet.H"
@@ -60,11 +145,6 @@ Description
 #include "regionProperties.H"
 #include "polyMeshTools.H"
 
-#ifdef HAVE_ZOLTAN
-    #include "zoltanRenumber.H"
-#endif
-
-
 using namespace Foam;
 
 
@@ -80,16 +160,17 @@ tmp<volScalarField> createScalarField
     (
         name,
         mesh,
-        dimensionedScalar(dimless, Zero),
+        dimensionedScalar(word::null, dimless, -1),
         fvPatchFieldBase::zeroGradientType()
     );
     auto& fld = tfld.ref();
 
-    forAll(fld, celli)
+    forAll(elems, celli)
     {
        fld[celli] = elems[celli];
     }
 
+    fld.correctBoundaryConditions();
     return tfld;
 }
 
@@ -140,10 +221,10 @@ void getBand
     bandwidth = max(cellBandwidth);
 
     // Do not use field algebra because of conversion label to scalar
-    profile = 0.0;
-    forAll(cellBandwidth, celli)
+    profile = 0;
+    for (const label width : cellBandwidth)
     {
-        profile += 1.0*cellBandwidth[celli];
+        profile += scalar(width);
     }
 
     sumSqrIntersect = 0.0;
@@ -175,54 +256,46 @@ labelList getFaceOrder
 
     label newFacei = 0;
 
-    labelList nbr;
-    labelList order;
+    DynamicList<label> nbr(64);
+    DynamicList<label> order(64);
 
     forAll(cellOrder, newCelli)
     {
-        label oldCelli = cellOrder[newCelli];
+        const label oldCelli = cellOrder[newCelli];
 
         const cell& cFaces = mesh.cells()[oldCelli];
 
         // Neighbouring cells
-        nbr.setSize(cFaces.size());
+        nbr.clear();
 
-        forAll(cFaces, i)
+        for (const label facei : cFaces)
         {
-            label facei = cFaces[i];
+            label nbrCelli = -1;
 
             if (mesh.isInternalFace(facei))
             {
                 // Internal face. Get cell on other side.
-                label nbrCelli = reverseCellOrder[mesh.faceNeighbour()[facei]];
+                nbrCelli = reverseCellOrder[mesh.faceNeighbour()[facei]];
                 if (nbrCelli == newCelli)
                 {
                     nbrCelli = reverseCellOrder[mesh.faceOwner()[facei]];
                 }
 
-                if (newCelli < nbrCelli)
+                // The nbrCell is actually the master (let it handle the face)
+                if (nbrCelli <= newCelli)
                 {
-                    // Celli is master
-                    nbr[i] = nbrCelli;
-                }
-                else
-                {
-                    // nbrCell is master. Let it handle this face.
-                    nbr[i] = -1;
+                    nbrCelli = -1;
                 }
             }
-            else
-            {
-                // External face. Do later.
-                nbr[i] = -1;
-            }
+
+            nbr.push_back(nbrCelli);
         }
 
-        sortedOrder(nbr, order);
+        Foam::sortedOrder(nbr, order);
 
         for (const label index : order)
         {
-            if (nbr[index] != -1)
+            if (nbr[index] >= 0)
             {
                 oldToNewFace[cFaces[index]] = newFacei++;
             }
@@ -266,73 +339,71 @@ labelList getRegionFaceOrder
 
     label newFacei = 0;
 
-    label prevRegion = -1;
+    DynamicList<label> nbr(64);
+    DynamicList<label> order(64);
 
     forAll(cellOrder, newCelli)
     {
-        label oldCelli = cellOrder[newCelli];
-
-        if (cellToRegion[oldCelli] != prevRegion)
-        {
-            prevRegion = cellToRegion[oldCelli];
-        }
+        const label oldCelli = cellOrder[newCelli];
 
         const cell& cFaces = mesh.cells()[oldCelli];
 
-        SortableList<label> nbr(cFaces.size());
+        // Neighbouring cells
+        nbr.clear();
 
-        forAll(cFaces, i)
+        for (const label facei : cFaces)
         {
-            label facei = cFaces[i];
+            label nbrCelli = -1;
 
             if (mesh.isInternalFace(facei))
             {
                 // Internal face. Get cell on other side.
-                label nbrCelli = reverseCellOrder[mesh.faceNeighbour()[facei]];
+                nbrCelli = reverseCellOrder[mesh.faceNeighbour()[facei]];
                 if (nbrCelli == newCelli)
                 {
                     nbrCelli = reverseCellOrder[mesh.faceOwner()[facei]];
                 }
 
-                if (cellToRegion[oldCelli] != cellToRegion[cellOrder[nbrCelli]])
+                // The nbrCell is actually the master (let it handle the face)
+                if (nbrCelli <= newCelli)
                 {
-                    // Treat like external face. Do later.
-                    nbr[i] = -1;
-                }
-                else if (newCelli < nbrCelli)
-                {
-                    // Celli is master
-                    nbr[i] = nbrCelli;
+                    nbrCelli = -1;
                 }
                 else
                 {
-                    // nbrCell is master. Let it handle this face.
-                    nbr[i] = -1;
+                    // A region boundary? - treat like an external face
+                    label ownRegion = cellToRegion[oldCelli];
+                    label neiRegion = cellToRegion[cellOrder[nbrCelli]];
+
+                    if (ownRegion != neiRegion)
+                    {
+                        nbrCelli = -1;
+                    }
                 }
             }
-            else
-            {
-                // External face. Do later.
-                nbr[i] = -1;
-            }
+
+            nbr.push_back(nbrCelli);
         }
 
-        nbr.sort();
+        Foam::sortedOrder(nbr, order);
 
-        forAll(nbr, i)
+        for (const label index : order)
         {
-            if (nbr[i] != -1)
+            if (nbr[index] >= 0)
             {
-                oldToNewFace[cFaces[nbr.indices()[i]]] = newFacei++;
+                oldToNewFace[cFaces[index]] = newFacei++;
             }
         }
     }
 
+    // This seems to be broken...
+
     // Do region interfaces
-    label nRegions = max(cellToRegion)+1;
     {
+        const label nRegions = max(cellToRegion)+1;
+
         // Sort in increasing region
-        SortableList<label> sortKey(mesh.nFaces(), labelMax);
+        SortableList<label> sortKey(mesh.nInternalFaces(), labelMax);
 
         for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
         {
@@ -346,10 +417,10 @@ labelList getRegionFaceOrder
                    +max(ownRegion, neiRegion);
             }
         }
+
         sortKey.sort();
 
         // Extract.
-        label prevKey = -1;
         forAll(sortKey, i)
         {
             label key = sortKey[i];
@@ -359,11 +430,6 @@ labelList getRegionFaceOrder
                 break;
             }
 
-            if (prevKey != key)
-            {
-                prevKey = key;
-            }
-
             oldToNewFace[sortKey.indices()[i]] = newFacei++;
         }
     }
@@ -381,8 +447,7 @@ labelList getRegionFaceOrder
         if (oldToNewFace[facei] == -1)
         {
             FatalErrorInFunction
-                << "Did not determine new position"
-                << " for face " << facei
+                << "Did not determine new position for face " << facei
                 << abort(FatalError);
         }
     }
@@ -426,13 +491,10 @@ autoPtr<mapPolyMesh> reorderMesh
     labelHashSet flipFaceFlux(newOwner.size());
     forAll(newNeighbour, facei)
     {
-        label own = newOwner[facei];
-        label nei = newNeighbour[facei];
-
-        if (nei < own)
+        if (newOwner[facei] > newNeighbour[facei])
         {
-            newFaces[facei].flip();
             std::swap(newOwner[facei], newNeighbour[facei]);
+            newFaces[facei].flip();
             flipFaceFlux.insert(facei);
         }
     }
@@ -543,55 +605,94 @@ autoPtr<mapPolyMesh> reorderMesh
 }
 
 
-// Return new to old cell numbering
-labelList regionRenumber
+// Return new to old cell numbering, region-wise
+CompactListList<label> regionRenumber
 (
     const renumberMethod& method,
     const fvMesh& mesh,
-    const labelList& cellToRegion
+    const labelList& cellToRegion,
+    label nRegions = -1   // Number of regions or auto-detect
 )
 {
-    Info<< "Determining cell order:" << endl;
+    // Info<< "Determining cell order:" << endl;
 
-    labelList cellOrder(cellToRegion.size());
-
-    label nRegions = max(cellToRegion)+1;
-
-    labelListList regionToCells(invertOneToMany(nRegions, cellToRegion));
-
-    label celli = 0;
-
-    forAll(regionToCells, regioni)
+    if (nRegions < 0)
     {
-        Info<< "    region " << regioni << " starts at " << celli << endl;
+        nRegions = Foam::max(cellToRegion)+1;
+    }
 
-        // Make sure no parallel comms
-        const bool oldParRun = UPstream::parRun(false);
+    // Initially the per-region cell selection
+    CompactListList<label> regionCellOrder
+    (
+        invertOneToManyCompact(nRegions, cellToRegion)
+    );
 
-        // Per region do a reordering.
-        fvMeshSubset subsetter(mesh, regioni, cellToRegion);
-
-        const fvMesh& subMesh = subsetter.subMesh();
-
-        labelList subCellOrder = method.renumber
-        (
-            subMesh,
-            subMesh.cellCentres()
-        );
-
-        UPstream::parRun(oldParRun);  // Restore parallel state
-
-
-        const labelList& cellMap = subsetter.cellMap();
-
-        forAll(subCellOrder, i)
+    if (method.needs_mesh())
+    {
+        forAll(regionCellOrder, regioni)
         {
-            cellOrder[celli++] = cellMap[subCellOrder[i]];
+            // Info<< "    region " << regioni
+            //     << " starts at " << regionCellOrder.localStart(regioni)
+            //     << nl;
+
+            // No parallel communication
+            const bool oldParRun = UPstream::parRun(false);
+
+            // Get local sub-mesh
+            fvMeshSubset subsetter(mesh, regionCellOrder[regioni]);
+
+            // Note: cellMap will be identical to regionToCells[regioni]
+            // (assuming they are properly sorted!)
+            const labelList& cellMap = subsetter.cellMap();
+
+            labelList subCellOrder = method.renumber
+            (
+                subsetter.subMesh(),
+                subsetter.subMesh().cellCentres()
+            );
+
+            UPstream::parRun(oldParRun);  // Restore parallel state
+
+            // Per region reordering
+            regionCellOrder[regioni] = labelUIndList(cellMap, subCellOrder);
         }
     }
-    Info<< endl;
+    else
+    {
+        forAll(regionCellOrder, regioni)
+        {
+            // Info<< "    region " << regioni
+            //     << " starts at " << regionCellOrder.localStart(regioni)
+            //     << nl;
 
-    return cellOrder;
+            // 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
+            );
+
+            // Note: cellCentres not needed by every renumber method
+            labelList subCellOrder = method.renumber
+            (
+                cellCells,
+                pointField(mesh.cellCentres(), cellMap)
+            );
+
+            UPstream::parRun(oldParRun);  // Restore parallel state
+
+            // Per region reordering
+            regionCellOrder[regioni] = labelUIndList(cellMap, subCellOrder);
+        }
+    }
+    // Info<< endl;
+
+    return regionCellOrder;
 }
 
 
@@ -601,48 +702,138 @@ int main(int argc, char *argv[])
 {
     argList::addNote
     (
-        "Renumber mesh cells to reduce the bandwidth"
+        "Renumber mesh cells to reduce the bandwidth. Use the -lib option or"
+        " dictionary 'libs' entry to load additional libraries"
     );
 
-    #include "addAllRegionOptions.H"
-    #include "addOverwriteOption.H"
-    #include "addTimeOptions.H"
+    argList::addDryRunOption
+    (
+        "Test without writing. "
+        "Changes -write-maps to write VTK output."
+    );
+    argList::addVerboseOption();
 
     argList::addOption("dict", "file", "Alternative renumberMeshDict");
 
     argList::addBoolOption
     (
         "frontWidth",
-        "Calculate the rms of the front-width"
+        "Calculate the RMS of the front-width"
+    );
+
+    argList::addBoolOption
+    (
+        "decompose",
+        "Aggregate initially with a decomposition method (serial only)"
+    );
+
+    argList::addBoolOption
+    (
+        "write-maps",
+        "Write renumber mappings"
+    );
+
+    argList::addBoolOption
+    (
+        "no-fields",
+        "Suppress renumbering of fields (eg, when they are only uniform)"
+    );
+
+    argList::addBoolOption
+    (
+        "list-renumber",
+        "List available renumbering methods"
+    );
+
+    argList::addOption
+    (
+        "renumber-method",
+        "name",
+        "Specify renumber method (default: CuthillMcKee) without dictionary"
+    );
+
+    argList::addOption
+    (
+        "renumber-coeffs",
+        "string-content",
+        "Specify renumber coefficients (dictionary content) as string. "
+        "eg, 'reverse true;'"
     );
 
     argList::noFunctionObjects();  // Never use function objects
 
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "getAllRegionOptions.H"
+    #include "addAllRegionOptions.H"
+    #include "addOverwriteOption.H"
+    #include "addTimeOptions.H"
 
-    // Force linker to include zoltan symbols. This section is only needed since
-    // Zoltan is a static library
-    #ifdef HAVE_ZOLTAN
-        Info<< "renumberMesh built with zoltan support." << nl << endl;
-        (void)zoltanRenumber::typeName;
-    #endif
+    // -------------------------
+
+    #include "setRootCase.H"
+
+    {
+        bool listOptions = false;
+
+        if (args.found("list-renumber"))
+        {
+            listOptions = true;
+            Info<< nl
+                << "Available renumber methods:" << nl
+                << "    "
+                << flatOutput(renumberMethod::supportedMethods()) << nl
+                << nl;
+        }
+
+        if (listOptions)
+        {
+            return 0;
+        }
+    }
+
+    const bool dryrun = args.dryRun();
+
+    const bool readDict = args.found("dict");
+    const bool doFrontWidth = args.found("frontWidth");
+    const bool overwrite = args.found("overwrite");
+    const bool doFields = !args.found("no-fields");
+    const bool doDecompose = args.found("decompose");
+
+    word renumberMethodName;
+    args.readIfPresent("renumber-method", renumberMethodName);
+
+    if (doDecompose && UPstream::parRun())
+    {
+        FatalErrorIn(args.executable())
+            << "Cannot use -decompose option in parallel"
+            << " ... giving up" << nl
+            << exit(FatalError);
+    }
+
+    #include "createTime.H"
+
+    // Allow override of decomposeParDict location
+    fileName decompDictFile;
+    if
+    (
+        args.readIfPresent("decomposeParDict", decompDictFile)
+     && !decompDictFile.empty() && !decompDictFile.isAbsolute()
+    )
+    {
+        decompDictFile = runTime.globalPath()/decompDictFile;
+    }
+
+    // Get region names
+    #include "getAllRegionOptions.H"
 
     // Get times list
     instantList Times = runTime.times();
 
-    // Set startTime and endTime depending on -time and -latestTime options
+    // Set startTime depending on -time and -latestTime options
     #include "checkTimeOptions.H"
 
     runTime.setTime(Times[startTime], startTime);
 
     #include "createNamedMeshes.H"
 
-    const bool readDict = args.found("dict");
-    const bool doFrontWidth = args.found("frontWidth");
-    const bool overwrite = args.found("overwrite");
-
 
     for (fvMesh& mesh : meshes)
     {
@@ -667,13 +858,11 @@ int main(int argc, char *argv[])
 
         reduce(band, maxOp<label>());
         reduce(profile, sumOp<scalar>());
+        reduce(sumSqrIntersect, sumOp<scalar>());
+
         scalar rmsFrontwidth = Foam::sqrt
         (
-            returnReduce
-            (
-                sumSqrIntersect,
-                sumOp<scalar>()
-            )/mesh.globalData().nTotalCells()
+            sumSqrIntersect/mesh.globalData().nTotalCells()
         );
 
         Info<< "Mesh " << mesh.name()
@@ -690,12 +879,12 @@ int main(int argc, char *argv[])
         Info<< endl;
 
         bool sortCoupledFaceCells = false;
-        bool writeMaps = false;
+        bool writeMaps = args.found("write-maps");
         bool orderPoints = false;
         label blockSize = 0;
 
         // Construct renumberMethod
-        autoPtr<IOdictionary> renumberDictPtr;
+        dictionary renumberDict;
         autoPtr<renumberMethod> renumberPtr;
 
         if (readDict)
@@ -705,38 +894,31 @@ int main(int argc, char *argv[])
 
             Info<< "Renumber according to " << dictIO.name() << nl << endl;
 
-            renumberDictPtr.reset(new IOdictionary(dictIO));
-            const IOdictionary& renumberDict = renumberDictPtr();
+            renumberDict = IOdictionary::readContents(dictIO);
 
             renumberPtr = renumberMethod::New(renumberDict);
 
-            sortCoupledFaceCells = renumberDict.getOrDefault
-            (
-                "sortCoupledFaceCells",
-                false
-            );
-            if (sortCoupledFaceCells)
+            // This options are likely orthogonal to -decompose mode
+            if (!doDecompose)
             {
-                Info<< "Sorting cells on coupled boundaries to be last." << nl
-                    << endl;
-            }
+                sortCoupledFaceCells =
+                    renumberDict.getOrDefault("sortCoupledFaceCells", false);
 
-            blockSize = renumberDict.getOrDefault("blockSize", 0);
-            if (blockSize > 0)
-            {
-                Info<< "Ordering cells into regions of size " << blockSize
-                    << " (using decomposition);"
-                    << " ordering faces into region-internal"
-                    << " and region-external."
-                    << nl << endl;
-
-                if (blockSize < 0 || blockSize >= mesh.nCells())
+                if (sortCoupledFaceCells)
                 {
-                    FatalErrorInFunction
-                        << "Block size " << blockSize
-                        << " should be positive integer"
-                        << " and less than the number of cells in the mesh."
-                        << exit(FatalError);
+                    Info<< "Sorting cells on coupled boundaries to be last."
+                        << nl << endl;
+                }
+
+                blockSize = renumberDict.getOrDefault<label>("blockSize", 0);
+
+                if (blockSize > 0)
+                {
+                    Info<< "Ordering cells into regions of size " << blockSize
+                        << " (using decomposition);"
+                        << " ordering faces into region-internal"
+                        << " and region-external."
+                        << nl << endl;
                 }
             }
 
@@ -744,27 +926,56 @@ int main(int argc, char *argv[])
             if (orderPoints)
             {
                 Info<< "Ordering points into internal and boundary points."
-                    << nl
-                    << endl;
+                    << nl << endl;
             }
 
-            renumberDict.readEntry("writeMaps", writeMaps);
-            if (writeMaps)
+            if
+            (
+                renumberDict.readIfPresent("writeMaps", writeMaps)
+             && writeMaps
+            )
             {
-                Info<< "Writing renumber maps (new to old) to polyMesh." << nl
-                    << endl;
+                Info<< "Writing renumber maps (new to old) to polyMesh."
+                    << nl << endl;
             }
         }
         else
         {
-            Info<< "Using default renumberMethod." << nl << endl;
-            dictionary renumberDict;
-            renumberPtr.reset(new CuthillMcKeeRenumber(renumberDict));
+            if (args.found("renumber-coeffs"))
+            {
+                ITstream is = args.lookup("renumber-coeffs");
+
+                is >> renumberDict;
+
+                Info<< "Specified renumber coefficients:" << nl
+                    << renumberDict << nl;
+            }
+
+            if (!renumberMethodName.empty())
+            {
+                // Specify/override the "method" within dictionary
+                renumberDict.set("method", renumberMethodName);
+                renumberPtr = renumberMethod::New(renumberDict);
+            }
+            else if (renumberDict.found("method"))
+            {
+                // Use the "method" type within dictionary
+                renumberPtr = renumberMethod::New(renumberDict);
+            }
         }
 
-        Info<< "Selecting renumberMethod " << renumberPtr().type() << nl
-            << endl;
-
+        // Default renumbering method
+        if (!renumberPtr)
+        {
+            renumberPtr.reset(new CuthillMcKeeRenumber(renumberDict));
+            Info<< "Using renumber-method: " << renumberPtr().type()
+                << " [default]" << nl << endl;
+        }
+        else
+        {
+            Info<< "Using renumber-method: " << renumberPtr().type()
+                << nl << endl;
+        }
 
 
         // Read parallel reconstruct maps
@@ -776,10 +987,10 @@ int main(int argc, char *argv[])
                 mesh.facesInstance(),
                 polyMesh::meshSubDir,
                 mesh,
-                IOobject::READ_IF_PRESENT,
+                (dryrun ? IOobject::NO_READ : IOobject::READ_IF_PRESENT),
                 IOobject::AUTO_WRITE
             ),
-            labelList(0)
+            labelList()
         );
 
         labelIOList faceProcAddressing
@@ -790,10 +1001,10 @@ int main(int argc, char *argv[])
                 mesh.facesInstance(),
                 polyMesh::meshSubDir,
                 mesh,
-                IOobject::READ_IF_PRESENT,
+                (dryrun ? IOobject::NO_READ : IOobject::READ_IF_PRESENT),
                 IOobject::AUTO_WRITE
             ),
-            labelList(0)
+            labelList()
         );
         labelIOList pointProcAddressing
         (
@@ -803,10 +1014,10 @@ int main(int argc, char *argv[])
                 mesh.pointsInstance(),
                 polyMesh::meshSubDir,
                 mesh,
-                IOobject::READ_IF_PRESENT,
+                (dryrun ? IOobject::NO_READ : IOobject::READ_IF_PRESENT),
                 IOobject::AUTO_WRITE
             ),
-            labelList(0)
+            labelList()
         );
         labelIOList boundaryProcAddressing
         (
@@ -816,75 +1027,77 @@ int main(int argc, char *argv[])
                 mesh.pointsInstance(),
                 polyMesh::meshSubDir,
                 mesh,
-                IOobject::READ_IF_PRESENT,
+                (dryrun ? IOobject::NO_READ : IOobject::READ_IF_PRESENT),
                 IOobject::AUTO_WRITE
             ),
-            labelList(0)
+            labelList()
         );
 
 
-        // Read objects in time directory
-        IOobjectList objects(mesh, runTime.timeName());
+        // List of objects read from time directory
+        // List of stored objects to clear from mesh
+
+        IOobjectList objects;
+        DynamicList<regIOobject*> storedObjects;
+
+        if (!dryrun && doFields)
+        {
+            objects = IOobjectList(mesh, runTime.timeName());
+            storedObjects.reserve(objects.size());
+
+            const predicates::always nameMatcher;
+
+            // Read GeometricFields
+
+            #undef  ReadFields
+            #define ReadFields(FieldType)                                     \
+            readFields<FieldType>(mesh, objects, nameMatcher, storedObjects);
+
+            // Read volume fields
+            ReadFields(volScalarField);
+            ReadFields(volVectorField);
+            ReadFields(volSphericalTensorField);
+            ReadFields(volSymmTensorField);
+            ReadFields(volTensorField);
+
+            // Read internal fields
+            ReadFields(volScalarField::Internal);
+            ReadFields(volVectorField::Internal);
+            ReadFields(volSphericalTensorField::Internal);
+            ReadFields(volSymmTensorField::Internal);
+            ReadFields(volTensorField::Internal);
+
+            // Read surface fields
+            ReadFields(surfaceScalarField);
+            ReadFields(surfaceVectorField);
+            ReadFields(surfaceSphericalTensorField);
+            ReadFields(surfaceSymmTensorField);
+            ReadFields(surfaceTensorField);
 
 
-        // Read vol fields.
+            // Read point fields
+            const pointMesh& pMesh = pointMesh::New(mesh);
+            #undef  ReadPointFields
+            #define ReadPointFields(FieldType)                                \
+            readFields<FieldType>(pMesh, objects, nameMatcher, storedObjects);
 
-        PtrList<volScalarField> vsFlds;
-        ReadFields(mesh, objects, vsFlds);
+            ReadPointFields(pointScalarField);
+            ReadPointFields(pointVectorField);
+            ReadPointFields(pointSphericalTensorField);
+            ReadPointFields(pointSymmTensorField);
+            ReadPointFields(pointTensorField);
 
-        PtrList<volVectorField> vvFlds;
-        ReadFields(mesh, objects, vvFlds);
-
-        PtrList<volSphericalTensorField> vstFlds;
-        ReadFields(mesh, objects, vstFlds);
-
-        PtrList<volSymmTensorField> vsymtFlds;
-        ReadFields(mesh, objects, vsymtFlds);
-
-        PtrList<volTensorField> vtFlds;
-        ReadFields(mesh, objects, vtFlds);
-
-
-        // Read surface fields.
-
-        PtrList<surfaceScalarField> ssFlds;
-        ReadFields(mesh, objects, ssFlds);
-
-        PtrList<surfaceVectorField> svFlds;
-        ReadFields(mesh, objects, svFlds);
-
-        PtrList<surfaceSphericalTensorField> sstFlds;
-        ReadFields(mesh, objects, sstFlds);
-
-        PtrList<surfaceSymmTensorField> ssymtFlds;
-        ReadFields(mesh, objects, ssymtFlds);
-
-        PtrList<surfaceTensorField> stFlds;
-        ReadFields(mesh, objects, stFlds);
-
-
-        // Read point fields.
-
-        PtrList<pointScalarField> psFlds;
-        ReadFields(pointMesh::New(mesh), objects, psFlds);
-
-        PtrList<pointVectorField> pvFlds;
-        ReadFields(pointMesh::New(mesh), objects, pvFlds);
-
-        PtrList<pointSphericalTensorField> pstFlds;
-        ReadFields(pointMesh::New(mesh), objects, pstFlds);
-
-        PtrList<pointSymmTensorField> psymtFlds;
-        ReadFields(pointMesh::New(mesh), objects, psymtFlds);
-
-        PtrList<pointTensorField> ptFlds;
-        ReadFields(pointMesh::New(mesh), objects, ptFlds);
+            #undef ReadFields
+            #undef ReadPointFields
+        }
 
 
         // Read sets
         PtrList<cellSet> cellSets;
         PtrList<faceSet> faceSets;
         PtrList<pointSet> pointSets;
+
+        if (!dryrun)
         {
             // Read sets
             IOobjectList objects(mesh, mesh.facesInstance(), "polyMesh/sets");
@@ -901,7 +1114,7 @@ int main(int argc, char *argv[])
         labelList cellOrder;
         labelList faceOrder;
 
-        if (blockSize > 0)
+        if (blockSize > 0 && !doDecompose)
         {
             // Renumbering in two phases. Should be done in one so mapping of
             // fields is done correctly!
@@ -910,15 +1123,14 @@ int main(int argc, char *argv[])
             Info<< "nBlocks   = " << nBlocks << endl;
 
             // Read decompositionMethod dictionary
-            dictionary decomposeDict(renumberDictPtr().subDict("blockCoeffs"));
+            dictionary decomposeDict(renumberDict.subDict("blockCoeffs"));
             decomposeDict.set("numberOfSubdomains", nBlocks);
 
+            // No parallel communication
             const bool oldParRun = UPstream::parRun(false);
 
-            autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
-            (
-                decomposeDict
-            );
+            autoPtr<decompositionMethod> decomposePtr =
+                decompositionMethod::New(decomposeDict);
 
             labelList cellToRegion
             (
@@ -945,7 +1157,10 @@ int main(int argc, char *argv[])
                 << nl << endl;
 
 
-            cellOrder = regionRenumber(renumberPtr(), mesh, cellToRegion);
+            CompactListList<label> regionCellOrder =
+                regionRenumber(renumberPtr(), mesh, cellToRegion);
+
+            cellOrder = regionCellOrder.values();
 
             // Determine new to old face order with new cell numbering
             faceOrder = getRegionFaceOrder
@@ -957,12 +1172,77 @@ int main(int argc, char *argv[])
         }
         else
         {
-            // Determines sorted back to original cell ordering
-            cellOrder = renumberPtr().renumber
-            (
-                mesh,
-                mesh.cellCentres()
-            );
+            if (doDecompose)
+            {
+                // Two-step renumbering.
+                // 1. decompose into regions (like decomposePar)
+                // 2. renumber each sub-region
+
+                // Read decompositionMethod dictionary
+                IOdictionary decomposeDict
+                (
+                    IOobject::selectIO
+                    (
+                        IOobject
+                        (
+                            decompositionModel::canonicalName,
+                            runTime.time().system(),
+                            mesh.regionName(),  // region (if non-default)
+                            runTime,
+                            IOobject::MUST_READ,
+                            IOobject::NO_WRITE,
+                            IOobject::NO_REGISTER
+                        ),
+                        decompDictFile
+                    )
+                );
+
+                // No parallel communication
+                const bool oldParRun = UPstream::parRun(false);
+
+                autoPtr<decompositionMethod> decomposePtr
+                    = decompositionMethod::New(decomposeDict);
+
+                labelList cellToRegion
+                (
+                    decomposePtr().decompose
+                    (
+                        mesh,
+                        mesh.cellCentres()
+                    )
+                );
+
+                UPstream::parRun(oldParRun);  // Restore parallel state
+
+                CompactListList<label> regionCellOrder =
+                    regionRenumber
+                    (
+                        renumberPtr(),
+                        mesh,
+                        cellToRegion,
+                        decomposePtr().nDomains()
+                    );
+
+                cellOrder = regionCellOrder.values();
+
+                // HACK - retain partitioning information for possible use
+                // at a later stage
+                mesh.data().set
+                (
+                    "requested.partition-offsets",
+                    regionCellOrder.offsets()
+                );
+            }
+            else
+            {
+                // Determines sorted back to original cell ordering
+                cellOrder = renumberPtr().renumber
+                (
+                    mesh,
+                    mesh.cellCentres()
+                );
+            }
+
 
             if (sortCoupledFaceCells)
             {
@@ -1003,11 +1283,11 @@ int main(int argc, char *argv[])
                         }
                     }
                 }
-                bndCells.setSize(nBndCells);
-                bndCellMap.setSize(nBndCells);
+                bndCells.resize(nBndCells);
+                bndCellMap.resize(nBndCells);
 
                 // Sort
-                labelList order(sortedOrder(bndCellMap));
+                labelList order(Foam::sortedOrder(bndCellMap));
 
                 // Redo newReverseCellOrder
                 labelList newReverseCellOrder(mesh.nCells(), -1);
@@ -1228,13 +1508,11 @@ int main(int argc, char *argv[])
             );
             reduce(band, maxOp<label>());
             reduce(profile, sumOp<scalar>());
+            reduce(sumSqrIntersect, sumOp<scalar>());
+
             scalar rmsFrontwidth = Foam::sqrt
             (
-                returnReduce
-                (
-                    sumSqrIntersect,
-                    sumOp<scalar>()
-                )/mesh.globalData().nTotalCells()
+                sumSqrIntersect/mesh.globalData().nTotalCells()
             );
 
             Info<< "After renumbering :" << nl
@@ -1243,7 +1521,6 @@ int main(int argc, char *argv[])
 
             if (doFrontWidth)
             {
-
                 Info<< "    rms frontwidth : " << rmsFrontwidth << nl;
             }
 
@@ -1306,110 +1583,191 @@ int main(int argc, char *argv[])
         }
 
 
-        if (overwrite)
+        if (dryrun)
         {
-            mesh.setInstance(oldInstance);
+            if (writeMaps)
+            {
+                fileName file = mesh.time().globalPath()/"renumberMesh";
+
+                const word& regionDir = mesh.regionName();
+
+                if (!regionDir.empty())
+                {
+                    file += "_" + regionDir;
+                }
+
+                // VTK output
+                {
+                    const vtk::vtuCells cells(mesh);
+
+                    vtk::internalMeshWriter writer
+                    (
+                        mesh,
+                        cells,
+                        file,
+                        UPstream::parRun()
+                    );
+
+                    writer.writeGeometry();
+                    writer.beginCellData();
+                    writer.writeCellIDs();
+
+                    labelList ids;
+
+                    if (UPstream::parRun())
+                    {
+                        const label cellOffset =
+                            mesh.globalData().globalMeshCellAddr().localStart();
+
+                        ids.resize(mesh.nCells());
+                        std::transform
+                        (
+                            map().cellMap().cbegin(),
+                            map().cellMap().cend(),
+                            ids.begin(),
+                            [=](const label id) { return (id + cellOffset); }
+                        );
+
+                        writer.writeCellData("origCellID", ids);
+                        writer.writeProcIDs();
+                    }
+                    else
+                    {
+                        writer.writeCellData("origCellID", map().cellMap());
+
+                        // Recover any partition information (in serial)
+                        globalIndex partitionOffsets;
+                        if
+                        (
+                            mesh.data().readIfPresent
+                            (
+                                "requested.partition-offsets",
+                                partitionOffsets.offsets()
+                            )
+                        )
+                        {
+                            if (partitionOffsets.totalSize() != mesh.nCells())
+                            {
+                                WarningInFunction
+                                    << "Requested partition total-size: "
+                                    << partitionOffsets.totalSize()
+                                    << " mesh total-size: "
+                                    << mesh.nCells()
+                                    << " ... ignoring" << endl;
+
+                                partitionOffsets.clear();
+                            }
+                        }
+
+                        ids.resize(partitionOffsets.totalSize());
+
+                        for (const label proci : partitionOffsets.allProcs())
+                        {
+                            ids.slice(partitionOffsets.range(proci)) = proci;
+                        }
+
+                        if (!partitionOffsets.empty())
+                        {
+                            writer.writeCellData("procID", ids);
+                        }
+                    }
+
+                    Info<< "Wrote renumbered mesh to "
+                        << mesh.time().relativePath(writer.output())
+                        << " for visualization."
+                        << endl;
+                }
+            }
         }
         else
         {
-            mesh.setInstance(runTime.timeName());
+            if (overwrite)
+            {
+                mesh.setInstance(oldInstance);
+            }
+            else
+            {
+                mesh.setInstance(runTime.timeName());
+            }
+
+            Info<< "Writing mesh to " << mesh.facesInstance() << endl;
+
+            // Remove old procAddressing files
+            processorMeshes::removeFiles(mesh);
+
+            // Update refinement data
+
+            hexRef8Data refData
+            (
+                IOobject
+                (
+                    "dummy",
+                    mesh.facesInstance(),
+                    polyMesh::meshSubDir,
+                    mesh,
+                    (dryrun ? IOobject::NO_READ : IOobject::READ_IF_PRESENT),
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
+                )
+            );
+            refData.updateMesh(map());
+            refData.write();
+
+            // Update sets
+            topoSet::updateMesh(mesh.facesInstance(), map(), cellSets);
+            topoSet::updateMesh(mesh.facesInstance(), map(), faceSets);
+            topoSet::updateMesh(mesh.facesInstance(), map(), pointSets);
+
+            mesh.write();
+
+            if (writeMaps)
+            {
+                // For debugging: write out region
+                createScalarField
+                (
+                    mesh,
+                    "origCellID",
+                    map().cellMap()
+                )().write();
+
+                createScalarField
+                (
+                    mesh,
+                    "cellID",
+                    identity(mesh.nCells())
+                )().write();
+
+                Info<< nl
+                    << "Wrote current cellID and origCellID as volScalarField"
+                    << " for use in postprocessing." << nl << endl;
+
+                IOobject meshMapIO
+                (
+                    "map-name",
+                    mesh.facesInstance(),
+                    polyMesh::meshSubDir,
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
+                );
+
+                meshMapIO.resetHeader("cellMap");
+                IOListRef<label>(meshMapIO, map().cellMap()).write();
+
+                meshMapIO.resetHeader("faceMap");
+                IOListRef<label>(meshMapIO, map().faceMap()).write();
+
+                meshMapIO.resetHeader("pointMap");
+                IOListRef<label>(meshMapIO, map().pointMap()).write();
+            }
         }
 
-
-        Info<< "Writing mesh to " << mesh.facesInstance() << endl;
-
-        // Remove old procAddressing files
-        processorMeshes::removeFiles(mesh);
-
-        // Update refinement data
-        hexRef8Data refData
-        (
-            IOobject
-            (
-                "dummy",
-                mesh.facesInstance(),
-                polyMesh::meshSubDir,
-                mesh,
-                IOobject::READ_IF_PRESENT,
-                IOobject::NO_WRITE,
-                IOobject::NO_REGISTER
-            )
-        );
-        refData.updateMesh(map());
-        refData.write();
-
-        // Update sets
-        topoSet::updateMesh(mesh.facesInstance(), map(), cellSets);
-        topoSet::updateMesh(mesh.facesInstance(), map(), faceSets);
-        topoSet::updateMesh(mesh.facesInstance(), map(), pointSets);
-
-        mesh.write();
-
-        if (writeMaps)
+        // Cleanup loaded fields
+        while (!storedObjects.empty())
         {
-            // For debugging: write out region
-            createScalarField
-            (
-                mesh,
-                "origCellID",
-                map().cellMap()
-            )().write();
-
-            createScalarField
-            (
-                mesh,
-                "cellID",
-                identity(mesh.nCells())
-            )().write();
-
-            Info<< nl
-                << "Written current cellID and origCellID as volScalarField"
-                << " for use in postprocessing." << nl << endl;
-
-            labelIOList
-            (
-                IOobject
-                (
-                    "cellMap",
-                    mesh.facesInstance(),
-                    polyMesh::meshSubDir,
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    IOobject::NO_REGISTER
-                ),
-                map().cellMap()
-            ).write();
-
-            labelIOList
-            (
-                IOobject
-                (
-                    "faceMap",
-                    mesh.facesInstance(),
-                    polyMesh::meshSubDir,
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    IOobject::NO_REGISTER
-                ),
-                map().faceMap()
-            ).write();
-
-            labelIOList
-            (
-                IOobject
-                (
-                    "pointMap",
-                    mesh.facesInstance(),
-                    polyMesh::meshSubDir,
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    IOobject::NO_REGISTER
-                ),
-                map().pointMap()
-            ).write();
+            storedObjects.back()->checkOut();
+            storedObjects.pop_back();
         }
     }
 
diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRunWrite.C b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRunWrite.C
index 56924ebad3..86abe22562 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRunWrite.C
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecompositionDryRunWrite.C
@@ -83,6 +83,7 @@ void Foam::domainDecompositionDryRun::writeVTK
     writer.writeGeometry();
     writer.beginCellData();
     writer.writeCellData("procID", cellToProc);
+    writer.writeCellIDs();
 
     Info<< "Wrote decomposition to "
         << this->mesh().time().relativePath(writer.output())
diff --git a/bin/tools/CleanFunctions b/bin/tools/CleanFunctions
index 4b54ca4b34..bda0428f45 100644
--- a/bin/tools/CleanFunctions
+++ b/bin/tools/CleanFunctions
@@ -181,7 +181,7 @@ cleanCase()
     # Debug output (blockMesh, decomposePar)
     rm -f \
         blockTopology.vtu blockFaces.vtp blockTopology.obj blockCentres.obj \
-        cellDist.vtu \
+        cellDist.vtu decomposePar.vtu renumberMesh.vtu \
         0/cellDist
 
     # From mpirunDebug
diff --git a/etc/caseDicts/annotated/renumberMeshDict b/etc/caseDicts/annotated/renumberMeshDict
index 00ce52e885..80f960d014 100644
--- a/etc/caseDicts/annotated/renumberMeshDict
+++ b/etc/caseDicts/annotated/renumberMeshDict
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2312                                 |
+|  \\    /   O peration     | Version:  v2406                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -15,8 +15,12 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Write maps from renumbered back to original mesh
-writeMaps false;
+// Additional libraries
+//libs        (SloanRenumber);
+//libs        (zoltanRenumber);
+
+// Optional entry: write maps from renumbered back to original mesh
+writeMaps   false;
 
 // Optional entry: sort cells on coupled boundaries to last for use with
 // e.g. nonBlockingGaussSeidel.
@@ -36,61 +40,66 @@ sortCoupledFaceCells false;
 
 
 method          CuthillMcKee;
-//method          Sloan;
+//method          RCM;  // == reverseCuthillMcKee;
+//method          Sloan;        //<-  libs (zoltanRenumber);
 //method          manual;
 //method          random;
 //method          structured;
 //method          spring;
-//method          zoltan;             // only if compiled with zoltan support
+//method          zoltan;       //<-  libs (zoltanRenumber);
 
 //CuthillMcKeeCoeffs
 //{
-//    // Reverse CuthillMcKee (RCM) or plain
-//    reverse true;
+//    // Plain or reverse CuthillMcKee (RCM)
+//    reverse     true;
 //}
 
 manualCoeffs
 {
     // In system directory: new-to-original (i.e. order) labelIOList
-    dataFile "cellMap";
+    dataFile    "cellMap";
 }
 
 
 // For extruded (i.e. structured in one direction) meshes
 structuredCoeffs
 {
-    // Patches that mesh was extruded from. These determine the starting
-    // layer of cells
-    patches (movingWall);
+    // Patches that mesh was extruded from.
+    // These determine the starting layer of cells
+    patches     (movingWall);
+
     // Method to renumber the starting layer of cells
-    method  random;
+    method      random;
 
     // Renumber in columns (depthFirst) or in layers
-    depthFirst true;
+    depthFirst  true;
 
     // Reverse ordering
-    reverse false;
+    reverse     false;
 }
 
 
 springCoeffs
 {
     // Maximum jump of cell indices. Is fraction of number of cells
-    maxCo 0.01;
+    maxCo       0.01;
 
     // Limit the amount of movement; the fraction maxCo gets decreased
     // with every iteration
     freezeFraction 0.999;
 
     // Maximum number of iterations
-    maxIter 1000;
+    maxIter    1000;
+
+    // Enable/disable verbosity
+    verbose    true;
 }
 
 
 blockCoeffs
 {
-    method          scotch;
-    //method hierarchical;
+    method      scotch;
+    //method      hierarchical;
     //hierarchicalCoeffs
     //{
     //    n           (1 2 1);
diff --git a/etc/config.sh/zoltan b/etc/config.sh/zoltan
new file mode 100644
index 0000000000..f93cb8a4b7
--- /dev/null
+++ b/etc/config.sh/zoltan
@@ -0,0 +1,34 @@
+#----------------------------------*-sh-*--------------------------------------
+# =========                 |
+# \\      /  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.
+#
+# File
+#     etc/config.sh/zoltan
+#     - sourced during wmake process only.
+#
+# Description
+#     Setup for ZOLTAN include/libraries (usually ThirdParty installation).
+#
+#     To disable its use:               ZOLTAN_VERSION=zoltan-none
+#     For system-wide installations:    ZOLTAN_VERSION=zoltan-system
+#
+# Note
+#     A csh version is not needed, since the values here are only sourced
+#     during the wmake process
+#
+#------------------------------------------------------------------------------
+# USER EDITABLE PART: Changes made here may be lost with the next upgrade
+
+export ZOLTAN_VERSION=Zoltan-3.901
+export ZOLTAN_ARCH_PATH=$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER$WM_PRECISION_OPTION$WM_LABEL_OPTION/$ZOLTAN_VERSION
+
+# END OF (NORMAL) USER EDITABLE PART
+#------------------------------------------------------------------------------
diff --git a/src/Allwmake b/src/Allwmake
index 96071901af..3d2e2ac222 100755
--- a/src/Allwmake
+++ b/src/Allwmake
@@ -72,6 +72,7 @@ wmakeLnInclude -u fvOptions
 wmake $targetType lagrangian/basic
 wmake $targetType lagrangian/distributionModels
 
+renumber/Allwmake $targetType $*
 parallel/Allwmake $targetType $*
 
 wmake $targetType dynamicFvMesh
@@ -100,7 +101,6 @@ wmake $targetType fvMotionSolver
 
 # snappyHexMesh uses overset voxelMesh
 mesh/Allwmake $targetType $*
-renumber/Allwmake $targetType $*
 fvAgglomerationMethods/Allwmake $targetType $*
 wmake $targetType waveModels
 
diff --git a/src/renumber/Allwmake b/src/renumber/Allwmake
index 0e5edabb3e..a367a8322e 100755
--- a/src/renumber/Allwmake
+++ b/src/renumber/Allwmake
@@ -21,6 +21,7 @@ fi
 warning="==> skip zoltanRenumber"
 if have_zoltan
 then
+    echo "zoltan - $ZOLTAN_ARCH_PATH"
     wmake $targetType zoltanRenumber || echo "$warning (build issues detected)"
 else
     echo "$warning (no library)"
diff --git a/src/renumber/SloanRenumber/SloanRenumber.C b/src/renumber/SloanRenumber/SloanRenumber.C
index f3dea4f762..9ef6ca08ed 100644
--- a/src/renumber/SloanRenumber/SloanRenumber.C
+++ b/src/renumber/SloanRenumber/SloanRenumber.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2017 OpenFOAM Foundation
-    Copyright (C) 2020-2023 OpenCFD Ltd.
+    Copyright (C) 2020-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -74,6 +74,9 @@ typedef adjacency_list
 typedef graph_traits<Graph>::vertex_descriptor Vertex;
 typedef graph_traits<Graph>::vertices_size_type size_type;
 
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
 namespace Foam
 {
     defineTypeNameAndDebug(SloanRenumber, 0);
@@ -89,6 +92,13 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::SloanRenumber::SloanRenumber(const bool reverse)
+:
+    renumberMethod(),
+    reverse_(reverse)
+{}
+
+
 Foam::SloanRenumber::SloanRenumber(const dictionary& dict)
 :
     renumberMethod(dict),
diff --git a/src/renumber/SloanRenumber/SloanRenumber.H b/src/renumber/SloanRenumber/SloanRenumber.H
index 7f61b3912b..46aa185356 100644
--- a/src/renumber/SloanRenumber/SloanRenumber.H
+++ b/src/renumber/SloanRenumber/SloanRenumber.H
@@ -55,16 +55,8 @@ class SloanRenumber
 {
     // Private Data
 
-        const bool reverse_;
-
-
-    // Private Member Functions
-
-        //- No copy construct
-        SloanRenumber(const SloanRenumber&) = delete;
-
-        //- No copy assignment
-        void operator=(const SloanRenumber&) = delete;
+        //- Use reverse indexing
+        bool reverse_;
 
 
 public:
@@ -75,6 +67,9 @@ public:
 
     // Constructors
 
+        //- Default construct, optionally with reverse
+        explicit SloanRenumber(const bool reverse = false);
+
         //- Construct given the renumber dictionary
         explicit SloanRenumber(const dictionary& dict);
 
@@ -85,6 +80,13 @@ public:
 
     // Member Functions
 
+        //- Toggle reverse on/off
+        void reverse(bool on) noexcept { reverse_ = on; }
+
+        //- The renumbering method does not require a polyMesh
+        virtual bool needs_mesh() const { return false; }
+
+
         //- Return the order in which cells need to be visited
         //- (ie. from ordered back to original cell label).
         //  This is only defined for geometric renumberMethods.
diff --git a/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C b/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
index c5a867decd..241be2cbe3 100644
--- a/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
+++ b/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
@@ -30,11 +30,12 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "bandCompression.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(CuthillMcKeeRenumber, 0);
+    defineTypeName(CuthillMcKeeRenumber);
+    defineTypeName(reverseCuthillMcKeeRenumber);
 
     addToRunTimeSelectionTable
     (
@@ -42,11 +43,34 @@ namespace Foam
         CuthillMcKeeRenumber,
         dictionary
     );
+
+    addToRunTimeSelectionTable
+    (
+        renumberMethod,
+        reverseCuthillMcKeeRenumber,
+        dictionary
+    );
+
+    // Select under the name "RCM"
+    addNamedToRunTimeSelectionTable
+    (
+        renumberMethod,
+        reverseCuthillMcKeeRenumber,
+        dictionary,
+        RCM
+    );
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::CuthillMcKeeRenumber::CuthillMcKeeRenumber(const bool reverse)
+:
+    renumberMethod(),
+    reverse_(reverse)
+{}
+
+
 Foam::CuthillMcKeeRenumber::CuthillMcKeeRenumber(const dictionary& dict)
 :
     renumberMethod(dict),
@@ -58,19 +82,45 @@ Foam::CuthillMcKeeRenumber::CuthillMcKeeRenumber(const dictionary& dict)
 {}
 
 
+Foam::CuthillMcKeeRenumber::CuthillMcKeeRenumber
+(
+    const dictionary& dict,
+    const bool reverse
+)
+:
+    renumberMethod(dict),
+    reverse_(reverse)
+{}
+
+
+Foam::reverseCuthillMcKeeRenumber::reverseCuthillMcKeeRenumber()
+:
+    CuthillMcKeeRenumber(true)  // reverse = true
+{}
+
+
+Foam::reverseCuthillMcKeeRenumber::reverseCuthillMcKeeRenumber
+(
+    const dictionary& dict
+)
+:
+    CuthillMcKeeRenumber(dict, true)  // reverse = true
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::labelList Foam::CuthillMcKeeRenumber::renumber
 (
     const polyMesh& mesh,
-    const pointField& points
+    const pointField&
 ) const
 {
     labelList orderedToOld = meshTools::bandCompression(mesh);
 
     if (reverse_)
     {
-        reverse(orderedToOld);
+        Foam::reverse(orderedToOld);
     }
 
     return orderedToOld;
@@ -88,7 +138,7 @@ Foam::labelList Foam::CuthillMcKeeRenumber::renumber
 
     if (reverse_)
     {
-        reverse(orderedToOld);
+        Foam::reverse(orderedToOld);
     }
 
     return orderedToOld;
@@ -98,14 +148,14 @@ Foam::labelList Foam::CuthillMcKeeRenumber::renumber
 Foam::labelList Foam::CuthillMcKeeRenumber::renumber
 (
     const CompactListList<label>& cellCells,
-    const pointField& cc
+    const pointField&
 ) const
 {
     labelList orderedToOld = meshTools::bandCompression(cellCells);
 
     if (reverse_)
     {
-        reverse(orderedToOld);
+        Foam::reverse(orderedToOld);
     }
 
     return orderedToOld;
@@ -115,14 +165,14 @@ Foam::labelList Foam::CuthillMcKeeRenumber::renumber
 Foam::labelList Foam::CuthillMcKeeRenumber::renumber
 (
     const labelListList& cellCells,
-    const pointField& points
+    const pointField&
 ) const
 {
     labelList orderedToOld = meshTools::bandCompression(cellCells);
 
     if (reverse_)
     {
-        reverse(orderedToOld);
+        Foam::reverse(orderedToOld);
     }
 
     return orderedToOld;
diff --git a/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.H b/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.H
index fd6df4cc27..45a07be355 100644
--- a/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.H
+++ b/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2022-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,7 @@ Class
     Foam::CuthillMcKeeRenumber
 
 Description
-    Cuthill-McKee renumbering
+    Cuthill-McKee renumbering (CM or RCM)
 
 SourceFiles
     CuthillMcKeeRenumber.C
@@ -53,29 +53,28 @@ class CuthillMcKeeRenumber
 {
     // Private Data
 
-        const bool reverse_;
-
-
-    // Private Member Functions
-
-        //- No copy construct
-        CuthillMcKeeRenumber(const CuthillMcKeeRenumber&) = delete;
-
-        //- No copy assignment
-        void operator=(const CuthillMcKeeRenumber&) = delete;
+        //- Use reverse indexing
+        bool reverse_;
 
 
 public:
 
     //- Runtime type information
-    TypeName("CuthillMcKee");
+    TypeNameNoDebug("CuthillMcKee");
 
 
     // Constructors
 
+        //- Default construct, optionally with reverse
+        explicit CuthillMcKeeRenumber(const bool reverse = false);
+
         //- Construct given the renumber dictionary
         explicit CuthillMcKeeRenumber(const dictionary& dict);
 
+        //- Construct given the renumber dictionary (ignored)
+        //- and specified reverse handling
+        CuthillMcKeeRenumber(const dictionary& dict, const bool reverse);
+
 
     //- Destructor
     virtual ~CuthillMcKeeRenumber() = default;
@@ -83,6 +82,13 @@ public:
 
     // Member Functions
 
+        //- Toggle reverse on/off
+        void reverse(bool on) noexcept { reverse_ = on; }
+
+        //- The renumbering method does not require a polyMesh
+        virtual bool needs_mesh() const { return false; }
+
+
         //- Return the order in which cells need to be visited
         //- (ie. from ordered back to original cell label).
         //  This is only defined for geometric renumberMethods.
@@ -97,8 +103,10 @@ public:
         //  Use the mesh connectivity (if needed)
         virtual labelList renumber
         (
+            //! Mesh connectivity (see globalMeshData::calcCellCells)
             const polyMesh& mesh,
-            const pointField& cc
+            //! \em ignored
+            const pointField& cellCentres
         ) const;
 
         //- Return the order in which cells need to be visited
@@ -109,14 +117,17 @@ public:
         (
             const labelList& cellCells,
             const labelList& offsets,
-            const pointField& cc
+            //! \em ignored
+            const pointField& cellCentres
         ) const;
 
         //- Return the order in which cells need to be visited
         //- (ie. from ordered back to original cell label).
         virtual labelList renumber
         (
+            //! Mesh connectivity
             const CompactListList<label>& cellCells,
+            //! \em ignored
             const pointField& cellCentres
         ) const;
 
@@ -126,12 +137,43 @@ public:
         //  - the connections are across coupled patches
         virtual labelList renumber
         (
+            //! Mesh connectivity
             const labelListList& cellCells,
+            //! \em ignored
             const pointField& cellCentres
         ) const;
 };
 
 
+/*---------------------------------------------------------------------------*\
+                   Class reverseCuthillMcKeeRenumber Declaration
+\*---------------------------------------------------------------------------*/
+
+//- Reverse Cuthill-McKee renumbering
+class reverseCuthillMcKeeRenumber
+:
+    public CuthillMcKeeRenumber
+{
+public:
+
+    //- Runtime type information
+    TypeNameNoDebug("reverseCuthillMcKee");
+
+
+    // Constructors
+
+        //- Default construct
+        reverseCuthillMcKeeRenumber();
+
+        //- Construct given the renumber dictionary (ignored)
+        explicit reverseCuthillMcKeeRenumber(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~reverseCuthillMcKeeRenumber() = default;
+};
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/renumber/renumberMethods/Make/files b/src/renumber/renumberMethods/Make/files
index 6ee4924d33..863d3d2100 100644
--- a/src/renumber/renumberMethods/Make/files
+++ b/src/renumber/renumberMethods/Make/files
@@ -1,7 +1,10 @@
 renumberMethod/renumberMethod.C
+
 manualRenumber/manualRenumber.C
-CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
+noRenumber/noRenumber.C
 randomRenumber/randomRenumber.C
+CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
+
 springRenumber/springRenumber.C
 structuredRenumber/structuredRenumber.C
 structuredRenumber/OppositeFaceCellWaveBase.C
diff --git a/src/renumber/renumberMethods/manualRenumber/manualRenumber.C b/src/renumber/renumberMethods/manualRenumber/manualRenumber.C
index bce6268cef..61dc497239 100644
--- a/src/renumber/renumberMethods/manualRenumber/manualRenumber.C
+++ b/src/renumber/renumberMethods/manualRenumber/manualRenumber.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,14 +28,13 @@ License
 
 #include "manualRenumber.H"
 #include "addToRunTimeSelectionTable.H"
-#include "IFstream.H"
 #include "labelIOList.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(manualRenumber, 0);
+    defineTypeName(manualRenumber);
 
     addToRunTimeSelectionTable
     (
@@ -48,12 +47,19 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::manualRenumber::manualRenumber(const fileName& file)
+:
+    renumberMethod(),
+    dataFile_(file)
+{}
+
+
 Foam::manualRenumber::manualRenumber(const dictionary& dict)
 :
     renumberMethod(dict),
     dataFile_
     (
-        dict.optionalSubDict(typeName+"Coeffs").get<fileName>("dataFile")
+        dict.optionalSubDict(typeName + "Coeffs").get<fileName>("dataFile")
     )
 {}
 
diff --git a/src/renumber/renumberMethods/manualRenumber/manualRenumber.H b/src/renumber/renumberMethods/manualRenumber/manualRenumber.H
index f6dbb1abf5..208a6d959b 100644
--- a/src/renumber/renumberMethods/manualRenumber/manualRenumber.H
+++ b/src/renumber/renumberMethods/manualRenumber/manualRenumber.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2022-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,23 +55,17 @@ class manualRenumber
 
         fileName dataFile_;
 
-
-    // Private Member Functions
-
-        //- No copy construct
-        manualRenumber(const manualRenumber&) = delete;
-
-        //- No copy assignment
-        void operator=(const manualRenumber&) = delete;
-
 public:
 
     //- Runtime type information
-    TypeName("manual");
+    TypeNameNoDebug("manual");
 
 
     // Constructors
 
+        //- Construct with given data file
+        explicit manualRenumber(const fileName& file);
+
         //- Construct given the renumber dictionary
         explicit manualRenumber(const dictionary& dict);
 
@@ -82,6 +76,10 @@ public:
 
     // Member Functions
 
+        //- The renumbering method needs a polyMesh (for its IOobject)
+        virtual bool needs_mesh() const { return true; }
+
+
         //- Return the order in which cells need to be visited
         //- (ie. from ordered back to original cell label).
         virtual labelList renumber(const pointField&) const
diff --git a/src/renumber/renumberMethods/noRenumber/noRenumber.C b/src/renumber/renumberMethods/noRenumber/noRenumber.C
new file mode 100644
index 0000000000..6387830399
--- /dev/null
+++ b/src/renumber/renumberMethods/noRenumber/noRenumber.C
@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "noRenumber.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeName(noRenumber);
+
+    addNamedToRunTimeSelectionTable
+    (
+        renumberMethod,
+        noRenumber,
+        dictionary,
+        none
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::noRenumber::noRenumber()
+:
+    renumberMethod()
+{}
+
+
+Foam::noRenumber::noRenumber(const dictionary& dict)
+:
+    renumberMethod(dict)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::noRenumber::renumber
+(
+    const pointField& cellCentres
+) const
+{
+    return Foam::identity(cellCentres.size());
+}
+
+
+Foam::labelList Foam::noRenumber::renumber
+(
+    const polyMesh& mesh,
+    const pointField&
+) const
+{
+    return Foam::identity(mesh.nCells());
+}
+
+
+Foam::labelList Foam::noRenumber::renumber
+(
+    const CompactListList<label>& cellCells,
+    const pointField&
+) const
+{
+    return Foam::identity(cellCells.size());
+}
+
+
+Foam::labelList Foam::noRenumber::renumber
+(
+    const labelListList& cellCells,
+    const pointField&
+) const
+{
+    return Foam::identity(cellCells.size());
+}
+
+
+// ************************************************************************* //
diff --git a/src/renumber/renumberMethods/noRenumber/noRenumber.H b/src/renumber/renumberMethods/noRenumber/noRenumber.H
new file mode 100644
index 0000000000..98273362da
--- /dev/null
+++ b/src/renumber/renumberMethods/noRenumber/noRenumber.H
@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::noRenumber
+
+Description
+    A dummy renumber method, selected as \c none.
+
+    Method coefficients: \a none
+
+SourceFiles
+    noRenumber.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_noRenumber_H
+#define Foam_noRenumber_H
+
+#include "renumberMethod.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class noRenumber Declaration
+\*---------------------------------------------------------------------------*/
+
+class noRenumber
+:
+    public renumberMethod
+{
+public:
+
+    //- Runtime type information
+    TypeNameNoDebug("none");
+
+
+    // Constructors
+
+        //- Default construct
+        noRenumber();
+
+        //- Construct given the renumber dictionary (unused)
+        explicit noRenumber(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~noRenumber() = default;
+
+
+    // Member Functions
+
+        //- The renumbering method does not require a polyMesh
+        virtual bool needs_mesh() const { return false; }
+
+
+        //- Return the order in which cells need to be visited
+        //- (ie. from ordered back to original cell label).
+        virtual labelList renumber(const pointField&) const;
+
+        //- Return the order in which cells need to be visited
+        //- (ie. from ordered back to original cell label).
+        virtual labelList renumber
+        (
+            //! Mesh number of cells
+            const polyMesh& mesh,
+            //! \em ignored
+            const pointField& cellCentres
+        ) const;
+
+        //- Return the order in which cells need to be visited
+        //- (ie. from ordered back to original cell label).
+        virtual labelList renumber
+        (
+            //! Mesh connectivity for number of cells
+            const CompactListList<label>& cellCells,
+            //! \em ignored
+            const pointField& cellCentres
+        ) const;
+
+        //- Return the order in which cells need to be visited
+        //- (ie. from ordered back to original cell label).
+        virtual labelList renumber
+        (
+            //! Mesh connectivity for number of cells
+            const labelListList& cellCells,
+            //! \em ignored
+            const pointField& cellCentres
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/renumber/renumberMethods/randomRenumber/randomRenumber.C b/src/renumber/renumberMethods/randomRenumber/randomRenumber.C
index a18be4d3c8..4ae111352d 100644
--- a/src/renumber/renumberMethods/randomRenumber/randomRenumber.C
+++ b/src/renumber/renumberMethods/randomRenumber/randomRenumber.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,11 +30,11 @@ License
 #include "Random.H"
 #include "addToRunTimeSelectionTable.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(randomRenumber, 0);
+    defineTypeName(randomRenumber);
 
     addToRunTimeSelectionTable
     (
@@ -45,8 +45,39 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+labelList randomMap(const label nCells)
+{
+    Random rndGen(0);
+
+    // Full coverage
+    labelList newToOld(Foam::identity(nCells));
+
+    // Fisher-Yates shuffle algorithm
+    for (label i = nCells - 1; i > 0; --i)
+    {
+        label j = rndGen.position<label>(0, i);
+        std::swap(newToOld[i], newToOld[j]);
+    }
+
+    return newToOld;
+}
+
+} // End namespace Foam
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::randomRenumber::randomRenumber()
+:
+    renumberMethod()
+{}
+
+
 Foam::randomRenumber::randomRenumber(const dictionary& dict)
 :
     renumberMethod(dict)
@@ -57,52 +88,40 @@ Foam::randomRenumber::randomRenumber(const dictionary& dict)
 
 Foam::labelList Foam::randomRenumber::renumber
 (
-    const pointField& points
+    const pointField& cellCentres
 ) const
 {
-    Random rndGen(0);
-
-    labelList newToOld(identity(points.size()));
-
-    for (label iter = 0; iter < 10; iter++)
-    {
-        forAll(newToOld, i)
-        {
-            label j = rndGen.position<label>(0, newToOld.size()-1);
-            std::swap(newToOld[i], newToOld[j]);
-        }
-    }
-    return newToOld;
+    return randomMap(cellCentres.size());
 }
 
 
 Foam::labelList Foam::randomRenumber::renumber
 (
     const polyMesh& mesh,
-    const pointField& points
+    const pointField&
 ) const
 {
-    return renumber(points);
+    return randomMap(mesh.nCells());
 }
 
 
 Foam::labelList Foam::randomRenumber::renumber
 (
     const CompactListList<label>& cellCells,
-    const pointField& points
+    const pointField&
 ) const
 {
-    return renumber(points);
+    return randomMap(cellCells.size());
 }
 
 
 Foam::labelList Foam::randomRenumber::renumber
 (
     const labelListList& cellCells,
-    const pointField& points
+    const pointField&
 ) const
 {
-    return renumber(points);
+    return randomMap(cellCells.size());
 }
 
 
diff --git a/src/renumber/renumberMethods/randomRenumber/randomRenumber.H b/src/renumber/renumberMethods/randomRenumber/randomRenumber.H
index 233a98638c..5769be0ce7 100644
--- a/src/renumber/renumberMethods/randomRenumber/randomRenumber.H
+++ b/src/renumber/renumberMethods/randomRenumber/randomRenumber.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2022-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -51,24 +51,18 @@ class randomRenumber
 :
     public renumberMethod
 {
-    // Private Member Functions
-
-        //- No copy construct
-        randomRenumber(const randomRenumber&) = delete;
-
-        //- No copy assignment
-        void operator=(const randomRenumber&) = delete;
-
-
 public:
 
     //- Runtime type information
-    TypeName("random");
+    TypeNameNoDebug("random");
 
 
     // Constructors
 
-        //- Construct given the renumber dictionary
+        //- Default construct
+        randomRenumber();
+
+        //- Construct given the renumber dictionary (unused)
         explicit randomRenumber(const dictionary& dict);
 
 
@@ -78,19 +72,21 @@ public:
 
     // Member Functions
 
+        //- The renumbering method does not require a polyMesh
+        virtual bool needs_mesh() const { return false; }
+
+
         //- Return the order in which cells need to be visited
         //- (ie. from ordered back to original cell label).
         virtual labelList renumber(const pointField&) const;
 
-        //- Return the order in which cells need to be visited
-        //- (ie. from ordered back to original cell label).
-        virtual labelList renumber(const polyMesh&, const pointField&) const;
-
         //- Return the order in which cells need to be visited
         //- (ie. from ordered back to original cell label).
         virtual labelList renumber
         (
-            const CompactListList<label>& cellCells,
+            //! Mesh number of cells
+            const polyMesh& mesh,
+            //! \em ignored
             const pointField& cellCentres
         ) const;
 
@@ -98,7 +94,19 @@ public:
         //- (ie. from ordered back to original cell label).
         virtual labelList renumber
         (
+            //! Mesh connectivity for number of cells
+            const CompactListList<label>& cellCells,
+            //! \em ignored
+            const pointField& cellCentres
+        ) const;
+
+        //- Return the order in which cells need to be visited
+        //- (ie. from ordered back to original cell label).
+        virtual labelList renumber
+        (
+            //! Mesh connectivity for number of cells
             const labelListList& cellCells,
+            //! \em ignored
             const pointField& cellCentres
         ) const;
 };
diff --git a/src/renumber/renumberMethods/renumberMethod/renumberMethod.C b/src/renumber/renumberMethods/renumberMethod/renumberMethod.C
index 7be4654d92..7cb556bd8f 100644
--- a/src/renumber/renumberMethods/renumberMethod/renumberMethod.C
+++ b/src/renumber/renumberMethods/renumberMethod/renumberMethod.C
@@ -27,18 +27,30 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "renumberMethod.H"
+#include "dlLibraryTable.H"
 #include "globalMeshData.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(renumberMethod, 0);
+    defineTypeName(renumberMethod);
     defineRunTimeSelectionTable(renumberMethod, dictionary);
 }
 
+
 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
 
+Foam::wordList Foam::renumberMethod::supportedMethods()
+{
+    if (dictionaryConstructorTablePtr_)
+    {
+        return dictionaryConstructorTablePtr_->sortedToc();
+    }
+    return wordList();
+}
+
+
 Foam::autoPtr<Foam::renumberMethod> Foam::renumberMethod::New
 (
     const dictionary& dict
@@ -48,6 +60,9 @@ Foam::autoPtr<Foam::renumberMethod> Foam::renumberMethod::New
 
     //Info<< "Selecting renumberMethod " << methodType << endl;
 
+    // Load any additional libs (verbose = false)
+    dlLibraryTable::libs().open("libs", dict, false);
+
     auto* ctorPtr = dictionaryConstructorTable(methodType);
 
     if (!ctorPtr)
@@ -67,6 +82,16 @@ Foam::autoPtr<Foam::renumberMethod> Foam::renumberMethod::New
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::labelList Foam::renumberMethod::renumber
+(
+    const pointField&
+) const
+{
+    NotImplemented;
+    return labelList();
+}
+
+
 Foam::labelList Foam::renumberMethod::renumber
 (
     const polyMesh& mesh,
diff --git a/src/renumber/renumberMethods/renumberMethod/renumberMethod.H b/src/renumber/renumberMethods/renumberMethod/renumberMethod.H
index bc23c85404..9bfa6dd6c5 100644
--- a/src/renumber/renumberMethods/renumberMethod/renumberMethod.H
+++ b/src/renumber/renumberMethods/renumberMethod/renumberMethod.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2022-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,6 +41,7 @@ SourceFiles
 #include "polyMesh.H"
 #include "pointField.H"
 #include "CompactListList.H"
+#include "wordList.H"
 
 namespace Foam
 {
@@ -51,26 +52,10 @@ namespace Foam
 
 class renumberMethod
 {
-protected:
-
-    // Protected Data
-
-        const dictionary& renumberDict_;
-
-
-    // Protected Member Functions
-
-        //- No copy construct
-        renumberMethod(const renumberMethod&) = delete;
-
-        //- No copy assignment
-        void operator=(const renumberMethod&) = delete;
-
-
 public:
 
     //- Runtime type information
-    TypeName("renumberMethod");
+    TypeNameNoDebug("renumberMethod");
 
 
     // Declare run-time constructor selection tables
@@ -87,23 +72,25 @@ public:
         );
 
 
-    // Selectors
-
-        //- Return a reference to the selected renumbering method
-        static autoPtr<renumberMethod> New
-        (
-            const dictionary& renumberDict
-        );
-
-
     // Constructors
 
-        //- Construct given the renumber dictionary
-        explicit renumberMethod(const dictionary& dict)
-        :
-            renumberDict_(dict)
+        //- Default construct
+        renumberMethod()
         {}
 
+        //- Construct with renumber dictionary (which is currently unused)
+        explicit renumberMethod(const dictionary&)
+        {}
+
+
+    // Selectors
+
+        //- Construct/select a renumbering method
+        static autoPtr<renumberMethod> New(const dictionary& dict);
+
+        //- Return a list of the known methods
+        static wordList supportedMethods();
+
 
     //- Destructor
     virtual ~renumberMethod() = default;
@@ -111,31 +98,25 @@ public:
 
     // Member Functions
 
+        //- Does renumbering method require a polyMesh?
+        virtual bool needs_mesh() const { return false; }
+
+
+    // No topology
+
         //- Return the order in which cells need to be visited
         //- (ie. from ordered back to original cell label).
-        //  This is only defined for geometric renumberMethods.
-        virtual labelList renumber(const pointField&) const
-        {
-            NotImplemented;
-            return labelList();
-        }
+        //  This is only defined for geometric renumber methods.
+        virtual labelList renumber(const pointField&) const;
+
+
+    // Topology provided by mesh
 
         //- Return the order in which cells need to be visited
         //- (ie. from ordered back to original cell label).
         //  Use the mesh connectivity (if needed)
         virtual labelList renumber(const polyMesh&, const pointField&) const;
 
-        //- Return the order in which cells need to be visited
-        //- (ie. from ordered back to original cell label).
-        //  Addressing in losort addressing (= neighbour + offsets into
-        //  neighbour)
-        virtual labelList renumber
-        (
-            const labelList& cellCells,
-            const labelList& offsets,
-            const pointField&
-        ) const;
-
         //- Return the order in which cells need to be visited
         //- (ie. from ordered back to original cell label).
         //  Gets passed agglomeration map (from fine to coarse cells) and coarse
@@ -151,6 +132,20 @@ public:
             const pointField& coarsePoints
         ) const;
 
+
+    // Topology provided explicitly
+
+        //- Return the order in which cells need to be visited
+        //- (ie. from ordered back to original cell label).
+        //  Addressing in losort addressing (= neighbour + offsets into
+        //  neighbour)
+        virtual labelList renumber
+        (
+            const labelList& cellCells,
+            const labelList& offsets,
+            const pointField&
+        ) const;
+
         //- Return the order in which cells need to be visited
         //- (ie. from ordered back to original cell label).
         //  Uses 'unpack' internally, so should be overloaded when possible
diff --git a/src/renumber/renumberMethods/springRenumber/springRenumber.C b/src/renumber/renumberMethods/springRenumber/springRenumber.C
index c331b23ad8..b0063d6b59 100644
--- a/src/renumber/renumberMethods/springRenumber/springRenumber.C
+++ b/src/renumber/renumberMethods/springRenumber/springRenumber.C
@@ -50,10 +50,11 @@ namespace Foam
 Foam::springRenumber::springRenumber(const dictionary& dict)
 :
     renumberMethod(dict),
-    coeffsDict_(dict.optionalSubDict(typeName+"Coeffs")),
+    coeffsDict_(dict.optionalSubDict(typeName + "Coeffs")),
     maxIter_(coeffsDict_.get<label>("maxIter")),
     maxCo_(coeffsDict_.get<scalar>("maxCo")),
-    freezeFraction_(coeffsDict_.get<scalar>("freezeFraction"))
+    freezeFraction_(coeffsDict_.get<scalar>("freezeFraction")),
+    verbose_(coeffsDict_.getOrDefault("verbose", true))
 {}
 
 
@@ -71,10 +72,7 @@ Foam::labelList Foam::springRenumber::renumberImpl
     // Move cells to the average 'position' of their neighbour.
 
     scalarField position(nOldCells);
-    forAll(position, celli)
-    {
-        position[celli] = celli;
-    }
+    std::iota(position.begin(), position.end(), 0);
 
     // Sum force per cell. Also reused for the displacement
     scalarField sumForce(nOldCells);
@@ -90,9 +88,9 @@ Foam::labelList Foam::springRenumber::renumberImpl
         //    << endl;
 
         //Pout<< "Position :" << nl
-        //    << "    min : " << min(position) << nl
-        //    << "    max : " << max(position) << nl
-        //    << "    avg : " << average(position) << nl
+        //    << "    min : " << Foam::min(position) << nl
+        //    << "    max : " << Foam::max(position) << nl
+        //    << "    avg : " << Foam::average(position) << nl
         //    << endl;
 
         // Sum force per cell.
@@ -100,9 +98,8 @@ Foam::labelList Foam::springRenumber::renumberImpl
         for (label oldCelli = 0; oldCelli < nOldCells; ++oldCelli)
         {
             const label celli = oldToNew[oldCelli];
-            const auto& neighbours = cellCells[oldCelli];
 
-            for (const label nbr : neighbours)
+            for (const label nbr : cellCells[oldCelli])
             {
                 const label nbrCelli = oldToNew[nbr];
 
@@ -111,35 +108,38 @@ Foam::labelList Foam::springRenumber::renumberImpl
         }
 
         //Pout<< "Force :" << nl
-        //    << "    min    : " << min(sumForce) << nl
-        //    << "    max    : " << max(sumForce) << nl
-        //    << "    avgMag : " << average(mag(sumForce)) << nl
+        //    << "    min    : " << Foam::min(sumForce) << nl
+        //    << "    max    : " << Foam::max(sumForce) << nl
+        //    << "    avgMag : " << Foam::average(mag(sumForce)) << nl
         //    << "DeltaT : " << deltaT << nl
         //    << endl;
 
         // Limit displacement
         scalar deltaT = maxCo / max(mag(sumForce));
 
-        Info<< "Iter:" << iter
-            << "  maxCo:" << maxCo
-            << "  deltaT:" << deltaT
-            << "  average force:" << average(mag(sumForce)) << endl;
+        if (verbose_)
+        {
+            Info<< "Iter:" << iter
+                << "  maxCo:" << maxCo
+                << "  deltaT:" << deltaT
+                << "  average force:" << average(mag(sumForce)) << endl;
+        }
 
         // Determine displacement
         scalarField& displacement = sumForce;
         displacement *= deltaT;
 
         //Pout<< "Displacement :" << nl
-        //    << "    min    : " << min(displacement) << nl
-        //    << "    max    : " << max(displacement) << nl
-        //    << "    avgMag : " << average(mag(displacement)) << nl
+        //    << "    min    : " << Foam::min(displacement) << nl
+        //    << "    max    : " << Foam::max(displacement) << nl
+        //    << "    avgMag : " << Foam::average(mag(displacement)) << nl
         //    << endl;
 
         // Calculate new position and scale to be within original range
         // (0..nCells-1) for ease of postprocessing.
         position += displacement;
-        position -= min(position);
-        position *= (position.size()-1)/max(position);
+        position -= Foam::min(position);
+        position *= (position.size()-1)/Foam::max(position);
 
         // Slowly freeze.
         maxCo *= freezeFraction_;
@@ -148,7 +148,7 @@ Foam::labelList Foam::springRenumber::renumberImpl
     //writeOBJ("endPosition.obj", cellCells, position);
 
     // Move cells to new position
-    labelList shuffle(sortedOrder(position));
+    labelList shuffle(Foam::sortedOrder(position));
 
     // Reorder oldToNew
     inplaceReorder(shuffle, oldToNew);
diff --git a/src/renumber/renumberMethods/springRenumber/springRenumber.H b/src/renumber/renumberMethods/springRenumber/springRenumber.H
index c7e6ae4f9e..bf0d8ac8b5 100644
--- a/src/renumber/renumberMethods/springRenumber/springRenumber.H
+++ b/src/renumber/renumberMethods/springRenumber/springRenumber.H
@@ -72,6 +72,9 @@ class springRenumber
 
         const scalar freezeFraction_;
 
+        //- Verbose reporting of progress
+        bool verbose_;
+
 
     // Private Member Functions
 
diff --git a/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C b/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C
index e94aa0406b..f5a13168a9 100644
--- a/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C
+++ b/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C
@@ -32,7 +32,7 @@ License
 #include "fvMeshSubset.H"
 #include "OppositeFaceCellWave.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
@@ -59,7 +59,7 @@ Foam::structuredRenumber::structuredRenumber
     patches_(coeffsDict_.get<wordRes>("patches")),
     nLayers_(coeffsDict_.getOrDefault<label>("nLayers", labelMax)),
     depthFirst_(coeffsDict_.get<bool>("depthFirst")),
-    reverse_(coeffsDict_.get<bool>("reverse")),
+    reverse_(coeffsDict_.getOrDefault("reverse", false)),
     method_(renumberMethod::New(coeffsDict_))
 {}
 
@@ -70,10 +70,10 @@ bool Foam::structuredRenumber::layerLess::operator()
 (
     const label a,
     const label b
-)
+) const
 {
-    const topoDistanceData<label>& ta = distance_[a];
-    const topoDistanceData<label>& tb = distance_[b];
+    const auto& ta = distance_[a];
+    const auto& tb = distance_[b];
 
     int dummy;
 
@@ -273,7 +273,7 @@ Foam::labelList Foam::structuredRenumber::renumber
     // Return furthest away cell first
     if (reverse_)
     {
-        reverse(orderedToOld);
+        Foam::reverse(orderedToOld);
     }
 
     return orderedToOld;
diff --git a/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.H b/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.H
index c79ce09a46..f436862832 100644
--- a/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.H
+++ b/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -52,7 +52,7 @@ namespace Foam
 template<class Type> class topoDistanceData;
 
 /*---------------------------------------------------------------------------*\
-                   Class structuredRenumber Declaration
+                     Class structuredRenumber Declaration
 \*---------------------------------------------------------------------------*/
 
 class structuredRenumber
@@ -63,8 +63,7 @@ public:
 
     // Public Classes
 
-        //- Less function class that can be used for sorting according to
-        //  column and layer
+        //- Function class for sorting according to column and layer
         class layerLess
         {
             const bool depthFirst_;
@@ -85,10 +84,12 @@ public:
                 distance_(distance)
             {}
 
-            bool operator()(const label a, const label b);
+            bool operator()(const label a, const label b) const;
         };
 
 
+private:
+
     // Private Data
 
         const dictionary& coeffsDict_;
@@ -131,6 +132,10 @@ public:
 
     // Member Functions
 
+        //- Renumbering method requires a polyMesh
+        virtual bool needs_mesh() const { return true; }
+
+
         //- Return the order in which cells need to be visited
         //- (ie. from ordered back to original cell label).
         //  This is only defined for geometric renumberMethods.
diff --git a/src/renumber/zoltanRenumber/Make/files b/src/renumber/zoltanRenumber/Make/files
index 15659a52d0..91cd6da72b 100644
--- a/src/renumber/zoltanRenumber/Make/files
+++ b/src/renumber/zoltanRenumber/Make/files
@@ -1,3 +1,4 @@
 zoltanRenumber.C
 
-LIB = $(FOAM_LIBBIN)/libzoltanRenumber
+/* Build with MPI dependency */
+LIB = $(FOAM_MPI_LIBBIN)/libzoltanRenumber
diff --git a/src/renumber/zoltanRenumber/Make/options b/src/renumber/zoltanRenumber/Make/options
index 589a3edc8b..54063354a0 100644
--- a/src/renumber/zoltanRenumber/Make/options
+++ b/src/renumber/zoltanRenumber/Make/options
@@ -8,5 +8,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
-    -L$(ZOLTAN_LIB_DIR) -lzoltan \
-    -lmeshTools
+    -lmeshTools \
+    -lrenumberMethods \
+    $(PLIBS) \
+    -L$(ZOLTAN_LIB_DIR) -lzoltan
diff --git a/src/renumber/zoltanRenumber/zoltanRenumber.C b/src/renumber/zoltanRenumber/zoltanRenumber.C
index 7b50f620d0..5cd44e1193 100644
--- a/src/renumber/zoltanRenumber/zoltanRenumber.C
+++ b/src/renumber/zoltanRenumber/zoltanRenumber.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,7 +58,10 @@ SourceFiles
 #include "polyMesh.H"
 #include "globalMeshData.H"
 #include "globalIndex.H"
+#include "CStringList.H"
 #include "uint.H"
+#include <algorithm>
+#include <numeric>
 
 // Include MPI without any C++ bindings
 #ifndef MPICH_SKIP_MPICXX
@@ -71,7 +75,7 @@ SourceFiles
 #include "zoltan.h"
 #include <mpi.h>
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
@@ -86,43 +90,80 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+// [ZOLTAN_NUM_OBJ_FN]
+// The number of graph vertices (locally)
 static int get_number_of_vertices(void *data, int *ierr)
 {
-    const Foam::polyMesh& mesh = *static_cast<const Foam::polyMesh*>(data);
     *ierr = ZOLTAN_OK;
+    const auto& mesh = *static_cast<const Foam::polyMesh*>(data);
+
     return mesh.nCells();
 }
 
 
-static void get_vertex_list(void *data, int sizeGID, int sizeLID,
-            ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID,
-                  int wgt_dim, float *obj_wgts, int *ierr)
+// [ZOLTAN_OBJ_LIST_FN]
+// The ids for the graph vertices
+static void get_vertex_list
+(
+    void *data,
+    int sizeGID,                // (currently unused)
+    int sizeLID,                // (currently unused)
+    ZOLTAN_ID_PTR global_ids,
+    ZOLTAN_ID_PTR local_ids,
+    int wgt_dim,                // (currently unused)
+    float *obj_wgts,            // (currently unused)
+    int *ierr
+)
 {
-    const Foam::polyMesh& mesh = *static_cast<const Foam::polyMesh*>(data);
-
     *ierr = ZOLTAN_OK;
+    const auto& mesh = *static_cast<const Foam::polyMesh*>(data);
 
-  /* In this example, return the IDs of our vertices, but no weights.
-   * Zoltan will assume equally weighted vertices.
-   */
+    const auto nCells = mesh.nCells();
 
-    wgt_dim = 0;
-    obj_wgts = nullptr;
+    // Offset for globally unique IDs
+    const auto myProci = Foam::UPstream::myProcNo();
+    const auto myProcOffset =
+        mesh.globalData().globalMeshCellAddr().localStart(myProci);
 
-    for (Foam::label i=0; i<mesh.nCells(); i++)
-    {
-        globalID[i] = i;        // should be global
-        localID[i] = i;
-    }
+    // Global indices
+    std::iota(global_ids, (global_ids + nCells), ZOLTAN_ID_TYPE(myProcOffset));
+
+    // Local indices
+    std::iota(local_ids, (local_ids + nCells), ZOLTAN_ID_TYPE(0));
+
+    // No weights?
+    // Zoltan will assume equally weighted vertices.
+
+    // if (obj_wgts && wgt_dim > 0)
+    // {
+    //     std::fill_n(obj_wgts, wgt_dim, float(1));
+    // }
 }
 
 
-static void get_num_edges_list(void *data, int sizeGID, int sizeLID,
-                      int num_obj,
-             ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID,
-             int *numEdges, int *ierr)
+// [ZOLTAN_NUM_EDGES_MULTI_FN]
+// query function returns the number of edges in the communication
+// graph of the application for each object in a list of objects.
+// That is, for each object in the global_ids/local_ids arrays, the
+// number of objects with which the given object must share
+// information is returned.
+
+static void get_num_edges_list
+(
+    void *data,
+    int sizeGID,
+    int sizeLID,
+    int num_obj,                // Number of graph vertices (mesh cells)
+    ZOLTAN_ID_PTR global_ids,   // (currently unused)
+    ZOLTAN_ID_PTR local_ids,    // rank-local vertex id (mesh cell id)
+    int *numEdges,
+    int *ierr
+)
 {
-  const Foam::polyMesh& mesh = *static_cast<const Foam::polyMesh*>(data);
+    *ierr = ZOLTAN_OK;
+    const auto& mesh = *static_cast<const Foam::polyMesh*>(data);
 
     if ((sizeGID != 1) || (sizeLID != 1) || (num_obj != mesh.nCells()))
     {
@@ -130,32 +171,46 @@ static void get_num_edges_list(void *data, int sizeGID, int sizeLID,
         return;
     }
 
-    for (Foam::label i=0; i < num_obj ;i++)
-    {
-        Foam::label celli = localID[i];
-        const Foam::cell& cFaces = mesh.cells()[celli];
-        forAll(cFaces, cFacei)
-        {
-            Foam::label n = 0;
-            if (mesh.isInternalFace(cFaces[cFacei]))
-            {
-                n++;
-            }
-            numEdges[i] = n;
-        }
-    }
+    const Foam::label nCells = num_obj;
 
-    *ierr = ZOLTAN_OK;
+    for (Foam::label i=0; i < nCells; ++i)
+    {
+        const Foam::label celli = local_ids[i];
+        int numNbr = 0;
+
+        for (const auto facei : mesh.cells()[celli])
+        {
+            if (mesh.isInternalFace(facei))
+            {
+                ++numNbr;
+            }
+            // TBD: check coupled etc
+        }
+
+        numEdges[i] = numNbr;
+    }
 }
 
 
-static void get_edge_list(void *data, int sizeGID, int sizeLID,
-        int num_obj, ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID,
-        int *num_edges,
-        ZOLTAN_ID_PTR nborGID, int *nborProc,
-        int wgt_dim, float *ewgts, int *ierr)
+// [ZOLTAN_EDGE_LIST_MULTI_FN]
+static void get_edge_list
+(
+    void *data,
+    int sizeGID,
+    int sizeLID,
+    int num_obj,
+    ZOLTAN_ID_PTR global_ids,   // (currently unused)
+    ZOLTAN_ID_PTR local_ids,    // rank-local vertex id (mesh cell id)
+    int *num_edges,
+    ZOLTAN_ID_PTR nborGID,
+    int *nborProc,
+    int wgt_dim,
+    float *ewgts,
+    int *ierr
+)
 {
-    const Foam::polyMesh& mesh = *static_cast<const Foam::polyMesh*>(data);
+    *ierr = ZOLTAN_OK;
+    const auto& mesh = *static_cast<const Foam::polyMesh*>(data);
 
     if
     (
@@ -169,20 +224,24 @@ static void get_edge_list(void *data, int sizeGID, int sizeLID,
         return;
     }
 
-    ZOLTAN_ID_TYPE* nextNbor = nborGID;
-    int* nextProc = nborProc;
-    float* nextWgt = ewgts;
+    // Offset for globally unique IDs
+    const auto myProci = Foam::UPstream::myProcNo();
+    const auto myProcOffset =
+        mesh.globalData().globalMeshCellAddr().localStart(myProci);
 
-    for (Foam::label i=0; i < num_obj; i++)
+    auto* nextNbor = nborGID;
+    auto* nextProc = nborProc;
+    auto* nextWgt = ewgts;
+
+    const Foam::label nCells = num_obj;
+
+    for (Foam::label i=0; i < nCells; ++i)
     {
-        Foam::label celli = localID[i];
+        const Foam::label celli = local_ids[i];
+        int numNbr = 0;
 
-        const Foam::cell& cFaces = mesh.cells()[celli];
-        forAll(cFaces, cFacei)
+        for (const auto facei : mesh.cells()[celli])
         {
-            Foam::label n = 0;
-
-            Foam::label facei = cFaces[cFacei];
             if (mesh.isInternalFace(facei))
             {
                 Foam::label nbr = mesh.faceOwner()[facei];
@@ -191,32 +250,37 @@ static void get_edge_list(void *data, int sizeGID, int sizeLID,
                     nbr = mesh.faceNeighbour()[facei];
                 }
 
-                // Note: global index
-                *nextNbor++ = nbr;
-                *nextProc++ = 0;
+                *nextNbor++ = (nbr + myProcOffset); // global id
+                *nextProc++ = myProci;              // rank-local connection
                 *nextWgt++ = 1.0;
 
-                n++;
-            }
-            if (n != num_edges[i])
-            {
-                *ierr = ZOLTAN_FATAL;
-                return;
+                ++numNbr;
             }
         }
+
+        // Sanity check
+        if (numNbr != num_edges[i])
+        {
+            *ierr = ZOLTAN_FATAL;
+            return;
+        }
     }
-    *ierr = ZOLTAN_OK;
 }
 
 
+// [ZOLTAN_NUM_GEOM_FN]
+// The dimensionality of the mesh geometry (2D,3D, etc)
 static int get_mesh_dim(void *data, int *ierr)
 {
-    const Foam::polyMesh& mesh = *static_cast<const Foam::polyMesh*>(data);
+    *ierr = ZOLTAN_OK;
+    const auto& mesh = *static_cast<const Foam::polyMesh*>(data);
 
     return mesh.nSolutionD();
 }
 
 
+// [ZOLTAN_GEOM_MULTI_FN]
+// The geometric location of the graph vertices (mesh cellCentres)
 static void get_geom_list
 (
     void *data,
@@ -225,11 +289,12 @@ static void get_geom_list
     int num_obj,
     ZOLTAN_ID_PTR global_ids,
     ZOLTAN_ID_PTR local_ids,
-    int num_dim,
-    double *geom_vec,
+    int num_dim,                // dimensionality (2|3)
+    double *geom_vec,           // [out] cellCentres
     int *ierr
 )
 {
+    *ierr = ZOLTAN_OK;
     const Foam::polyMesh& mesh = *static_cast<const Foam::polyMesh*>(data);
 
     if
@@ -244,26 +309,46 @@ static void get_geom_list
         return;
     }
 
-    double* p = geom_vec;
-
-
-    const Foam::Vector<Foam::label>& sol = mesh.solutionD();
-
     const Foam::pointField& cc = mesh.cellCentres();
 
-    for (Foam::label celli = 0; celli < num_obj; celli++)
+    if (num_dim == 3)
     {
-        const Foam::point& pt = cc[celli];
+        // Fast path for 3D - assumes that local_ids are still ordered
+        std::copy_n
+        (
+            reinterpret_cast<const Foam::point::cmptType*>(cc.cdata()),
+            (3*num_obj),
+            geom_vec
+        );
+    }
+    else
+    {
+        // [out] cellCentres
+        double* p = geom_vec;
 
-        for (Foam::direction cmpt = 0; cmpt < Foam::vector::nComponents; cmpt++)
+        const Foam::Vector<Foam::label>& sol = mesh.solutionD();
+
+        const Foam::label nCells = num_obj;
+
+        // Assumes that local_ids are still ordered
+        for (Foam::label celli = 0; celli < nCells; ++celli)
         {
-            if (sol[cmpt] == 1)
+            const Foam::point& pt = cc[celli];
+
+            for
+            (
+                Foam::direction cmpt = 0;
+                cmpt < Foam::vector::nComponents;
+                ++cmpt
+            )
             {
-                *p++ = pt[cmpt];
+                if (sol[cmpt] == 1)
+                {
+                    *p++ = pt[cmpt];
+                }
             }
         }
     }
-    *ierr = ZOLTAN_OK;
 }
 
 
@@ -272,7 +357,7 @@ static void get_geom_list
 Foam::zoltanRenumber::zoltanRenumber(const dictionary& dict)
 :
     renumberMethod(dict),
-    coeffsDict_(dict.optionalSubDict(typeName+"Coeffs"))
+    coeffsDict_(dict.optionalSubDict(typeName + "Coeffs"))
 {}
 
 
@@ -284,20 +369,15 @@ Foam::labelList Foam::zoltanRenumber::renumber
     const pointField& points
 ) const
 {
-    stringList args(1);
-    args[0] = "zoltanRenumber";
+    // Zoltan_Initialize will trigger MPI_Init() if not already done.
+    // - use UPstream::initNull() so that OpenFOAM also knows about MPI
 
-    int argc = args.size();
-    char* argv[argc];
-    for (label i = 0; i < argc; i++)
-    {
-        argv[i] = strdup(args[i].c_str());
-    }
+    UPstream::initNull();
+
+    CStringList args({"zoltan-renumber"});
 
     float ver;
-    int rc = Zoltan_Initialize(argc, argv, &ver);
-
-    Foam::Pout<< "Initialised to " << ver << Foam::endl;
+    int rc = Zoltan_Initialize(args.size(), args.strings(), &ver);
 
     if (rc != ZOLTAN_OK)
     {
@@ -307,24 +387,55 @@ Foam::labelList Foam::zoltanRenumber::renumber
 
     struct Zoltan_Struct *zz = Zoltan_Create(MPI_COMM_WORLD);
 
-    polyMesh& mesh = const_cast<polyMesh&>(pMesh);
+    // const bool verbose = coeffsDict_.getOrDefault(verbose, false);
+    const bool verbose = true;
 
+    // Default order method
+    Zoltan_Set_Param(zz, "ORDER_METHOD", "LOCAL_HSFC");
+
+    if (false)
+    {
+        Info<< typeName << " : default ORDER_METHOD = LOCAL_HSFC" << nl
+            << typeName << " : default ORDER_TYPE = LOCAL" << nl;
+    }
 
     for (const entry& dEntry : coeffsDict_)
     {
+        const word& key = dEntry.keyword();
+
+        // Internal keywords
+        if
+        (
+            key == "method"
+         || key == "verbose"
+        )
+        {
+            continue;
+        }
+
         if (!dEntry.isDict())
         {
-            const word& key = dEntry.keyword();
             const word value(dEntry.get<word>());
 
-            Info<< typeName << " : setting parameter " << key
-                << " to " << value << endl;
+            if (verbose)
+            {
+                Info<< typeName
+                    << " : setting parameter "
+                    << key << " = " << value << nl;
+            }
 
             Zoltan_Set_Param(zz, key.c_str(), value.c_str());
         }
     }
 
+    // Always use rank LOCAL ordering
+    Zoltan_Set_Param(zz, "ORDER_TYPE", "LOCAL");
 
+
+    // Set callbacks
+    polyMesh& mesh = const_cast<polyMesh&>(pMesh);
+
+    // Callbacks for graph vertex IDs
     Zoltan_Set_Num_Obj_Fn(zz, get_number_of_vertices, &mesh);
     Zoltan_Set_Obj_List_Fn(zz, get_vertex_list, &mesh);
 
@@ -337,25 +448,26 @@ Foam::labelList Foam::zoltanRenumber::renumber
     Zoltan_Set_Edge_List_Multi_Fn(zz, get_edge_list, &mesh);
 
 
+    const auto nCells = mesh.nCells();
 
-    //Note: !global indices
-    List<ZOLTAN_ID_TYPE> wantedCells(mesh.nCells());
+    List<ZOLTAN_ID_TYPE> globalIds(nCells);
+    List<ZOLTAN_ID_TYPE> oldToNew(nCells);
 
-    globalIndex globalCells(mesh.nCells());
-    forAll(wantedCells, i)
-    {
-        //wantedCells[i] = i;
-        wantedCells[i] = globalCells.toGlobal(i);
-    }
+    // Offset for globally unique IDs
+    const label myProci = UPstream::myProcNo();
+    const label myProcOffset =
+        mesh.globalData().globalMeshCellAddr().localStart(myProci);
 
-    List<ZOLTAN_ID_TYPE> oldToNew(mesh.nCells());
+
+    // Global indices
+    std::iota(globalIds.begin(), globalIds.end(), ZOLTAN_ID_TYPE(myProcOffset));
 
     int err = Zoltan_Order
     (
         zz,
-        1,                                //int num_gid_entries,
-        mesh.globalData().nTotalCells(),  //int num_obj,
-        wantedCells.begin(),
+        1,                          //int num_gid_entries,
+        nCells,                     //int num_obj,
+        globalIds.begin(),
         oldToNew.begin()
     );
 
@@ -365,18 +477,22 @@ Foam::labelList Foam::zoltanRenumber::renumber
             << "Failed Zoltan_Order" << exit(FatalError);
     }
 
-
-    for (label i = 0; i < argc; i++)
-    {
-        free(argv[i]);
-    }
+    Zoltan_Destroy(&zz);
 
 
+    // From global number to rank-local numbering
+    globalIds.clear();
     labelList order(oldToNew.size());
-    forAll(order, i)
-    {
-        order[i] = oldToNew[i];
-    }
+
+    // From global to local (without checks)
+    std::transform
+    (
+        oldToNew.begin(),
+        oldToNew.end(),
+        order.begin(),
+        [=](const ZOLTAN_ID_TYPE id) -> label { return (id - myProcOffset); }
+    );
+
     return order;
 }
 
diff --git a/src/renumber/zoltanRenumber/zoltanRenumber.H b/src/renumber/zoltanRenumber/zoltanRenumber.H
index e2574b057a..b6f7c625ef 100644
--- a/src/renumber/zoltanRenumber/zoltanRenumber.H
+++ b/src/renumber/zoltanRenumber/zoltanRenumber.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2022-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -83,6 +83,10 @@ public:
 
     // Member Functions
 
+        //- Renumbering method requires a polyMesh
+        virtual bool needs_mesh() const { return true; }
+
+
         //- Return the order in which cells need to be visited
         //- (ie. from ordered back to original cell label).
         //  This is only defined for geometric renumberMethods.
diff --git a/tutorials/mesh/parallel/cavity/Allrun b/tutorials/mesh/parallel/cavity/Allrun
index 370935618e..057db8e73e 100755
--- a/tutorials/mesh/parallel/cavity/Allrun
+++ b/tutorials/mesh/parallel/cavity/Allrun
@@ -18,8 +18,7 @@ runParallel -s CuthillMcKee renumberMesh -overwrite $fileHandler
 runParallel -s CuthillMcKee icoFoam $fileHandler
 
 # Bit of bad renumbering and running
-runParallel -s random renumberMesh \
-    -overwrite -dict system/renumberMeshDict-random $fileHandler
+runParallel -s random renumberMesh -renumber-method random -overwrite $fileHandler
 runParallel -s random icoFoam $fileHandler
 
 # Pick up last result
diff --git a/tutorials/mesh/parallel/cavity/system/renumberMeshDict-random b/tutorials/mesh/parallel/cavity/system/renumberMeshDict-random
deleted file mode 100644
index 80b09ee2c3..0000000000
--- a/tutorials/mesh/parallel/cavity/system/renumberMeshDict-random
+++ /dev/null
@@ -1,109 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2312                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      renumberMeshDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Write maps from renumbered back to original mesh
-writeMaps true;
-
-// Optional entry: sort cells on coupled boundaries to last for use with
-// e.g. nonBlockingGaussSeidel.
-sortCoupledFaceCells false;
-
-// Optional entry: renumber on a block-by-block basis. It uses a
-// blockCoeffs dictionary to construct a decompositionMethod to do
-// a block subdivision) and then applies the renumberMethod to each
-// block in turn. This can be used in large cases to keep the blocks
-// fitting in cache with all the cache misses bunched at the end.
-// This number is the approximate size of the blocks - this gets converted
-// to a number of blocks that is the input to the decomposition method.
-//blockSize 1000;
-
-// Optional entry: sort points into internal and boundary points
-//orderPoints false;
-
-// Optional: suppress renumbering cellSets,faceSets,pointSets
-//renumberSets false;
-
-
-//method          CuthillMcKee;
-//method          Sloan;
-//method          manual;
-method          random;
-//method          structured;
-//method          spring;
-//method          zoltan;             // only if compiled with zoltan support
-
-//CuthillMcKeeCoeffs
-//{
-//    // Reverse CuthillMcKee (RCM) or plain
-//    reverse true;
-//}
-
-manualCoeffs
-{
-    // In system directory: new-to-original (i.e. order) labelIOList
-    dataFile "cellMap";
-}
-
-
-// For extruded (i.e. structured in one direction) meshes
-structuredCoeffs
-{
-    // Patches that mesh was extruded from. These determine the starting
-    // layer of cells
-    patches (movingWall);
-    // Method to renumber the starting layer of cells
-    method  random;
-
-    // Renumber in columns (depthFirst) or in layers
-    depthFirst true;
-
-    // Reverse ordering
-    reverse false;
-}
-
-
-springCoeffs
-{
-    // Maximum jump of cell indices. Is fraction of number of cells
-    maxCo 0.01;
-
-    // Limit the amount of movement; the fraction maxCo gets decreased
-    // with every iteration
-    freezeFraction 0.999;
-
-    // Maximum number of iterations
-    maxIter 1000;
-}
-
-
-blockCoeffs
-{
-    method          scotch;
-    //method hierarchical;
-    //hierarchicalCoeffs
-    //{
-    //    n           (1 2 1);
-    //}
-}
-
-
-zoltanCoeffs
-{
-    ORDER_METHOD    LOCAL_HSFC;
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/preProcessing/decompositionConstraints/geometric/Allclean b/tutorials/preProcessing/decompositionConstraints/geometric/Allclean
index b23011b1c2..4afce75b5c 100755
--- a/tutorials/preProcessing/decompositionConstraints/geometric/Allclean
+++ b/tutorials/preProcessing/decompositionConstraints/geometric/Allclean
@@ -6,6 +6,5 @@ cd "${0%/*}" || exit                                # Run from this directory
 cleanCase0
 
 rm -rf constant/triSurface
-rm -f cellDist.vtu
 
 #------------------------------------------------------------------------------
diff --git a/wmake/scripts/have_zoltan b/wmake/scripts/have_zoltan
index b78fcbcdaf..5222fa5365 100644
--- a/wmake/scripts/have_zoltan
+++ b/wmake/scripts/have_zoltan
@@ -5,7 +5,7 @@
 #   \\  /    A nd           | www.openfoam.com
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
-#     Copyright (C) 2018-2020 OpenCFD Ltd.
+#     Copyright (C) 2018-2024 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -64,6 +64,9 @@ search_zoltan()
     local prefix="${1:-system}"
     local header library
 
+    local mpiPrefix="$MPI_ARCH_PATH"
+    local mpiName="${MPI_ARCH_PATH##*/}"
+
     # ----------------------------------
     if isNone "$prefix"
     then
@@ -71,11 +74,29 @@ search_zoltan()
         return 1
     elif hasAbsdir "$prefix"
     then
-        header=$(findFirstFile "$prefix/include/$incName")
-        library=$(findExtLib "$libName")
+        header=$(findFirstFile \
+            "$prefix/include/$FOAM_MPI/$incName" \
+            "$prefix/include/zoltan/$incName" \
+            "$prefix/include/$incName" \
+            "$mpiPrefix/include/$incName" \
+            "$prefix/include/$mpiName/$incName" \
+            "$prefix/include/${mpiName}-$(uname -m)/$incName" \
+        )
+        library=$(findExtLib \
+            "$FOAM_MPI/$libName" \
+            "$libName" \
+        )
     elif isSystem "$prefix"
     then
-        header=$(findSystemInclude -name="$incName")
+        header=$(findFirstFile \
+            "/usr/local/include/zoltan/$incName" \
+            "/usr/local/include/$incName" \
+            "/usr/include/zoltan/$incName" \
+            "/usr/include/$incName" \
+            "$mpiPrefix/include/$incName" \
+            "/usr/include/$mpiName/$incName" \
+            "$prefix/include/${mpiName}-$(uname -m)/$incName" \
+        )
         prefix=$(sysPrefix "$header")
     else
         unset prefix
@@ -91,6 +112,7 @@ search_zoltan()
     # Library
     [ -n "$library" ] \
     || library=$(findLibrary -prefix="$prefix" -name="$libName") \
+    || library=$(findLibrary -prefix="$mpiPrefix" -name="$libName") \
     || {
         [ -n "$warn" ] && echo "$warn (no library)"
         return 2
@@ -99,6 +121,7 @@ search_zoltan()
     # ----------------------------------
 
     # OK
+    # echo "zoltan - $prefix"
     export HAVE_ZOLTAN=true
     export ZOLTAN_ARCH_PATH="$prefix"
     export ZOLTAN_INC_DIR="${header%/*}"     # Basename