WIP: SUBMODULE update

- previously read the cloud directly without any of the
  passiveParticle fields (origProcId, origId), which meant they would
  not actually be converted.

.gitlab/issue_templates/bug.md

@@ -49,7 +49,7 @@
 
 <!--
   Providing details of your set-up can help us identify any issues, e.g.
-  OpenFOAM version : v2306|v2212|v2206|v2112|v2106 etc
+  OpenFOAM version : v2312|v2306|v2212|v2206|v2112 etc
   Operating system : ubuntu|openSUSE|centos etc
   Hardware info    : any info that may help?
   Compiler         : gcc|intel|clang etc

META-INFO/api-info

@@ -1,2 +1,2 @@
-api=2309
+api=2310
 patch=0

README.md

@@ -40,9 +40,9 @@ Violations of the Trademark are monitored, and will be duly prosecuted.
 
 If OpenFOAM has already been compiled on your system, simply source
 the appropriate `etc/bashrc` or `etc/cshrc` file and get started.
-For example, for the OpenFOAM-v2306 version:
+For example, for the OpenFOAM-v2312 version:
 ```
-source /installation/path/OpenFOAM-v2306/etc/bashrc
+source /installation/path/OpenFOAM-v2312/etc/bashrc
 ```
 
 ## Compiling OpenFOAM
@@ -127,8 +127,8 @@ These 3rd-party sources are normally located in a directory parallel
 to the OpenFOAM directory. For example,
 ```
 /path/parent
-|-- OpenFOAM-v2306
+|-- OpenFOAM-v2312
-\-- ThirdParty-v2306
+\-- ThirdParty-v2312
 ```
 There are, however, many cases where this simple convention is inadequate:
 
@@ -136,7 +136,7 @@ There are, however, many cases where this simple convention is inadequate:
   operating system or cluster installation provides it)
 
 * When we have changed the OpenFOAM directory name to some arbitrary
-  directory name, e.g. openfoam-sandbox2306, etc..
+  directory name, e.g. openfoam-sandbox2312, etc..
 
 * When we would like any additional 3rd party software to be located
   inside of the OpenFOAM directory to ensure that the installation is
@@ -156,9 +156,9 @@ when locating the ThirdParty directory with the following precedence:
 2. PREFIX/ThirdParty-VERSION
    * this corresponds to the traditional approach
 3. PREFIX/ThirdParty-vAPI
-   * allows for an updated value of VERSION, *eg*, `v2306-myCustom`,
+   * allows for an updated value of VERSION, *eg*, `v2312-myCustom`,
      without requiring a renamed ThirdParty. The API value would still
-     be `2306` and the original `ThirdParty-v2306/` would be found.
+     be `2312` and the original `ThirdParty-v2312/` would be found.
 4. PREFIX/ThirdParty-API
    * same as the previous example, but using an unadorned API value.
 5. PREFIX/ThirdParty-common
@@ -213,4 +213,4 @@ ThirdParty directory will contain either an `Allwmake` file or a
 - [Governance](http://www.openfoam.com/governance/), [Governance Projects](https://www.openfoam.com/governance/projects)
 - [Contacting OpenCFD](http://www.openfoam.com/contact/)
 
-Copyright 2016-2022 OpenCFD Ltd
+Copyright 2016-2023 OpenCFD Ltd

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C

@@ -329,7 +329,7 @@ Foam::radiation::laserDTRM::laserDTRM(const volScalarField& T)
 :
     radiationModel(typeName, T),
     mode_(powerDistNames_.get("mode", *this)),
-    DTRMCloud_(mesh_, "DTRMCloud", IDLList<DTRMParticle>()),
+    DTRMCloud_(mesh_, Foam::zero{}, "DTRMCloud"),  // Empty cloud
     nParticles_(0),
     ndTheta_(get<label>("nTheta")),
     ndr_(get<label>("nr")),
@@ -427,7 +427,7 @@ Foam::radiation::laserDTRM::laserDTRM
 :
     radiationModel(typeName, dict, T),
     mode_(powerDistNames_.get("mode", *this)),
-    DTRMCloud_(mesh_, "DTRMCloud", IDLList<DTRMParticle>()),
+    DTRMCloud_(mesh_, Foam::zero{}, "DTRMCloud"),  // Empty cloud
     nParticles_(0),
     ndTheta_(get<label>("nTheta")),
     ndr_(get<label>("nr")),

applications/test/objectRegistry2/Test-objectRegistry2.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,10 +45,31 @@ Description
 
 using namespace Foam;
 
+enum loadTestTypes
+{
+    PLAIN_PTR, AUTO_PTR, REF_PTR, TMP_PTR, CACHE_PTR
+};
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+void printIOobject(const regIOobject* io)
+{
+    if (io)
+    {
+        Info<< io->name() << ' ' << "type=" << io->type()
+            << " registered=" << io->registered()
+            << " owned=" << io->ownedByRegistry() << endl;
+    }
+}
+
+
 template<class Type>
-bool loadField(fvMesh& mesh, const word& fieldName)
+bool loadField
+(
+    fvMesh& mesh,
+    const word& fieldName,
+    enum loadTestTypes wrapper = loadTestTypes::PLAIN_PTR
+)
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
@@ -66,20 +87,67 @@ bool loadField(fvMesh& mesh, const word& fieldName)
         mesh,
         IOobject::MUST_READ,
         IOobject::NO_WRITE
+        // Value of register is fairly irrelevant
     );
 
     if (fieldHeader.typeHeaderOk<VolFieldType>(true, true, false))
     {
         // Store field on mesh database
-        VolFieldType* ptr = new VolFieldType(fieldHeader, mesh);
+        switch (wrapper)
-        mesh.objectRegistry::store(ptr);
+        {
+            case loadTestTypes::PLAIN_PTR :
+            {
+                auto* ptr = new VolFieldType(fieldHeader, mesh);
+                printIOobject(ptr);
+                regIOobject::store(ptr);
+                break;
+            }
+
+            case loadTestTypes::AUTO_PTR :
+            {
+                auto ptr = autoPtr<VolFieldType>::New(fieldHeader, mesh);
+                printIOobject(ptr.get());
+                regIOobject::store(ptr);
+                break;
+            }
+
+            case loadTestTypes::REF_PTR :
+            {
+                auto ptr = refPtr<VolFieldType>::New(fieldHeader, mesh);
+                printIOobject(ptr.get());
+                regIOobject::store(ptr);
+                break;
+            }
+
+            case loadTestTypes::TMP_PTR :
+            {
+                auto ptr = tmp<VolFieldType>::New(fieldHeader, mesh);
+                printIOobject(ptr.get());
+                Info<< "pointer:" << ptr.is_pointer()
+                    << " movable:" << ptr.movable() << nl;
+                regIOobject::store(ptr);
+                break;
+            }
+
+            case loadTestTypes::CACHE_PTR :
+            {
+                auto ptr = tmp<VolFieldType>::New(fieldHeader, mesh);
+                ptr.protect(true);
+                printIOobject(ptr.get());
+                Info<< "pointer:" << ptr.is_pointer()
+                    << " movable:" << ptr.movable() << nl;
+                regIOobject::store(ptr);
+                break;
+            }
+        }
+
         return true;
     }
     else if (fieldHeader.typeHeaderOk<SurfaceFieldType>(true, true, false))
     {
         // Store field on mesh database
         SurfaceFieldType* ptr = new SurfaceFieldType(fieldHeader, mesh);
-        mesh.objectRegistry::store(ptr);
+        regIOobject::store(ptr);
         return true;
     }
 
@@ -87,28 +155,38 @@ bool loadField(fvMesh& mesh, const word& fieldName)
 }
 
 
-bool loadField(fvMesh& mesh, const word& fieldName)
+bool loadField
+(
+    fvMesh& mesh,
+    const word& fieldName,
+    enum loadTestTypes wrapper = loadTestTypes::PLAIN_PTR
+)
 {
     return
     (
         !mesh.objectRegistry::found(fieldName)
     &&
         (
-            loadField<scalar>(mesh, fieldName)
+            loadField<scalar>(mesh, fieldName, wrapper)
-         || loadField<vector>(mesh, fieldName)
+         || loadField<vector>(mesh, fieldName, wrapper)
-         || loadField<sphericalTensor>(mesh, fieldName)
+         || loadField<sphericalTensor>(mesh, fieldName, wrapper)
-         || loadField<symmTensor>(mesh, fieldName)
+         || loadField<symmTensor>(mesh, fieldName, wrapper)
-         || loadField<tensor>(mesh, fieldName)
+         || loadField<tensor>(mesh, fieldName, wrapper)
         )
     );
 }
 
 
-void loadFields(fvMesh& mesh, const IOobjectList& objects)
+void loadFields
+(
+    fvMesh& mesh,
+    const IOobjectList& objects,
+    enum loadTestTypes wrapper = loadTestTypes::PLAIN_PTR
+)
 {
     for (const word& fieldName : objects.names())
     {
-        loadField(mesh, fieldName);
+        loadField(mesh, fieldName, wrapper);
     }
 }
 
@@ -286,8 +364,46 @@ int main(int argc, char *argv[])
     // timeSelector::addOptions();
     timeSelector::addOptions(true, true);
 
+    argList::addBoolOption("plain", "Store from raw pointer (default)");
+    argList::addBoolOption("autoPtr", "Store from autoPtr");
+    argList::addBoolOption("refPtr", "Store from refPtr");
+    argList::addBoolOption("cacheTmp", "Store from tmp (cached)");
+    argList::addBoolOption("tmp",    "Store from tmp (regular)");
+
+    argList::addVerboseOption("increase debug value");
+
     #include "setRootCase.H"
 
+    const int optVerbose = args.verbose();
+
+    enum loadTestTypes loadWrapper = loadTestTypes::PLAIN_PTR;
+
+    if (args.found("autoPtr"))
+    {
+        Info<< "loading via autoPtr" << nl;
+        loadWrapper = loadTestTypes::AUTO_PTR;
+    }
+    else if (args.found("refPtr"))
+    {
+        Info<< "loading via refPtr" << nl;
+        loadWrapper = loadTestTypes::REF_PTR;
+    }
+    else if (args.found("cacheTmp"))
+    {
+        Info<< "loading via tmp (cached)" << nl;
+        loadWrapper = loadTestTypes::CACHE_PTR;
+    }
+    else if (args.found("tmp"))
+    {
+        Info<< "loading via tmp (regular)" << nl;
+        loadWrapper = loadTestTypes::TMP_PTR;
+    }
+    else
+    {
+        Info<< "loading via plain ptr" << nl;
+        loadWrapper = loadTestTypes::PLAIN_PTR;
+    }
+
 //    wordRes matcher;
 //    if (args.readListIfPresent<wordRe>("filter", matcher))
 //    {
@@ -299,6 +415,12 @@ int main(int argc, char *argv[])
 
     instantList timeDirs = timeSelector::select0(runTime, args);
 
+    if (optVerbose)
+    {
+        objectRegistry::debug = optVerbose;
+        regIOobject::debug = optVerbose;
+    }
+
     forAll(timeDirs, timeI)
     {
         runTime.setTime(timeDirs[timeI], timeI);
@@ -309,7 +431,7 @@ int main(int argc, char *argv[])
         IOobjectList objects(mesh, runTime.timeName());
 
         // Read volFields
-        loadFields(mesh, objects);
+        loadFields(mesh, objects, loadWrapper);
 
         printRegistry(Info, mesh);
 

applications/test/parallel-comm1/oneBlock_cyclic/0/T

@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    top
+    {
+        type            cyclic; //AMI;
+        //useImplicit     true;
+        value           $internalField;
+    }
+
+    bottom
+    {
+        type            cyclic; //AMI;
+        //useImplicit     true;
+        value           $internalField;
+    }
+
+    left
+    {
+        type            fixedValue;
+        value           uniform 1;
+    }
+
+    right
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

applications/test/parallel-comm1/oneBlock_cyclic/constant/transportProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+DT              4e-05;
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/oneBlock_cyclic/system/blockMeshDict

@@ -0,0 +1,93 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   0.1;
+
+vertices
+(
+    (0 0 0)
+    (1 0 0)
+    (1 1 0)
+    (0 1 0)
+    (0 0 0.1)
+    (1 0 0.1)
+    (1 1 0.1)
+    (0 1 0.1)
+);
+
+blocks
+(
+    //- Left block
+    hex (0 1 2 3 4 5 6 7) (60 50 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    top
+    {
+        type            cyclic; //AMI;
+        neighbourPatch  bottom;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    bottom
+    {
+        type            cyclic; //AMI;
+        neighbourPatch  top;
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+    left
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+
+    right
+    {
+        type patch;
+        faces
+        (
+            (1 2 5 6)
+        );
+    }
+
+    frontAndBack
+    {
+        type empty;
+        faces
+        (
+            (0 3 2 1)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

applications/test/parallel-comm1/oneBlock_cyclic/system/controlDict

@@ -0,0 +1,82 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+libs        (utilityFunctionObjects);
+
+DebugSwitches
+{
+//    mappedPatchBase 2;
+//    syncObjects     2;
+
+    disallowGenericFvPatchField 1;
+    disallowGenericFvsPatchField 1;
+    disallowGenericPointPatchField 1;
+    disallowGenericPolyPatch    1;
+    lduPrimitiveMesh            1;
+    GAMGAgglomeration           2;
+    GAMG                        2;
+}
+
+
+application     laplacianFoam;
+
+startFrom       startTime;  //latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1;
+
+deltaT          1;
+
+//writeControl    runTime;
+//writeInterval   0.1;
+writeControl    timeStep;
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+//functions
+//{
+//    syncObjects
+//    {
+//        type        syncObjects;
+//        libs        (utilityFunctionObjects);
+//
+//        // Where is data located relative to runTime. Given as a filename
+//        // with every '/' indicating a sub-objectRegistry w.r.t. runTime.
+//        // Local data is under <root>/send/processorXXX. After execution
+//        // data will be under the corresponding <root>/receive/processorYYY
+//        // objectRegistry
+//        //root        "level0/level1/level2";
+//    }
+//}
+
+// ************************************************************************* //

applications/test/parallel-comm1/oneBlock_cyclic/system/decomposeParDict

@@ -0,0 +1,35 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- The total number of domains (mandatory)
+numberOfSubdomains  20;
+
+//- The decomposition method (mandatory)
+method          scotch;
+
+constraints
+{
+    patches
+    {
+        //- Keep owner and neighbour on same processor for faces in patches
+        //  (only makes sense for cyclic patches and cyclicAMI)
+        type    preservePatches;
+        patches (cyclic);
+    }
+}
+
+// ************************************************************************* //

applications/test/parallel-comm1/oneBlock_cyclic/system/fvSchemes

@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;    //Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(T)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+}
+
+laplacianSchemes
+{
+    default         none;
+    laplacian(DT,T) Gauss linear corrected;
+    laplacian(DTV,T) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/oneBlock_cyclic/system/fvSolution

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    T
+    {
+        //solver          PCG;
+        //preconditioner  DIC;
+        solver          GAMG;
+        smoother        GaussSeidel;
+        nCellsInCoarsestLevel   1;
+        processorAgglomerator   masterCoarsest;
+
+        tolerance       1e-06;
+        relTol          0;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+}
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/0/T

@@ -0,0 +1,65 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    AMI1
+    {
+        type            cyclicAMI;
+        //useImplicit     true;
+        value           $internalField;
+    }
+
+    AMI2
+    {
+        type            cyclicAMI;
+        //useImplicit     true;
+        value           $internalField;
+    }
+
+    top
+    {
+        type            zeroGradient;
+    }
+
+    bottom
+    {
+        type            zeroGradient;
+    }
+
+    left
+    {
+        type            fixedValue;
+        value           uniform 1;
+    }
+
+    right
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/0/cellToRegion

@@ -0,0 +1,155 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2208                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       labelList;
+    location    "0";
+    object      cellToRegion;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+130
+(
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+1
+)
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/0/processorID

@@ -0,0 +1,62 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2208                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       volScalarField;
+    location    "0";
+    object      processorID;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    AMI1
+    {
+        type            cyclicAMI;
+        value           uniform 0;
+    }
+    AMI2
+    {
+        type            cyclicAMI;
+        value           uniform 0;
+    }
+    top
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    bottom
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    left
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    right
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/constant/transportProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+DT              4e-05;
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/system/blockMeshDict

@@ -0,0 +1,136 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   0.1;
+
+vertices
+(
+    (0 0 0)
+    (0.5 0 0)
+    (0.5 1 0)
+    (0 1 0)
+    (0 0 0.1)
+    (0.5 0 0.1)
+    (0.5 1 0.1)
+    (0 1 0.1)
+
+
+    (0.5 0 0)
+    (1 0 0)
+    (1 1 0)
+    (0.5 1 0)
+    (0.5 0 0.1)
+    (1 0 0.1)
+    (1 1 0.1)
+    (0.5 1 0.1)
+
+);
+
+blocks
+(
+    //- Left block
+    hex (0 1 2 3 4 5 6 7) left (5 2 1) simpleGrading (1 1 1)
+
+    //- Right block
+    hex (8 9 10 11 12 13 14 15) right (5 2 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    AMI1
+    {
+        type            cyclicAMI;
+        neighbourPatch  AMI2;
+        transform       noOrdering;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+
+    AMI2
+    {
+        type            cyclicAMI;
+        neighbourPatch  AMI1;
+        transform       noOrdering;
+        faces
+        (
+            (8 12 15 11)
+        );
+    }
+
+    top
+    {
+        type patch;
+        faces
+        (
+            (3 7 6 2)
+
+            (11 15 14 10)
+        );
+    }
+    bottom
+    {
+        type patch;
+        faces
+        (
+            (1 5 4 0)
+
+            (9 13 12 8)
+        );
+    }
+    left
+    {
+        type patch;
+        faces
+        (
+            (0 4 7 3)
+
+        );
+    }
+
+    right
+    {
+        type patch;
+        faces
+        (
+            (10 14 13 9)
+        );
+    }
+
+    frontAndBack
+    {
+        type empty;
+        faces
+        (
+            (0 3 2 1)
+            (4 5 6 7)
+
+            (8 11 10 9)
+            (12 13 14 15)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/system/controlDict

@@ -0,0 +1,82 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+libs        (utilityFunctionObjects);
+
+DebugSwitches
+{
+//    mappedPatchBase 2;
+//    syncObjects     2;
+
+    disallowGenericFvPatchField 1;
+    disallowGenericFvsPatchField 1;
+    disallowGenericPointPatchField 1;
+    disallowGenericPolyPatch    1;
+    lduPrimitiveMesh            1;
+    GAMGAgglomeration           2;
+    GAMG                        2;
+}
+
+
+application     laplacianFoam;
+
+startFrom       startTime;  //latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1;
+
+deltaT          1;
+
+//writeControl    runTime;
+//writeInterval   0.1;
+writeControl    timeStep;
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+//functions
+//{
+//    syncObjects
+//    {
+//        type        syncObjects;
+//        libs        (utilityFunctionObjects);
+//
+//        // Where is data located relative to runTime. Given as a filename
+//        // with every '/' indicating a sub-objectRegistry w.r.t. runTime.
+//        // Local data is under <root>/send/processorXXX. After execution
+//        // data will be under the corresponding <root>/receive/processorYYY
+//        // objectRegistry
+//        //root        "level0/level1/level2";
+//    }
+//}
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/system/decomposeParDict

@@ -0,0 +1,26 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- The total number of domains (mandatory)
+numberOfSubdomains  3;
+
+//- The decomposition method (mandatory)
+method          hierarchical;
+n               (3 1 1);
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/system/fvSchemes

@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;    //Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(T)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+}
+
+laplacianSchemes
+{
+    default         none;
+    laplacian(DT,T) Gauss linear corrected;
+    laplacian(DTV,T) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //

applications/test/parallel-comm1/twoBlocks_cyclicAMI/system/fvSolution

@@ -0,0 +1,40 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    T
+    {
+        //solver          PCG;
+        //preconditioner  DIC;
+        solver          GAMG;
+        smoother        GaussSeidel;
+        nCellsInCoarsestLevel   1;
+        processorAgglomerator   masterCoarsest;
+
+        tolerance       1e-06;
+        relTol          0;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+}
+
+
+// ************************************************************************* //

applications/test/passiveParticle/Make/files

@@ -1,3 +1,3 @@
-Test-passiveParticle.C
+Test-passiveParticle.cxx
 
 EXE = $(FOAM_USER_APPBIN)/Test-passiveParticle

applications/test/passiveParticle/Test-passiveParticle.cxx

@@ -24,7 +24,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    testPassiveParticle
+    Test-passiveParticle
 
 Description
     Test cloud of passive particles.
@@ -50,12 +50,8 @@ int main(int argc, char *argv[])
 
     {
         // Start with empty cloud
-        passiveParticleCloud particles
+        passiveParticleCloud particles(mesh, Foam::zero{}, cloudName);
-        (
+
-            mesh,
-            cloudName,
-            IDLList<passiveParticle>()
-        );
         Pout<< "Starting particles:" << particles.size() << endl;
 
         Pout<< "Adding a particle." << endl;

applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C

@@ -73,8 +73,7 @@ Usage
 #include "tensorIOField.H"
 #include "labelFieldIOField.H"
 #include "vectorFieldIOField.H"
-#include "Cloud.H"
+#include "passiveParticleCloud.H"
-#include "passiveParticle.H"
 #include "fieldDictionary.H"
 
 #include "writeMeshObject.H"
@@ -95,7 +94,7 @@ bool writeZones
     const word& name,
     const fileName& meshDir,
     Time& runTime,
-    const IOstreamOption::compressionType compression
+    IOstreamOption::compressionType compression
 )
 {
     IOobject io
@@ -123,11 +122,11 @@ bool writeZones
 
         IOPtrList<entry> meshObject(io);
 
-        forAll(meshObject, i)
+        for (entry& e : meshObject)
         {
-            if (meshObject[i].isDict())
+            if (e.isDict())
             {
-                dictionary& d = meshObject[i].dict();
+                dictionary& d = e.dict();
 
                 if (d.found("faceLabels"))
                 {
@@ -197,13 +196,13 @@ struct uniqueEqOp
 };
 
 
-template<class T>
+template<class Type>
 bool writeOptionalMeshObject
 (
     const word& name,
     const fileName& meshDir,
     Time& runTime,
-    const bool writeOnProc
+    bool writeOnProc
 )
 {
     IOobject io
@@ -217,22 +216,26 @@ bool writeOptionalMeshObject
         IOobject::NO_REGISTER
     );
 
+    // Check if available and the correct type.
+    // Done as typeHeaderOk<regIOobject> + isHeaderClass to ensure
+    // local-only reading and circumvent is_globalIOobject<Type> check
+    // in typeHeaderOk<Type>
+
+    bool readOnProc =
+    (
+        io.typeHeaderOk<regIOobject>(false)
+     && io.isHeaderClass<Type>()
+    );
+
     bool writeOk = false;
-    const bool haveFile = io.typeHeaderOk<regIOobject>(false);
 
-    // Make sure all know if there is a valid class name
+    if (returnReduceOr(readOnProc))
-    wordList classNames(1, io.headerClassName());
-    Pstream::combineReduce(classNames, uniqueEqOp<word>());
-
-    // Check for correct type
-    if (classNames[0] == T::typeName)
     {
-        Info<< "        Reading " << classNames[0]
+        Info<< "        Reading " << Type::typeName << " : " << name << endl;
-            << " : " << name << endl;
+        Type meshObject(io, readOnProc && writeOnProc);
-        T meshObject(io, writeOnProc && haveFile);
 
         Info<< "        Writing " << name << endl;
-        writeOk = meshObject.regIOobject::write(writeOnProc && haveFile);
+        writeOk = meshObject.regIOobject::write(readOnProc && writeOnProc);
     }
 
     return writeOk;
@@ -416,7 +419,6 @@ int main(int argc, char *argv[])
         }
 
 
-
         // Check for lagrangian
         fileNameList lagrangianDirs
         (
@@ -451,12 +453,14 @@ int main(int argc, char *argv[])
                             polyMesh::defaultRegion,
                             runTime.timeName(),
                             runTime,
-                            Foam::IOobject::MUST_READ
+                            IOobject::MUST_READ
                         )
                     )
                 );
             }
 
+            const auto& mesh = meshPtr();
+
             fileNameList cloudDirs
             (
                 fileHandler().readDir
@@ -468,11 +472,14 @@ int main(int argc, char *argv[])
 
             Pstream::combineReduce(cloudDirs, uniqueEqOp<fileName>());
 
-            forAll(cloudDirs, i)
+            for (const auto& cloudDir : cloudDirs)
             {
-                fileName dir(cloud::prefix/cloudDirs[i]);
+                fileName dir(cloud::prefix/cloudDir);
 
-                Cloud<passiveParticle> parcels(meshPtr(), cloudDirs[i], false);
+                // Read with origProcId,origId fields
+                passiveParticleCloud parcels(mesh, cloudDir, true);
+
+                const bool writeOnProc = parcels.size();
 
                 parcels.writeObject
                 (
@@ -481,7 +488,7 @@ int main(int argc, char *argv[])
                         runTime.writeFormat(),
                         runTime.writeCompression()
                     ),
-                    parcels.size()
+                    writeOnProc
                 );
 
 
@@ -496,10 +503,7 @@ int main(int argc, char *argv[])
 
                 for (const word& name : cloudFields)
                 {
-                    // Note: try the various field types. Make sure to
+                    // These ones already done by cloud itself
-                    //       exit once successful conversion to avoid re-read
-                    //       converted file.
-
                     if
                     (
                         name == "positions"
@@ -511,81 +515,32 @@ int main(int argc, char *argv[])
                         continue;
                     }
 
-                    bool writeOk = writeOptionalMeshObject<labelIOField>
+                    #undef  doLocalCode
-                    (
+                    #define doLocalCode(Type)                                 \
-                        name,
+                    if                                                        \
-                        dir,
+                    (                                                         \
-                        runTime,
+                        writeOptionalMeshObject<Type>                         \
-                        parcels.size() > 0
+                        (                                                     \
-                    );
+                            name, dir, runTime, writeOnProc                   \
-                    if (writeOk) continue;
+                        )                                                     \
-
+                    )                                                         \
-                    writeOk = writeOptionalMeshObject<scalarIOField>
+                    {                                                         \
-                    (
+                        continue;                                             \
-                        name,
-                        dir,
-                        runTime,
-                        parcels.size() > 0
-                    );
-                    if (writeOk) continue;
-
-                    writeOk = writeOptionalMeshObject<vectorIOField>
-                    (
-                        name,
-                        dir,
-                        runTime,
-                        parcels.size() > 0
-                    );
-                    if (writeOk) continue;
-
-                    writeOk = writeOptionalMeshObject<sphericalTensorIOField>
-                    (
-                        name,
-                        dir,
-                        runTime,
-                        parcels.size() > 0
-                    );
-                    if (writeOk) continue;
-
-                    writeOk = writeOptionalMeshObject<symmTensorIOField>
-                    (
-                        name,
-                        dir,
-                        runTime,
-                        parcels.size() > 0
-                    );
-                    if (writeOk) continue;
-
-                    writeOk = writeOptionalMeshObject<tensorIOField>
-                    (
-                        name,
-                        dir,
-                        runTime,
-                        parcels.size() > 0
-                    );
-                    if (writeOk) continue;
-
-                    writeOk = writeOptionalMeshObject<labelFieldIOField>
-                    (
-                        name,
-                        dir,
-                        runTime,
-                        parcels.size() > 0
-                    );
-                    if (writeOk) continue;
-
-                    writeOk = writeOptionalMeshObject<vectorFieldIOField>
-                    (
-                        name,
-                        dir,
-                        runTime,
-                        parcels.size() > 0
-                    );
-
-                    if (!writeOk)
-                    {
-                        Info<< "        Failed converting " << name << endl;
                     }
+
+                    doLocalCode(labelIOField);
+                    doLocalCode(scalarIOField);
+                    doLocalCode(vectorIOField);
+                    doLocalCode(sphericalTensorIOField);
+                    doLocalCode(symmTensorIOField);
+                    doLocalCode(tensorIOField);
+
+                    doLocalCode(labelFieldIOField);
+                    doLocalCode(vectorFieldIOField);
+
+                    #undef doLocalCode
+
+                    Info<< "        Failed converting " << name << endl;
                 }
             }
         }

applications/utilities/miscellaneous/foamFormatConvert/writeMeshObject.H

@@ -71,10 +71,9 @@ inline bool writeMeshObject
 
         // Switch off type checking (for reading e.g. faceZones as
         // generic list of dictionaries).
-        word oldTypeName;
+        const word oldTypeName = Type::typeName;
         if (disableHeaderChecking)
         {
-            oldTypeName = Type::typeName;
             const_cast<word&>(Type::typeName) = word::null;
         }
 

applications/utilities/parallelProcessing/redistributePar/parLagrangianDistributor.C

@@ -227,19 +227,21 @@ Foam::parLagrangianDistributor::distributeLagrangianPositions
     pBufs.finishedSends();
 
 
+    // The cloud name
+    const word cloudName = lpi.name();
+
     {
         // Temporarily rename original cloud so we can construct a new one
         // (to distribute the positions) without getting a duplicate
         // registration warning
-        const word cloudName = lpi.name();
         lpi.rename(cloudName + "_old");
 
-        // New cloud on tgtMesh
+        // New empty cloud on tgtMesh
         passivePositionParticleCloud lagrangianPositions
         (
             tgtMesh_,
-            cloudName,
+            Foam::zero{},
-            IDLList<passivePositionParticle>()
+            cloudName
         );
 
         // Retrieve from receive buffers
@@ -310,11 +312,11 @@ Foam::parLagrangianDistributor::distributeLagrangianPositions
         //        ).objectPath()
         //    );
         //}
-
-        // Restore cloud name
-        lpi.rename(cloudName);
     }
 
+    // Restore cloud name
+    lpi.rename(cloudName);
+
 
     // The constructMap is in linear (processor) order
     return autoPtr<mapDistributeBase>::New

applications/utilities/parallelProcessing/redistributePar/passivePositionParticleCloud.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,6 +34,7 @@ namespace Foam
     defineTemplateTypeNameAndDebug(Cloud<passivePositionParticle>, 0);
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::passivePositionParticleCloud::passivePositionParticleCloud
@@ -52,15 +53,4 @@ Foam::passivePositionParticleCloud::passivePositionParticleCloud
 }
 
 
-Foam::passivePositionParticleCloud::passivePositionParticleCloud
-(
-    const polyMesh& mesh,
-    const word& cloudName,
-    const IDLList<passivePositionParticle>& particles
-)
-:
-    Cloud<passivePositionParticle>(mesh, cloudName, particles)
-{}
-
-
 // ************************************************************************* //

applications/utilities/parallelProcessing/redistributePar/passivePositionParticleCloud.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2019 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,8 +34,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef passivePositionParticleCloud_H
+#ifndef Foam_passivePositionParticleCloud_H
-#define passivePositionParticleCloud_H
+#define Foam_passivePositionParticleCloud_H
 
 #include "Cloud.H"
 #include "passivePositionParticle.H"
@@ -53,23 +53,21 @@ class passivePositionParticleCloud
 :
     public Cloud<passivePositionParticle>
 {
-    // Private Member Functions
-
-        //- No copy construct
-        passivePositionParticleCloud
-        (
-            const passivePositionParticleCloud&
-        ) = delete;
-
-        //- No copy assignment
-        void operator=(const passivePositionParticleCloud&) = delete;
-
-
 public:
 
+    //- Type of parcel within the cloud
+    typedef passivePositionParticle parcelType;
+
+    //- No copy construct
+    passivePositionParticleCloud(const passivePositionParticleCloud&) = delete;
+
+    //- No copy assignment
+    void operator=(const passivePositionParticleCloud&) = delete;
+
+
     // Constructors
 
-        //- Construct given mesh
+        //- Read construct
         explicit passivePositionParticleCloud
         (
             const polyMesh& mesh,
@@ -77,13 +75,16 @@ public:
             bool readFields = true
         );
 
-        //- Construct from mesh, cloud name, and a list of particles
+        //- Construct without particles
         passivePositionParticleCloud
         (
             const polyMesh& mesh,
-            const word& cloudName,
+            const Foam::zero,
-            const IDLList<passivePositionParticle>& particles
+            const word& cloudName = cloud::defaultName
-        );
+        )
+        :
+            Cloud<passivePositionParticle>(mesh, Foam::zero{}, cloudName)
+        {}
 };
 
 

applications/utilities/parallelProcessing/redistributePar/redistributePar.C

@@ -112,7 +112,7 @@ using namespace Foam;
 
 // Use -verbose -verbose to switch on debug info. TBD.
 int debug(::Foam::debug::debugSwitch("redistributePar", 0));
-#define InfoOrPout (::debug ? Pout : Info())
+#define InfoOrPout (::debug ? Pout : Info.stream())
 
 
 // Allocate a new file handler on valid processors only
@@ -1216,7 +1216,10 @@ int main(int argc, char *argv[])
     #include "addOverwriteOption.H"
     argList::addBoolOption("decompose", "Decompose case");
     argList::addBoolOption("reconstruct", "Reconstruct case");
-    argList::addVerboseOption("Additional verbosity");
+    argList::addVerboseOption
+    (
+        "Additional verbosity. (Can be used multiple times for debugging)"
+    );
     argList::addDryRunOption
     (
         "Test without writing the decomposition. "
@@ -1267,9 +1270,9 @@ int main(int argc, char *argv[])
     // - reconstruct - reads parallel, write on master only and to parent
     //                 directory
 
-    // Detect if running data-distributed (processors cannot see all other
+    // Running data-distributed?
-    // processors' files)
+    // (processors cannot see all other processors' files)
-    const bool nfs = !fileHandler().distributed();
+    const bool hasDistributedFiles(fileHandler().distributed());
 
 
     // Set up loose processorsXXX directory matching (for collated) so e.g.
@@ -1278,21 +1281,6 @@ int main(int argc, char *argv[])
     // matching.
     fileOperation::nProcsFilter(0);
 
-
-    // Read handler on processors with a volMesh
-    refPtr<fileOperation> volMeshReadHandler;
-
-    // Read handler on processors with an areaMesh
-    refPtr<fileOperation> areaMeshReadHandler;
-
-    // Handler for master-only operation (read/writing from/to undecomposed)
-    refPtr<fileOperation> masterOnlyHandler;
-    if (UPstream::master(UPstream::worldComm))
-    {
-        masterOnlyHandler = fileOperation::NewUncollated();
-    }
-
-
     // Need this line since we don't include "setRootCase.H"
     #include "foamDlOpenLibs.H"
 
@@ -1306,19 +1294,36 @@ int main(int argc, char *argv[])
     bool decompose = args.found("decompose");
     bool overwrite = args.found("overwrite");
 
+    // Field restrictions...
+    const wordRes selectedFields;
+    const wordRes selectedLagrangianFields;
+
+
+    if (!UPstream::parRun())
+    {
+        FatalErrorInFunction
+            << ": This utility can only be run parallel"
+            << exit(FatalError);
+    }
+
+    if (decompose && reconstruct)
+    {
+        FatalErrorInFunction
+            << "Cannot specify both -decompose and -reconstruct"
+            << exit(FatalError);
+    }
+
     if (optVerbose)
     {
-        if (optVerbose == 2)
-        {
-            WarningInFunction
-                << "-verbose -verbose switches on full debugging"
-                << nl << endl;
-            ::debug = 1;
-        }
-
         // Report on output
         faMeshDistributor::verbose_ = 1;
         parPointFieldDistributor::verbose_ = 1;
+
+        if (optVerbose > 1)
+        {
+            Info<< "Additional debugging enabled" << nl << endl;
+            ::debug = 1;
+        }
     }
 
     // Disable NaN setting and floating point error trapping. This is to avoid
@@ -1332,19 +1337,33 @@ int main(int argc, char *argv[])
     // fvMeshDistribute.
     Foam::sigFpe::unset(true);
 
-    const wordRes selectedFields;
-    const wordRes selectedLagrangianFields;
 
+    // File handlers (read/write)
+    // ==========================
+
+    // Read handler on processors with a volMesh
+    refPtr<fileOperation> volMeshReadHandler;
+
+    // Read handler on processors with an areaMesh
+    refPtr<fileOperation> areaMeshReadHandler;
+
+    // Handler for master-only operation (read/writing from/to undecomposed)
+    // - only the 'real' master, not io-rank masters
+    refPtr<fileOperation> masterOnlyHandler;
+
+    if (UPstream::master(UPstream::worldComm))
+    {
+        const bool oldParRun = UPstream::parRun(false);
+
+        masterOnlyHandler = fileOperation::NewUncollated();
+
+        UPstream::parRun(oldParRun);
+    }
 
     if (decompose)
     {
-        InfoOrPout<< "Decomposing case (like decomposePar)" << nl << endl;
+        InfoOrPout<< "Decomposing case (like decomposePar)"
-        if (reconstruct)
+            << nl << endl;
-        {
-            FatalErrorInFunction
-                << "Cannot specify both -decompose and -reconstruct"
-                << exit(FatalError);
-        }
     }
     else if (reconstruct)
     {
@@ -1374,19 +1393,13 @@ int main(int argc, char *argv[])
     if ((decompose || reconstruct) && !overwrite)
     {
         overwrite = true;
-        WarningInFunction
-            << "Working in -decompose or -reconstruct mode:"
-               " automatically implies -overwrite" << nl << endl;
-    }
 
-    if (!Pstream::parRun())
+        Warning
-    {
+            << nl << "    "
-        FatalErrorInFunction
+            << "Added implicit -overwrite for decompose/reconstruct modes"
-            << ": This utility can only be run parallel"
+            << nl << endl;
-            << exit(FatalError);
     }
 
-
     if
     (
         fileHandler().ioRanks().contains(UPstream::myProcNo())
@@ -1401,7 +1414,7 @@ int main(int argc, char *argv[])
     }
 
 
-    if (!nfs)
+    if (hasDistributedFiles)
     {
         InfoOrPout<< "Detected multiple roots i.e. non-nfs running"
             << nl << endl;
@@ -1498,7 +1511,7 @@ int main(int argc, char *argv[])
     // This will read the same controlDict but might have a different
     // set of times so enforce same times
 
-    if (!nfs)
+    if (hasDistributedFiles)
     {
         InfoOrPout<< "Creating time directories for undecomposed Time"
             << " on all processors" << nl << endl;

applications/utilities/parallelProcessing/redistributePar/unmappedPassivePositionParticleCloud.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,8 +36,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef unmappedPassivePositionParticleCloud_H
+#ifndef Foam_unmappedPassivePositionParticleCloud_H
-#define unmappedPassivePositionParticleCloud_H
+#define Foam_unmappedPassivePositionParticleCloud_H
 
 #include "passivePositionParticleCloud.H"
 
@@ -54,12 +54,11 @@ class unmappedPassivePositionParticleCloud
 :
     public passivePositionParticleCloud
 {
-
 public:
 
     // Constructors
 
-        //- Construct given mesh
+        //- Read construct
         explicit unmappedPassivePositionParticleCloud
         (
             const polyMesh& mesh,
@@ -70,15 +69,15 @@ public:
             passivePositionParticleCloud(mesh, cloudName, readFields)
         {}
 
-        //- Construct from mesh, cloud name, and a list of particles
+        //- Construct without particles
         unmappedPassivePositionParticleCloud
         (
             const polyMesh& mesh,
-            const word& cloudName,
+            const Foam::zero,
-            const IDLList<passivePositionParticle>& particles
+            const word& cloudName = cloud::defaultName
         )
         :
-            passivePositionParticleCloud(mesh, cloudName, particles)
+            passivePositionParticleCloud(mesh, Foam::zero{}, cloudName)
         {}
 
 
@@ -89,7 +88,7 @@ public:
     // Member Functions
 
         //- Switch off remapping of cells of particles when
-        //  mesh topology changes
+        //- mesh topology changes
         virtual void autoMap(const mapPolyMesh&)
         {}
 

applications/utilities/preProcessing/mapFields/MapLagrangianFields.H

@@ -111,8 +111,6 @@ void MapLagrangianFields
     const labelList& addParticles
 )
 {
-    const fvMesh& meshTarget = meshToMesh0Interp.toMesh();
-
     MapLagrangianFields
     <
         IOField<Type>,

applications/utilities/preProcessing/mapFields/mapLagrangian.C

@@ -133,6 +133,7 @@ void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
                 cloudDir,
                 false
             );
+
             Info<< "    read " << sourceParcels.size()
                 << " parcels from source mesh." << endl;
 
@@ -140,8 +141,8 @@ void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
             passiveParticleCloud targetParcels
             (
                 meshTarget,
-                cloudDir,
+                Foam::zero{},
-                IDLList<passiveParticle>()
+                cloudDir
             );
 
             passiveParticle::trackingData td(targetParcels);

applications/utilities/preProcessing/mapFieldsPar/mapLagrangian.C

@@ -129,6 +129,7 @@ void mapLagrangian(const meshToMesh& interp)
                 cloudDir,
                 false
             );
+
             Info<< "    read " << sourceParcels.size()
                 << " parcels from source mesh." << endl;
 
@@ -136,8 +137,8 @@ void mapLagrangian(const meshToMesh& interp)
             passiveParticleCloud targetParcels
             (
                 meshTarget,
-                cloudDir,
+                Foam::zero{},
-                IDLList<passiveParticle>()
+                cloudDir
             );
 
             passiveParticle::trackingData td(targetParcels);

bin/tools/foamPackRelease

@@ -404,6 +404,68 @@ fi
 
 #------------------------------------------------------------------------------
 # Add in modules
+
+#
+# Recursive addition of submodule content.
+# NB: must be called from within the respective parent directory.
+# Example,
+#
+#    packModule abc  (implied cd)
+#       packModule abc/def
+#          packModule abc/def/hij
+#
+packModule()
+{
+    local parent="$1"
+
+    if [ -n "$parent" ]
+    then
+    (
+        cd "${parent##*/}" 2>/dev/null || exit
+
+        git ls-tree HEAD | \
+        while read mode gittype sha1 module
+        do
+            [ "$gittype" == commit ] || continue
+
+            case "$module" in
+            (. | ./)
+                echo "# module not initialized? : $parent"
+                continue
+                ;;
+            esac
+
+            echo
+            echo "# module"
+            echo "module=\""$parent${parent:+/}$module"\""
+            echo "commit=\""$sha1"\""
+            echo "tmpTarFile=\""$tarName-${module##*/}".tar\""
+            echo
+            echo '('
+            echo 'cd "$module" || exit'
+            echo 'newPrefix="$prefixDir${prefixDir:+/}$module"'
+            echo 'git -c tar.umask=user archive --format=tar --prefix="$newPrefix/" -o "$outputDir/$tmpTarFile" "$commit" || exit'
+            echo '# Without test, validation dirs (potentially large)'
+            echo 'tar --delete -f "$outputDir/$tmpTarFile" "$newPrefix/test" "$newPrefix/validation" 2>/dev/null'
+            echo 'tar -Af "$outputDir/$tarName.tar" "$outputDir/$tmpTarFile"'
+            echo 'rm -f "$outputDir/$tmpTarFile"'
+            echo '{'
+            echo '  echo'
+            echo '  echo "$module"'
+            echo '  echo commit="$commit"'
+            echo '  echo'
+            echo '  # Without test, validation dirs'
+            echo '  git ls-tree -r "$commit" | sed -e '"'"'/\ttest\//d;/\tvalidation\//d'"'"
+            echo '} >> "$outputDir/$manifest1"'
+            echo ')'
+
+            packModule "$parent/$module"
+        done
+    )
+    fi
+}
+
+
 if [ "$withModules" != false ]
 then
     echo
@@ -415,33 +477,7 @@ then
     echo '  echo head="$head"'
     echo '} > "$outputDir/$manifest1"'
 
-    git --git-dir="$gitbase/.git" ls-tree "$head" modules/ | \
+    packModule modules
-        while read mode gittype sha1 module
-    do
-        [ "$gittype" == commit ] || continue
-
-        echo
-        echo "module=\""$module"\""
-        echo "commit=\""$sha1"\""
-        echo "tarModule=\""$tarName-${module##*/}"\""
-        echo
-        echo 'if pushd "$module"; then'
-        echo 'moduleDir="$prefixDir${prefixDir:+/}$module"'
-        echo 'git -c tar.umask=user archive --format=tar --prefix="$moduleDir/" -o "$outputDir/$tarModule.tar" "$commit"'
-        echo '# Without test, validation dirs (potentially large)'
-        echo 'tar --delete -f "$outputDir/$tarModule.tar" "$moduleDir/test" "$moduleDir/validation" 2>/dev/null'
-        echo 'tar -Af "$outputDir/$tarName.tar" "$outputDir/$tarModule.tar"'
-        echo 'rm -f "$outputDir/$tarModule.tar"'
-        echo '{'
-        echo '  echo'
-        echo '  echo "$module"'
-        echo '  echo commit="$commit"'
-        echo '  echo'
-        echo '  # Without test, validation dirs'
-        echo '  git ls-tree -r "$commit" | sed -e '"'"'/\ttest\//d;/\tvalidation\//d'"'"
-        echo '} >> "$outputDir/$manifest1"'
-        echo 'popd; fi'
-    done
 
     echo
     echo '{ echo; echo "# End"; } >> "$outputDir/$manifest1"'
@@ -459,8 +495,9 @@ echo 'rm -f "$buildInfo" "$manifest0" "$manifest1"'
 echo 'popd; fi'
 
 echo
+echo "# -----------------------"
 echo "# End of creating archive"
-echo
+echo "# -----------------------"
 
 #------------------------------------------------------------------------------
 # Compression

doc/Build.md

@@ -1,23 +1,23 @@
 ## Getting the code
 
-Links to all code packs are available on https://dl.openfoam.com. For OpenFOAM-v2306:
+Links to all code packs are available on https://dl.openfoam.com. For OpenFOAM-v2312:
 
 - https://dl.openfoam.com/source/latest/
-- Source: https://dl.openfoam.com/source/v2306/OpenFOAM-v2306.tgz
+- Source: https://dl.openfoam.com/source/v2312/OpenFOAM-v2312.tgz
-- ThirdParty: https://dl.openfoam.com/source/v2306/ThirdParty-v2306.tgz
+- ThirdParty: https://dl.openfoam.com/source/v2312/ThirdParty-v2312.tgz
 
 ## OpenFOAM® Quick Build Guide
 
 Prior to building, ensure that the [system requirements][link openfoam-require]
 are satisfied (including any special [cross-compiling][wiki-cross-compile]
 considerations), and source the correct OpenFOAM environment.
-For example, for the OpenFOAM-v2306 version:
+For example, for the OpenFOAM-v2312 version:
 ```
-source <installation path>/OpenFOAM-v2306/etc/bashrc
+source <installation path>/OpenFOAM-v2312/etc/bashrc
 ```
 e.g. if installed under the `~/openfoam` directory
 ```
-source ~/openfoam/OpenFOAM-v2306/etc/bashrc
+source ~/openfoam/OpenFOAM-v2312/etc/bashrc
 ```
 
 

doc/BuildIssues.md

@@ -1,4 +1,4 @@
-## Known Build Issues (v2306, v2212, v2206, v2112, v2106, v2012)
+## Known Build Issues (v2312, v2306, v2212, v2206, v2112, v2106, v2012)
 
 ### Windows cross-compilation
 

etc/config.csh/adios2

@@ -5,7 +5,7 @@
 #   \\  /    A nd           | www.openfoam.com
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
-#     Copyright (C) 2017-2022 OpenCFD Ltd.
+#     Copyright (C) 2017-2023 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -20,7 +20,7 @@
 #------------------------------------------------------------------------------
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
 
-set adios2_version=ADIOS2-2.7.1
+set adios2_version=ADIOS2-2.9.2
 setenv ADIOS2_ARCH_PATH "$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$adios2_version"
 
 # END OF (NORMAL) USER EDITABLE PART

etc/config.sh/adios2

@@ -5,7 +5,7 @@
 #   \\  /    A nd           | www.openfoam.com
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
-#     Copyright (C) 2017-2022 OpenCFD Ltd.
+#     Copyright (C) 2017-2023 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -21,7 +21,7 @@
 #------------------------------------------------------------------------------
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
 
-adios2_version=ADIOS2-2.8.3
+adios2_version=ADIOS2-2.9.2
 export ADIOS2_ARCH_PATH="$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER/$adios2_version"
 
 # END OF (NORMAL) USER EDITABLE PART

src/OSspecific/MSwindows/Allwmake

@@ -19,6 +19,7 @@ fi
 unset COMP_FLAGS LINK_FLAGS
 
 # Make object (non-shared by default)
-wmake $targetType
+# Never want/need openmp, especially for static objects
+wmake -no-openmp $targetType
 
 #------------------------------------------------------------------------------

src/OSspecific/POSIX/Allwmake

@@ -14,6 +14,7 @@ then
 fi
 
 # Make object (non-shared by default)
-wmake $targetType
+# Never want/need openmp, especially for static objects
+wmake -no-openmp $targetType
 
 #------------------------------------------------------------------------------

src/OpenFOAM/db/IOobject/IOobject.H

@@ -656,21 +656,6 @@ public:
         template<class Type>
         void warnNoRereading() const;
 
-        //- Read header (localFilePath only) with optional searching
-        template<>
-        bool typeHeaderOk<void>
-        (
-            //! ignored for \c void
-            const bool checkType,
-            const bool search,
-            //! ignored for \c void
-            const bool verbose
-        );
-
-        //- Call localFilePath for \c void type
-        template<>
-        fileName typeFilePath<void>(const bool search) const;
-
 
     // Writing
 

src/OpenFOAM/db/IOobject/IOobjectReadHeader.C

@@ -272,37 +272,4 @@ bool Foam::IOobject::readAndCheckHeader
 }
 
 
-// * * * * * * * * * * * * Template Specializations  * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-template<>
-bool IOobject::typeHeaderOk<void>
-(
-    const bool checkType,
-    const bool search,
-    const bool verbose
-)
-{
-    return readAndCheckHeader
-    (
-        false,      // global = false
-        word::null,
-        false,      // checkType = false (not meaningful)
-        search,
-        false       // verbose = false (not meaningful)
-    );
-}
-
-
-template<>
-fileName IOobject::typeFilePath<void>(const bool search) const
-{
-    return this->localFilePath(word::null, search);
-}
-
-} // End namespace Foam
-
-
 // ************************************************************************* //

src/OpenFOAM/db/objectRegistry/objectRegistry.H

@@ -633,7 +633,13 @@ public:
         //FUTURE void addTemporaryObject(const word& name) const;
 
         //- True if given name is in the cacheTemporaryObjects set
-        bool cacheTemporaryObject(const word& name) const;
+        bool is_cacheTemporaryObject(const word& name) const;
+
+        //- True if name of object is in the cacheTemporaryObjects set
+        bool is_cacheTemporaryObject(const regIOobject* io) const;
+
+        //- True if name of object is in the cacheTemporaryObjects set
+        bool is_cacheTemporaryObject(const regIOobject& io) const;
 
         //- Cache the given object. Moves content and stores
         template<class Type>
@@ -643,7 +649,7 @@ public:
         void resetCacheTemporaryObject(const regIOobject* io) const;
 
         //- Reset the cache state of the given object
-        //  in the cacheTemporaryObjects set
+        //- in the cacheTemporaryObjects set
         void resetCacheTemporaryObject(const regIOobject& io) const;
 
         //- Check that all objects specified in the cacheTemporaryObjects

src/OpenFOAM/db/objectRegistry/objectRegistryCache.C

@@ -65,6 +65,12 @@ void Foam::objectRegistry::readCacheTemporaryObjects() const
         {
             cacheTemporaryObjects_.emplace(objName, false, false);
         }
+
+        if (objectRegistry::debug)
+        {
+            Info<< "objectRegistry::cacheTemporaryObjects : "
+                << flatOutput(objectNames) << endl;
+        }
     }
 }
 
@@ -75,6 +81,8 @@ void Foam::objectRegistry::deleteCachedObject(regIOobject* io) const
     {
         io->release();     // Relinquish any ownership by registry
         io->checkOut();
+        // Additional safety - not certain this is actually needed...
+        io->rename(io->name() + "-Cache");
         delete io;
     }
 }
@@ -90,12 +98,24 @@ void Foam::objectRegistry::deleteCachedObject(regIOobject* io) const
 // }
 
 
-bool Foam::objectRegistry::cacheTemporaryObject
+bool Foam::objectRegistry::is_cacheTemporaryObject
 (
     const word& name
 ) const
 {
-    return cacheTemporaryObjects_.found(name);
+    return cacheTemporaryObjects_.contains(name);
+}
+
+
+bool Foam::objectRegistry::is_cacheTemporaryObject(const regIOobject* io) const
+{
+    return io && cacheTemporaryObjects_.contains(io->name());
+}
+
+
+bool Foam::objectRegistry::is_cacheTemporaryObject(const regIOobject& io) const
+{
+    return cacheTemporaryObjects_.contains(io.name());
 }
 
 

src/OpenFOAM/db/regIOobject/regIOobject.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,7 +56,7 @@ Foam::regIOobject::regIOobject(const IOobject& io, const bool isTimeObject)
     metaDataPtr_(nullptr),
     isPtr_(nullptr)
 {
-    if (registerObject())
+    if (IOobject::registerObject())
     {
         // Register (check-in) with objectRegistry if requested
         checkIn();
@@ -137,7 +137,7 @@ Foam::regIOobject::regIOobject
     metaDataPtr_(rio.metaDataPtr_.clone()),
     isPtr_(nullptr)
 {
-    if (registerObject())
+    if (IOobject::registerObject())
     {
         checkIn();
     }
@@ -194,7 +194,7 @@ bool Foam::regIOobject::checkIn()
     {
         // multiple checkin of same object is disallowed - this would mess up
         // any mapping
-        registered_ = db().checkIn(*this);
+        registered_ = db().checkIn(this);
 
         // check-in on defaultRegion is allowed to fail, since subsetted meshes
         // are created with the same name as their originating mesh
@@ -205,17 +205,18 @@ bool Foam::regIOobject::checkIn()
                 // for ease of finding where attempted duplicate check-in
                 // originated
                 FatalErrorInFunction
-                    << "failed to register object " << objectPath()
+                    << "Failed to register: " << name() << ' '
-                    << " the name already exists in the objectRegistry" << endl
+                    << objectRelPath()
-                    << "Contents:" << db().sortedToc()
+                    << " : the name already exists in the registry" << nl
+                    << "Contents:" << db().sortedToc() << endl
                     << abort(FatalError);
             }
             else
             {
                 WarningInFunction
-                    << "failed to register object " << objectPath()
+                    << "Failed to register: " << name() << ' '
-                    << " the name already exists in the objectRegistry"
+                    << objectRelPath()
-                    << endl;
+                    << " : the name already exists in the registry" << endl;
             }
         }
     }
@@ -276,6 +277,23 @@ Foam::label Foam::regIOobject::addWatch(const fileName& f)
 
 void Foam::regIOobject::addWatch()
 {
+    // Everyone or just master
+    const bool masterOnly
+    (
+        global()
+     && (
+            IOobject::fileModificationChecking == IOobject::timeStampMaster
+         || IOobject::fileModificationChecking == IOobject::inotifyMaster
+        )
+    );
+
+    // if (debug)
+    // {
+    //     Pout<< "regIOobject::addWatch " << watchIndices_.size()
+    //         << " indices master-only:" << masterOnly
+    //         << " watchFiles: " << watchFiles_ << endl;
+    // }
+
     if
     (
         registered_
@@ -284,7 +302,7 @@ void Foam::regIOobject::addWatch()
     )
     {
         fileName f = filePath();
-        if (!f.size())
+        if (f.empty())
         {
             // We don't have this file but would like to re-read it.
             // Possibly if master-only reading mode.
@@ -302,14 +320,6 @@ void Foam::regIOobject::addWatch()
 
         // If master-only reading only the master will have all dependencies
         // so broadcast these to other ranks
-        const bool masterOnly
-        (
-            global()
-         && (
-                IOobject::fileModificationChecking == IOobject::timeStampMaster
-             || IOobject::fileModificationChecking == IOobject::inotifyMaster
-            )
-        );
 
         if (masterOnly && UPstream::parRun())
         {
@@ -327,8 +337,13 @@ void Foam::regIOobject::addWatch()
 
                 UPstream::parRun(oldParRun);
             }
-            Pstream::broadcast(filesToWatch);
 
+            Pstream::broadcasts
+            (
+                UPstream::worldComm,
+                filesToWatch,
+                watchFiles_
+            );
 
             // Add master files in same order
             if (!UPstream::master())
@@ -474,7 +489,7 @@ void Foam::regIOobject::rename(const word& newName)
 
     IOobject::rename(newName);
 
-    if (registerObject())
+    if (IOobject::registerObject())
     {
         // Re-register object with objectRegistry
         checkIn();
@@ -519,7 +534,7 @@ void Foam::regIOobject::operator=(const IOobject& io)
 
     IOobject::operator=(io);
 
-    if (registerObject())
+    if (IOobject::registerObject())
     {
         // Re-register object with objectRegistry
         checkIn();

src/OpenFOAM/db/regIOobject/regIOobject.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -116,9 +116,6 @@ private:
         //- Construct object stream, read header if not already constructed
         void readStream(const bool readOnProc);
 
-        //- No copy assignment
-        void operator=(const regIOobject&) = delete;
-
 
 public:
 
@@ -127,6 +124,12 @@ public:
         friend class fileOperations::uncollatedFileOperation;
 
 
+    // Generated Methods
+
+        //- No copy assignment
+        void operator=(const regIOobject&) = delete;
+
+
     // Static Data
 
         //- Runtime type information
@@ -366,6 +369,7 @@ public:
     // Member Operators
 
         //- Copy assignment
+        FOAM_DEPRECATED_STRICT(2023-12, "possibly remove in the future")
         void operator=(const IOobject& io);
 
 

src/OpenFOAM/db/regIOobject/regIOobjectI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -145,6 +145,8 @@ inline Type& Foam::regIOobject::store(tmp<Type>& ptr)
     if (ptr.is_pointer())
     {
         // Acquire ownership, pass management to objectRegistry
+
+        ptr.protect(false);  // Storing (ie, not cached/protected)
         p = ptr.ptr();
 
         store(p);

src/OpenFOAM/db/regIOobject/regIOobjectRead.C

@@ -100,7 +100,7 @@ void Foam::regIOobject::readStream(const bool readOnProc)
         if (watchIndices_.size())
         {
             // File is being watched. Read exact file that is being watched.
-            objPath = fileHandler().getFile(watchIndices_.last());
+            objPath = fileHandler().getFile(watchIndices_.back());
         }
         else
         {
@@ -190,26 +190,7 @@ bool Foam::regIOobject::read()
     // Note: cannot do anything in readStream itself since this is used by
     // e.g. GeometricField.
 
-
+    // Everyone or just master
-    // Save old watchIndices and clear (so the list of included files can
-    // change)
-    fileNameList oldWatchFiles;
-    if (watchIndices_.size())
-    {
-        oldWatchFiles.setSize(watchIndices_.size());
-        forAll(watchIndices_, i)
-        {
-            oldWatchFiles[i] = fileHandler().getFile(watchIndices_[i]);
-        }
-        forAllReverse(watchIndices_, i)
-        {
-            fileHandler().removeWatch(watchIndices_[i]);
-        }
-        watchIndices_.clear();
-    }
-
-
-    // Read
     const bool masterOnly
     (
         global()
@@ -219,12 +200,26 @@ bool Foam::regIOobject::read()
         )
     );
 
+    // Remove old watches (indices) and clear:
+    // so the list of included files can change
+
+    const bool needWatch(!watchIndices_.empty());
+
+    if (!watchIndices_.empty())
+    {
+        forAllReverse(watchIndices_, i)
+        {
+            fileHandler().removeWatch(watchIndices_[i]);
+        }
+        watchIndices_.clear();
+    }
+
     // Note: IOstream::binary flag is for all the processor comms. (Only for
     //       dictionaries should it be ascii)
     bool ok =
         fileHandler().read(*this, masterOnly, IOstreamOption::BINARY, type());
 
-    if (oldWatchFiles.size())
+    if (needWatch)
     {
         // Re-watch master file
         addWatch();
@@ -264,9 +259,9 @@ bool Foam::regIOobject::readIfModified()
 
     if (modified != -1)
     {
-        const fileName fName = fileHandler().getFile(watchIndices_.last());
+        const fileName fName = fileHandler().getFile(watchIndices_.back());
 
-        if (modified == watchIndices_.last())
+        if (modified == watchIndices_.back())
         {
             InfoInFunction
                 << "    Re-reading object " << name()

src/OpenFOAM/db/regIOobject/regIOobjectWrite.C

@@ -94,7 +94,7 @@ bool Foam::regIOobject::writeObject
     }
 
 
-    // Everyone check or just master
+    // Everyone or just master
     const bool masterOnly
     (
         isGlobal
@@ -122,9 +122,9 @@ bool Foam::regIOobject::writeObject
 
     // Only update the lastModified_ time if this object is re-readable,
     // i.e. lastModified_ is already set
-    if (watchIndices_.size())
+    if (!watchIndices_.empty())
     {
-        fileHandler().setUnmodified(watchIndices_.last());
+        fileHandler().setUnmodified(watchIndices_.back());
     }
 
     return osGood;

src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C

@@ -376,6 +376,16 @@ Foam::DimensionedField<Type, GeoMesh>::clone() const
 }
 
 
+// * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * * //
+
+template<class Type, class GeoMesh>
+Foam::DimensionedField<Type, GeoMesh>::~DimensionedField()
+{
+    // FUTURE: register cache field info
+    // // this->db().cacheTemporaryObject(*this);
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type, class GeoMesh>

src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H

@@ -117,6 +117,17 @@ private:
 
         void readIfPresent(const word& fieldDictEntry = "value");
 
+        //- Implementation for 'New' with specified registerObject preference.
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        template<class... Args>
+        static tmp<DimensionedField<Type, GeoMesh>> New_impl
+        (
+            IOobjectOption::registerOption regOpt,
+            const word& name,
+            const Mesh& mesh,
+            Args&&... args
+        );
 
 public:
 
@@ -298,10 +309,25 @@ public:
 
     // Factory Methods
 
-        //- Return tmp field from name, mesh, dimensions,
+        //- Return tmp field (NO_READ, NO_WRITE)
-        //- copy of internal field.
+        //- from name, mesh, dimensions, copy of internal field.
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- [Takes current timeName from the mesh registry].
-        //  current timeName from the mesh registry
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        static tmp<DimensionedField<Type, GeoMesh>> New
+        (
+            const word& name,
+            IOobjectOption::registerOption regOpt,
+            const Mesh& mesh,
+            const dimensionSet& dims,
+            const Field<Type>& iField
+        );
+
+        //- Return tmp field (NO_READ, NO_WRITE)
+        //- from name, mesh, dimensions, copy of internal field.
+        //- [Takes current timeName from the mesh registry].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<DimensionedField<Type, GeoMesh>> New
         (
             const word& name,
@@ -310,10 +336,25 @@ public:
             const Field<Type>& iField
         );
 
-        //- Return tmp field from name, mesh, dimensions,
+        //- Return tmp field (NO_READ, NO_WRITE)
-        //- moved internal field contents.
+        //- from name, mesh, dimensions, move internal field contents.
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- [Takes current timeName from the mesh registry].
-        //  current timeName from the mesh registry
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        static tmp<DimensionedField<Type, GeoMesh>> New
+        (
+            const word& name,
+            IOobjectOption::registerOption regOpt,
+            const Mesh& mesh,
+            const dimensionSet& dims,
+            Field<Type>&& iField
+        );
+
+        //- Return tmp field (NO_READ, NO_WRITE)
+        //- from name, mesh, dimensions, move internal field contents.
+        //- [Takes current timeName from the mesh registry].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<DimensionedField<Type, GeoMesh>> New
         (
             const word& name,
@@ -322,9 +363,24 @@ public:
             Field<Type>&& iField
         );
 
-        //- Return tmp field from name, mesh, dimensions.
+        //- Return tmp field (NO_READ, NO_WRITE)
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- from name, mesh, dimensions.
-        //  current timeName from the mesh registry
+        //- [Takes current timeName from the mesh registry].
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        static tmp<DimensionedField<Type, GeoMesh>> New
+        (
+            const word& name,
+            IOobjectOption::registerOption regOpt,
+            const Mesh& mesh,
+            const dimensionSet& dims
+        );
+
+        //- Return tmp field (NO_READ, NO_WRITE)
+        //- from name, mesh, dimensions.
+        //- [Takes current timeName from the mesh registry].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<DimensionedField<Type, GeoMesh>> New
         (
             const word& name,
@@ -332,9 +388,25 @@ public:
             const dimensionSet& dims
         );
 
-        //- Return uniform value tmp field from name, mesh, dimensions, value.
+        //- Return tmp field (NO_READ, NO_WRITE)
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- from name, mesh, value, dimensions.
-        //  current timeName from the mesh registry
+        //- [Takes current timeName from the mesh registry].
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        static tmp<DimensionedField<Type, GeoMesh>> New
+        (
+            const word& name,
+            IOobjectOption::registerOption regOpt,
+            const Mesh& mesh,
+            const Type& value,
+            const dimensionSet& dims
+        );
+
+        //- Return tmp field (NO_READ, NO_WRITE)
+        //- from name, mesh, value, dimensions.
+        //- [Takes current timeName from the mesh registry].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<DimensionedField<Type, GeoMesh>> New
         (
             const word& name,
@@ -343,9 +415,24 @@ public:
             const dimensionSet& dims
         );
 
-        //- Return tmp field from name, mesh, dimensioned\<Type\>.
+        //- Return tmp field (NO_READ, NO_WRITE)
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- from name, mesh, dimensioned-type.
-        //  current timeName from the mesh registry
+        //- [Takes current timeName from the mesh registry].
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        static tmp<DimensionedField<Type, GeoMesh>> New
+        (
+            const word& name,
+            IOobjectOption::registerOption regOpt,
+            const Mesh& mesh,
+            const dimensioned<Type>& dt
+        );
+
+        //- Return tmp field (NO_READ, NO_WRITE)
+        //- from name, mesh, dimensioned-type.
+        //- [Takes current timeName from the mesh registry].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<DimensionedField<Type, GeoMesh>> New
         (
             const word& name,
@@ -353,19 +440,21 @@ public:
             const dimensioned<Type>& dt
         );
 
-        //- Return renamed tmp field
+        //- Return renamed tmp field (NO_READ, NO_WRITE).
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- [Takes current timeName from the mesh registry].
-        //  current timeName from the mesh registry
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<DimensionedField<Type, GeoMesh>> New
         (
             const word& newName,
             const tmp<DimensionedField<Type, GeoMesh>>& tfld
         );
 
-        //- Construct tmp field based on mesh/registry information from
+        //- Construct tmp field (NO_READ, NO_WRITE)
-        //- an existing field.
+        //- based on mesh/registry information from an existing field.
-        //  Created NO_READ, NO_WRITE, NO_REGISTER, using the instance
+        //- [Takes instance from the field].
-        //  from the field
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         template<class AnyType>
         static tmp<DimensionedField<Type, GeoMesh>> New
         (
@@ -374,10 +463,11 @@ public:
             const dimensionSet& dims
         );
 
-        //- Construct tmp field based on mesh/registry information from
+        //- Construct tmp field (NO_READ, NO_WRITE)
-        //- an existing field and initialise with value.
+        //- based on mesh/registry information from an existing field
-        //  Created NO_READ, NO_WRITE, NO_REGISTER, using the instance
+        //- and initialise with value.
-        //  from the field
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         template<class AnyType>
         static tmp<DimensionedField<Type, GeoMesh>> New
         (
@@ -388,7 +478,7 @@ public:
 
 
     //- Destructor
-    virtual ~DimensionedField() = default;
+    virtual ~DimensionedField();
 
 
     // Member Functions
@@ -469,11 +559,15 @@ public:
         ) const;
 
 
-        // Write
+    // Write
 
-            bool writeData(Ostream& os, const word& fieldDictEntry) const;
+        //- Write dimensions, oriented flag (if valid) and the
+        //- field data as a dictionary entry with the specified name.
+        bool writeData(Ostream& os, const word& fieldDictEntry) const;
 
-            bool writeData(Ostream& os) const;
+        //- The writeData function (required by regIOobject),
+        //- call writeData with dictionary entry name = "value"
+        bool writeData(Ostream& os) const;
 
 
     // Member Operators

src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldNew.C

@@ -25,23 +25,21 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type, class GeoMesh>
+template<class... Args>
 Foam::tmp<Foam::DimensionedField<Type, GeoMesh>>
-Foam::DimensionedField<Type, GeoMesh>::New
+Foam::DimensionedField<Type, GeoMesh>::New_impl
 (
+    IOobjectOption::registerOption regOpt,
     const word& name,
     const Mesh& mesh,
-    const dimensionSet& dims,
+    Args&&... args
-    const Field<Type>& iField
 )
 {
-    const bool caching = mesh.thisDb().cacheTemporaryObject(name);
+    auto ptr = tmp<DimensionedField<Type, GeoMesh>>::New
-
-    return tmp<DimensionedField<Type, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
         IOobject
         (
             name,
@@ -49,9 +47,48 @@ Foam::DimensionedField<Type, GeoMesh>::New
             mesh.thisDb(),
             IOobjectOption::NO_READ,
             IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
+            IOobjectOption::NO_REGISTER
         ),
         mesh,
+        std::forward<Args>(args)...
+    );
+
+    if (IOobjectOption::REGISTER == regOpt)
+    {
+        ptr->checkIn();
+    }
+    else if
+    (
+        // LEGACY_REGISTER: detect if caching is desired
+        (IOobjectOption::LEGACY_REGISTER == regOpt)
+     && ptr->db().is_cacheTemporaryObject(ptr.get())
+    )
+    {
+        ptr.protect(true);
+        ptr->checkIn();
+    }
+    return ptr;
+}
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+template<class Type, class GeoMesh>
+Foam::tmp<Foam::DimensionedField<Type, GeoMesh>>
+Foam::DimensionedField<Type, GeoMesh>::New
+(
+    const word& name,
+    IOobjectOption::registerOption regOpt,
+    const Mesh& mesh,
+    const dimensionSet& dims,
+    const Field<Type>& iField
+)
+{
+    return DimensionedField<Type, GeoMesh>::New_impl
+    (
+        regOpt,
+        name,
+        mesh,
         dims,
         iField
     );
@@ -65,23 +102,35 @@ Foam::DimensionedField<Type, GeoMesh>::New
     const word& name,
     const Mesh& mesh,
     const dimensionSet& dims,
+    const Field<Type>& iField
+)
+{
+    return DimensionedField<Type, GeoMesh>::New_impl
+    (
+        IOobjectOption::LEGACY_REGISTER,
+        name,
+        mesh,
+        dims,
+        iField
+    );
+}
+
+
+template<class Type, class GeoMesh>
+Foam::tmp<Foam::DimensionedField<Type, GeoMesh>>
+Foam::DimensionedField<Type, GeoMesh>::New
+(
+    const word& name,
+    IOobjectOption::registerOption regOpt,
+    const Mesh& mesh,
+    const dimensionSet& dims,
     Field<Type>&& iField
 )
 {
-    const bool caching = mesh.thisDb().cacheTemporaryObject(name);
+    return DimensionedField<Type, GeoMesh>::New_impl
-
-    return tmp<DimensionedField<Type, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
+        regOpt,
-        IOobject
+        name,
-        (
-            name,
-            mesh.thisDb().time().timeName(),
-            mesh.thisDb(),
-            IOobjectOption::NO_READ,
-            IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
-        ),
         mesh,
         dims,
         std::move(iField)
@@ -95,26 +144,81 @@ Foam::DimensionedField<Type, GeoMesh>::New
 (
     const word& name,
     const Mesh& mesh,
+    const dimensionSet& dims,
+    Field<Type>&& iField
+)
+{
+    return DimensionedField<Type, GeoMesh>::New_impl
+    (
+        IOobjectOption::LEGACY_REGISTER,
+        name,
+        mesh,
+        dims,
+        std::move(iField)
+    );
+}
+
+
+template<class Type, class GeoMesh>
+Foam::tmp<Foam::DimensionedField<Type, GeoMesh>>
+Foam::DimensionedField<Type, GeoMesh>::New
+(
+    const word& name,
+    IOobjectOption::registerOption regOpt,
+    const Mesh& mesh,
     const dimensionSet& dims
 )
 {
-    const bool caching = mesh.thisDb().cacheTemporaryObject(name);
+    return DimensionedField<Type, GeoMesh>::New_impl
-
-    return tmp<DimensionedField<Type, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
+        regOpt,
-        IOobject
+        name,
-        (
-            name,
-            mesh.thisDb().time().timeName(),
-            mesh.thisDb(),
-            IOobjectOption::NO_READ,
-            IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
-        ),
         mesh,
         dims,
-        false  // checkIOFlags off
+        false  // No checkIOFlags (is NO_READ anyhow)
+    );
+}
+
+
+template<class Type, class GeoMesh>
+Foam::tmp<Foam::DimensionedField<Type, GeoMesh>>
+Foam::DimensionedField<Type, GeoMesh>::New
+(
+    const word& name,
+    const Mesh& mesh,
+    const dimensionSet& dims
+)
+{
+    return DimensionedField<Type, GeoMesh>::New_impl
+    (
+        IOobjectOption::LEGACY_REGISTER,
+        name,
+        mesh,
+        dims,
+        false  // No checkIOFlags (is NO_READ anyhow)
+    );
+}
+
+
+template<class Type, class GeoMesh>
+Foam::tmp<Foam::DimensionedField<Type, GeoMesh>>
+Foam::DimensionedField<Type, GeoMesh>::New
+(
+    const word& name,
+    IOobjectOption::registerOption regOpt,
+    const Mesh& mesh,
+    const Type& value,
+    const dimensionSet& dims
+)
+{
+    return DimensionedField<Type, GeoMesh>::New_impl
+    (
+        regOpt,
+        name,
+        mesh,
+        value,
+        dims,
+        false  // No checkIOFlags (is NO_READ anyhow)
     );
 }
 
@@ -129,24 +233,35 @@ Foam::DimensionedField<Type, GeoMesh>::New
     const dimensionSet& dims
 )
 {
-    const bool caching = mesh.thisDb().cacheTemporaryObject(name);
+    return DimensionedField<Type, GeoMesh>::New_impl
-
-    return tmp<DimensionedField<Type, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
+        IOobjectOption::LEGACY_REGISTER,
-        IOobject
+        name,
-        (
-            name,
-            mesh.thisDb().time().timeName(),
-            mesh.thisDb(),
-            IOobjectOption::NO_READ,
-            IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
-        ),
         mesh,
         value,
         dims,
-        false  // checkIOFlags off
+        false  // No checkIOFlags (is NO_READ anyhow)
+    );
+}
+
+
+template<class Type, class GeoMesh>
+Foam::tmp<Foam::DimensionedField<Type, GeoMesh>>
+Foam::DimensionedField<Type, GeoMesh>::New
+(
+    const word& name,
+    IOobjectOption::registerOption regOpt,
+    const Mesh& mesh,
+    const dimensioned<Type>& dt
+)
+{
+    return DimensionedField<Type, GeoMesh>::New
+    (
+        name,
+        regOpt,
+        mesh,
+        dt.value(),
+        dt.dimensions()
     );
 }
 
@@ -178,11 +293,8 @@ Foam::DimensionedField<Type, GeoMesh>::New
     const tmp<DimensionedField<Type, GeoMesh>>& tfld
 )
 {
-    const bool caching = tfld().db().cacheTemporaryObject(name);
+    auto ptr = tmp<DimensionedField<Type, GeoMesh>>::New
-
-    return tmp<DimensionedField<Type, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
         IOobject
         (
             name,
@@ -191,10 +303,20 @@ Foam::DimensionedField<Type, GeoMesh>::New
             tfld().db(),
             IOobjectOption::NO_READ,
             IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
+            IOobjectOption::NO_REGISTER
         ),
         tfld
     );
+
+    if
+    (
+        ptr->db().is_cacheTemporaryObject(ptr.get())
+    )
+    {
+        ptr.protect(true);
+        ptr->checkIn();
+    }
+    return ptr;
 }
 
 
@@ -208,11 +330,8 @@ Foam::DimensionedField<Type, GeoMesh>::New
     const dimensionSet& dims
 )
 {
-    const bool caching = fld.db().cacheTemporaryObject(name);
+    auto ptr = tmp<DimensionedField<Type, GeoMesh>>::New
-
-    return tmp<DimensionedField<Type, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
         IOobject
         (
             name,
@@ -220,11 +339,21 @@ Foam::DimensionedField<Type, GeoMesh>::New
             fld.db(),
             IOobjectOption::NO_READ,
             IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
+            IOobjectOption::NO_REGISTER
         ),
         fld.mesh(),
         dims
     );
+
+    if
+    (
+        ptr->db().is_cacheTemporaryObject(ptr.get())
+    )
+    {
+        ptr.protect(true);
+        ptr->checkIn();
+    }
+    return ptr;
 }
 
 
@@ -238,11 +367,8 @@ Foam::DimensionedField<Type, GeoMesh>::New
     const dimensioned<Type>& dt
 )
 {
-    const bool caching = fld.db().cacheTemporaryObject(name);
+    auto ptr = tmp<DimensionedField<Type, GeoMesh>>::New
-
-    return tmp<DimensionedField<Type, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
         IOobject
         (
             name,
@@ -250,12 +376,22 @@ Foam::DimensionedField<Type, GeoMesh>::New
             fld.db(),
             IOobjectOption::NO_READ,
             IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
+            IOobjectOption::NO_REGISTER
         ),
         fld.mesh(),
         dt.value(),
         dt.dimensions()
     );
+
+    if
+    (
+        ptr->db().is_cacheTemporaryObject(ptr.get())
+    )
+    {
+        ptr.protect(true);
+        ptr->checkIn();
+    }
+    return ptr;
 }
 
 

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

@@ -899,6 +899,11 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::~GeometricField()
 
     deleteDemandDrivenData(field0Ptr_);
     deleteDemandDrivenData(fieldPrevIterPtr_);
+
+    // FUTURE: register cache field info
+    // // this->db().cacheTemporaryObject(*this);
+
+    clearOldTimes();
 }
 
 
@@ -955,6 +960,14 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::boundaryFieldRef
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+Foam::label
+Foam::GeometricField<Type, PatchField, GeoMesh>::nOldTimes() const noexcept
+{
+    return (field0Ptr_ ? (field0Ptr_->nOldTimes() + 1) : 0);
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::storeOldTimes() const
 {
@@ -995,18 +1008,6 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::storeOldTime() const
 }
 
 
-template<class Type, template<class> class PatchField, class GeoMesh>
-Foam::label Foam::GeometricField<Type, PatchField, GeoMesh>::nOldTimes() const
-{
-    if (field0Ptr_)
-    {
-        return field0Ptr_->nOldTimes() + 1;
-    }
-
-    return 0;
-}
-
-
 template<class Type, template<class> class PatchField, class GeoMesh>
 const Foam::GeometricField<Type, PatchField, GeoMesh>&
 Foam::GeometricField<Type, PatchField, GeoMesh>::oldTime() const
@@ -1097,6 +1098,14 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::prevIter() const
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+void Foam::GeometricField<Type, PatchField, GeoMesh>::clearOldTimes()
+{
+    deleteDemandDrivenData(field0Ptr_);
+    deleteDemandDrivenData(fieldPrevIterPtr_);
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::
 correctBoundaryConditions()

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H

@@ -134,6 +134,18 @@ private:
         //- Read the field - create the field dictionary on-the-fly
         void readFields();
 
+        //- Implementation for 'New' with specified registerObject preference.
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        template<class... Args>
+        static tmp<GeometricField<Type, PatchField, GeoMesh>> New_impl
+        (
+            IOobjectOption::registerOption regOpt,
+            const word& name,
+            const Mesh& mesh,
+            Args&&... args
+        );
+
 public:
 
     //- Runtime type information
@@ -405,9 +417,25 @@ public:
 
     // Factory Methods
 
-        //- Return tmp field from name, mesh, dimensions and patch type.
+        //- Return tmp field (NO_READ, NO_WRITE)
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- from name, mesh, dimensions and patch type.
-        //  current timeName from the mesh registry
+        //- [Takes current timeName from the mesh registry].
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        static tmp<GeometricField<Type, PatchField, GeoMesh>> New
+        (
+            const word& name,
+            IOobjectOption::registerOption regOpt,
+            const Mesh& mesh,
+            const dimensionSet& dims,
+            const word& patchFieldType = PatchField<Type>::calculatedType()
+        );
+
+        //- Return tmp field (NO_READ, NO_WRITE)
+        //- from name, mesh, dimensions and patch type.
+        //- [Takes current timeName from the mesh registry].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<GeometricField<Type, PatchField, GeoMesh>> New
         (
             const word& name,
@@ -416,10 +444,26 @@ public:
             const word& patchFieldType = PatchField<Type>::calculatedType()
         );
 
-        //- Return tmp field from name, mesh, dimensions,
+        //- Return tmp field (NO_READ, NO_WRITE)
-        //- copy of internal field, with specified patch type.
+        //- from name, mesh, dimensions, copy of internal field and patch type.
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- [Takes current timeName from the mesh registry].
-        //  current timeName from the mesh registry
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        static tmp<GeometricField<Type, PatchField, GeoMesh>> New
+        (
+            const word& name,
+            IOobjectOption::registerOption regOpt,
+            const Mesh& mesh,
+            const dimensionSet& dims,
+            const Field<Type>& iField,
+            const word& patchFieldType = PatchField<Type>::calculatedType()
+        );
+
+        //- Return tmp field (NO_READ, NO_WRITE)
+        //- from name, mesh, dimensions, copy of internal field and patch type.
+        //- [Takes current timeName from the mesh registry].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<GeometricField<Type, PatchField, GeoMesh>> New
         (
             const word& name,
@@ -429,10 +473,28 @@ public:
             const word& patchFieldType = PatchField<Type>::calculatedType()
         );
 
-        //- Return tmp field from name, mesh, dimensions,
+        //- Return tmp field (NO_READ, NO_WRITE)
-        //- moved internal field contents, with specified patch type.
+        //- from name, mesh, dimensions, moved internal field contents
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- and patch type.
-        //  current timeName from the mesh registry
+        //- [Takes current timeName from the mesh registry].
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        static tmp<GeometricField<Type, PatchField, GeoMesh>> New
+        (
+            const word& name,
+            IOobjectOption::registerOption regOpt,
+            const Mesh& mesh,
+            const dimensionSet& dims,
+            Field<Type>&& iField,
+            const word& patchFieldType = PatchField<Type>::calculatedType()
+        );
+
+        //- Return tmp field (NO_READ, NO_WRITE)
+        //- from name, mesh, dimensions, moved internal field contents
+        //- and patch type.
+        //- [Takes current timeName from the mesh registry].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<GeometricField<Type, PatchField, GeoMesh>> New
         (
             const word& name,
@@ -442,10 +504,26 @@ public:
             const word& patchFieldType = PatchField<Type>::calculatedType()
         );
 
-        //- Return tmp field from name, mesh, dimensions, initial field value
+        //- Return tmp field (NO_READ, NO_WRITE)
-        //- and patch type.
+        //- from name, mesh, field value, dimensions and patch type.
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- [Takes current timeName from the mesh registry].
-        //  current timeName from the mesh registry
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        static tmp<GeometricField<Type, PatchField, GeoMesh>> New
+        (
+            const word& name,
+            IOobjectOption::registerOption regOpt,
+            const Mesh& mesh,
+            const Type& value,
+            const dimensionSet& dims,
+            const word& patchFieldType = PatchField<Type>::calculatedType()
+        );
+
+        //- Return tmp field (NO_READ, NO_WRITE)
+        //- from name, mesh, field value, dimensions and patch type.
+        //- [Takes current timeName from the mesh registry].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<GeometricField<Type, PatchField, GeoMesh>> New
         (
             const word& name,
@@ -455,10 +533,27 @@ public:
             const word& patchFieldType = PatchField<Type>::calculatedType()
         );
 
-        //- Return tmp field from name, mesh, dimensioned\<Type\>
+        //- Return tmp field (NO_READ, NO_WRITE)
-        //- and patch types.
+        //- from name, mesh, field value, dimensions and patch field types.
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- [Takes current timeName from the mesh registry].
-        //  current timeName from the mesh registry
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        static tmp<GeometricField<Type, PatchField, GeoMesh>> New
+        (
+            const word& name,
+            IOobjectOption::registerOption regOpt,
+            const Mesh& mesh,
+            const Type& value,
+            const dimensionSet& dims,
+            const wordList& patchFieldTypes,
+            const wordList& actualPatchTypes = wordList()
+        );
+
+        //- Return tmp field (NO_READ, NO_WRITE)
+        //- from name, mesh, field value, dimensions and patch field types.
+        //- [Takes current timeName from the mesh registry].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<GeometricField<Type, PatchField, GeoMesh>> New
         (
             const word& name,
@@ -469,10 +564,25 @@ public:
             const wordList& actualPatchTypes = wordList()
         );
 
-        //- Return tmp field from name, mesh, dimensioned\<Type\>
+        //- Return tmp field (NO_READ, NO_WRITE)
-        //- and patch type.
+        //- from name, mesh, dimensioned-type and patch type.
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- [Takes current timeName from the mesh registry].
-        //  current timeName from the mesh registry
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        static tmp<GeometricField<Type, PatchField, GeoMesh>> New
+        (
+            const word& name,
+            IOobjectOption::registerOption regOpt,
+            const Mesh& mesh,
+            const dimensioned<Type>& dt,
+            const word& patchFieldType = PatchField<Type>::calculatedType()
+        );
+
+        //- Return tmp field (NO_READ, NO_WRITE)
+        //- from name, mesh, dimensioned-type and patch type.
+        //- [Takes current timeName from the mesh registry].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<GeometricField<Type, PatchField, GeoMesh>> New
         (
             const word& name,
@@ -481,10 +591,26 @@ public:
             const word& patchFieldType = PatchField<Type>::calculatedType()
         );
 
-        //- Return tmp field from name, mesh, dimensioned\<Type\>
+        //- Return tmp field (NO_READ, NO_WRITE)
-        //- and patch types.
+        //- from name, mesh, dimensioned-type and patch field types.
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- [Takes current timeName from the mesh registry].
-        //  current timeName from the mesh registry
+        //  For LEGACY_REGISTER, registration is determined by
+        //  objectRegistry::is_cacheTemporaryObject().
+        static tmp<GeometricField<Type, PatchField, GeoMesh>> New
+        (
+            const word& name,
+            IOobjectOption::registerOption regOpt,
+            const Mesh& mesh,
+            const dimensioned<Type>& dt,
+            const wordList& patchFieldTypes,
+            const wordList& actualPatchTypes = wordList()
+        );
+
+        //- Return tmp field (NO_READ, NO_WRITE)
+        //- from name, mesh, dimensioned-type and patch field types.
+        //- [Takes current timeName from the mesh registry].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<GeometricField<Type, PatchField, GeoMesh>> New
         (
             const word& name,
@@ -494,18 +620,21 @@ public:
             const wordList& actualPatchTypes = wordList()
         );
 
-        //- Return renamed tmp field
+        //- Return renamed tmp field (NO_READ, NO_WRITE)
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- [Takes current timeName from the mesh registry].
-        //  current timeName from the mesh registry
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<GeometricField<Type, PatchField, GeoMesh>> New
         (
             const word& newName,
             const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
         );
 
-        //- Rename tmp field and reset patch field type
+        //- Return renamed tmp field (NO_READ, NO_WRITE)
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- with reset patch field type.
-        //  current timeName from the mesh registry
+        //- [Takes current timeName from the mesh registry].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<GeometricField<Type, PatchField, GeoMesh>> New
         (
             const word& newName,
@@ -513,9 +642,11 @@ public:
             const word& patchFieldType
         );
 
-        //- Rename tmp field and reset patch field types and return
+        //- Return renamed tmp field (NO_READ, NO_WRITE)
-        //  The field is NO_READ, NO_WRITE, unregistered and uses the
+        //- with reset patch field types.
-        //  current timeName from the mesh registry
+        //- [Takes instance from the field].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         static tmp<GeometricField<Type, PatchField, GeoMesh>> New
         (
             const word& newName,
@@ -524,10 +655,11 @@ public:
             const wordList& actualPatchTypes = wordList()
         );
 
-        //- Construct tmp field based on mesh/registry information from
+        //- Construct tmp field (NO_READ, NO_WRITE)
-        //- an existing field.
+        //- based on mesh/registry information from an existing field.
-        //  Created NO_READ, NO_WRITE, NO_REGISTER, using the instance
+        //- [Takes instance from the field].
-        //  from the field
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         template<class AnyType>
         static tmp<GeometricField<Type, PatchField, GeoMesh>> New
         (
@@ -537,10 +669,12 @@ public:
             const word& patchFieldType = PatchField<Type>::calculatedType()
         );
 
-        //- Construct tmp field based on mesh/registry information from
+        //- Construct tmp field (NO_READ, NO_WRITE)
-        //- an existing field and initialise with value.
+        //- based on mesh/registry information from an existing field.
-        //  Created NO_READ, NO_WRITE, NO_REGISTER, using the instance
+        //- and initialise with value.
-        //  from the field
+        //- [Takes instance from the field].
+        //  Registration/persistence determined by
+        //  objectRegistry::is_cacheTemporaryObject().
         template<class AnyType>
         static tmp<GeometricField<Type, PatchField, GeoMesh>> New
         (
@@ -610,15 +744,15 @@ public:
         //- Write-access to the time index of the field
         inline label& timeIndex() noexcept;
 
+        //- The number of old time fields stored
+        label nOldTimes() const noexcept;
+
         //- Store the old-time fields
         void storeOldTimes() const;
 
         //- Store the old-time field
         void storeOldTime() const;
 
-        //- Return the number of old time fields stored
-        label nOldTimes() const;
-
         //- Return old time field
         const GeometricField<Type, PatchField, GeoMesh>& oldTime() const;
 
@@ -632,6 +766,9 @@ public:
         //- Return previous iteration field
         const GeometricField<Type, PatchField, GeoMesh>& prevIter() const;
 
+        //- Remove old-time and prev-iter fields
+        void clearOldTimes();
+
         //- Correct boundary field
         void correctBoundaryConditions();
 
@@ -648,9 +785,6 @@ public:
             const direction
         ) const;
 
-        //- WriteData member function required by regIOobject
-        bool writeData(Ostream&) const;
-
         //- Return transpose (only if it is a tensor field)
         tmp<GeometricField<Type, PatchField, GeoMesh>> T() const;
 
@@ -774,6 +908,13 @@ public:
         bool operator!=(const dimensioned<Type>&) = delete;
 
 
+    // Write
+
+        //- The writeData function (required by regIOobject),
+        //- calls operator<<
+        bool writeData(Ostream& os) const;
+
+
     // Ostream Operators
 
         friend Ostream& operator<< <Type, PatchField, GeoMesh>

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricFieldNew.C

@@ -26,8 +26,76 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+template<class... Args>
+Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
+Foam::GeometricField<Type, PatchField, GeoMesh>::New_impl
+(
+    IOobjectOption::registerOption regOpt,
+    const word& name,
+    const Mesh& mesh,
+    Args&&... args
+)
+{
+    auto ptr = tmp<GeometricField<Type, PatchField, GeoMesh>>::New
+    (
+        IOobject
+        (
+            name,
+            mesh.thisDb().time().timeName(),
+            mesh.thisDb(),
+            IOobjectOption::NO_READ,
+            IOobjectOption::NO_WRITE,
+            IOobjectOption::NO_REGISTER
+        ),
+        mesh,
+        std::forward<Args>(args)...
+    );
+
+    if (IOobjectOption::REGISTER == regOpt)
+    {
+        ptr->checkIn();
+    }
+    else if
+    (
+        // LEGACY_REGISTER: detect if caching is desired
+        (IOobjectOption::LEGACY_REGISTER == regOpt)
+     && ptr->db().is_cacheTemporaryObject(ptr.get())
+    )
+    {
+        ptr.protect(true);
+        ptr->checkIn();
+    }
+    return ptr;
+}
+
+
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
+Foam::GeometricField<Type, PatchField, GeoMesh>::New
+(
+    const word& name,
+    IOobjectOption::registerOption regOpt,
+    const Mesh& mesh,
+    const dimensionSet& dims,
+    const word& patchFieldType
+)
+{
+    return GeometricField<Type, PatchField, GeoMesh>::New_impl
+    (
+        regOpt,
+        name,
+        mesh,
+        dims,
+        patchFieldType
+    );
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
 Foam::GeometricField<Type, PatchField, GeoMesh>::New
@@ -38,20 +106,10 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
     const word& patchFieldType
 )
 {
-    const bool caching = mesh.thisDb().cacheTemporaryObject(name);
+    return GeometricField<Type, PatchField, GeoMesh>::New_impl
-
-    return tmp<GeometricField<Type, PatchField, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
+        IOobjectOption::LEGACY_REGISTER,
-        IOobject
+        name,
-        (
-            name,
-            mesh.thisDb().time().timeName(),
-            mesh.thisDb(),
-            IOobjectOption::NO_READ,
-            IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
-        ),
         mesh,
         dims,
         patchFieldType
@@ -59,6 +117,30 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
+Foam::GeometricField<Type, PatchField, GeoMesh>::New
+(
+    const word& name,
+    IOobjectOption::registerOption regOpt,
+    const Mesh& mesh,
+    const dimensionSet& dims,
+    const Field<Type>& iField,
+    const word& patchFieldType
+)
+{
+    return GeometricField<Type, PatchField, GeoMesh>::New_impl
+    (
+        regOpt,
+        name,
+        mesh,
+        dims,
+        iField,
+        patchFieldType
+    );
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
 Foam::GeometricField<Type, PatchField, GeoMesh>::New
@@ -70,20 +152,10 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
     const word& patchFieldType
 )
 {
-    const bool caching = mesh.thisDb().cacheTemporaryObject(name);
+    return GeometricField<Type, PatchField, GeoMesh>::New_impl
-
-    return tmp<GeometricField<Type, PatchField, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
+        IOobjectOption::LEGACY_REGISTER,
-        IOobject
+        name,
-        (
-            name,
-            mesh.thisDb().time().timeName(),
-            mesh.thisDb(),
-            IOobjectOption::NO_READ,
-            IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
-        ),
         mesh,
         dims,
         iField,
@@ -92,6 +164,30 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
+Foam::GeometricField<Type, PatchField, GeoMesh>::New
+(
+    const word& name,
+    IOobjectOption::registerOption regOpt,
+    const Mesh& mesh,
+    const dimensionSet& dims,
+    Field<Type>&& iField,
+    const word& patchFieldType
+)
+{
+    return GeometricField<Type, PatchField, GeoMesh>::New_impl
+    (
+        regOpt,
+        name,
+        mesh,
+        dims,
+        std::move(iField),
+        patchFieldType
+    );
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
 Foam::GeometricField<Type, PatchField, GeoMesh>::New
@@ -103,20 +199,10 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
     const word& patchFieldType
 )
 {
-    const bool caching = mesh.thisDb().cacheTemporaryObject(name);
+    return GeometricField<Type, PatchField, GeoMesh>::New_impl
-
-    return tmp<GeometricField<Type, PatchField, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
+        IOobjectOption::LEGACY_REGISTER,
-        IOobject
+        name,
-        (
-            name,
-            mesh.thisDb().time().timeName(),
-            mesh.thisDb(),
-            IOobjectOption::NO_READ,
-            IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
-        ),
         mesh,
         dims,
         std::move(iField),
@@ -130,26 +216,17 @@ Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
 Foam::GeometricField<Type, PatchField, GeoMesh>::New
 (
     const word& name,
+    IOobjectOption::registerOption regOpt,
     const Mesh& mesh,
     const Type& value,
     const dimensionSet& dims,
     const word& patchFieldType
 )
 {
-    const bool caching = mesh.thisDb().cacheTemporaryObject(name);
+    return GeometricField<Type, PatchField, GeoMesh>::New_impl
-
-    return tmp<GeometricField<Type, PatchField, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
+        regOpt,
-        IOobject
+        name,
-        (
-            name,
-            mesh.thisDb().time().timeName(),
-            mesh.thisDb(),
-            IOobjectOption::NO_READ,
-            IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
-        ),
         mesh,
         value,
         dims,
@@ -166,24 +243,38 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
     const Mesh& mesh,
     const Type& value,
     const dimensionSet& dims,
+    const word& patchFieldType
+)
+{
+    return GeometricField<Type, PatchField, GeoMesh>::New_impl
+    (
+        IOobjectOption::LEGACY_REGISTER,
+        name,
+        mesh,
+        value,
+        dims,
+        patchFieldType
+    );
+}
+
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
+Foam::GeometricField<Type, PatchField, GeoMesh>::New
+(
+    const word& name,
+    IOobjectOption::registerOption regOpt,
+    const Mesh& mesh,
+    const Type& value,
+    const dimensionSet& dims,
     const wordList& patchFieldTypes,
     const wordList& actualPatchTypes
 )
 {
-    const bool caching = mesh.thisDb().cacheTemporaryObject(name);
+    return GeometricField<Type, PatchField, GeoMesh>::New_impl
-
-    return tmp<GeometricField<Type, PatchField, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
+        regOpt,
-        IOobject
+        name,
-        (
-            name,
-            mesh.thisDb().time().timeName(),
-            mesh.thisDb(),
-            IOobjectOption::NO_READ,
-            IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
-        ),
         mesh,
         value,
         dims,
@@ -193,6 +284,54 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
+Foam::GeometricField<Type, PatchField, GeoMesh>::New
+(
+    const word& name,
+    const Mesh& mesh,
+    const Type& value,
+    const dimensionSet& dims,
+    const wordList& patchFieldTypes,
+    const wordList& actualPatchTypes
+)
+{
+    return GeometricField<Type, PatchField, GeoMesh>::New_impl
+    (
+        IOobjectOption::LEGACY_REGISTER,
+        name,
+        mesh,
+        value,
+        dims,
+        patchFieldTypes,
+        actualPatchTypes
+    );
+}
+
+
+template<class Type, template<class> class PatchField, class GeoMesh>
+Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
+Foam::GeometricField<Type, PatchField, GeoMesh>::New
+(
+    const word& name,
+    IOobjectOption::registerOption regOpt,
+    const Mesh& mesh,
+    const dimensioned<Type>& dt,
+    const word& patchFieldType
+)
+{
+    return GeometricField<Type, PatchField, GeoMesh>::New
+    (
+        name,
+        regOpt,
+        mesh,
+        dt.value(),
+        dt.dimensions(),
+        patchFieldType
+    );
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
 Foam::GeometricField<Type, PatchField, GeoMesh>::New
@@ -214,6 +353,31 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
+Foam::GeometricField<Type, PatchField, GeoMesh>::New
+(
+    const word& name,
+    IOobjectOption::registerOption regOpt,
+    const Mesh& mesh,
+    const dimensioned<Type>& dt,
+    const wordList& patchFieldTypes,
+    const wordList& actualPatchTypes
+)
+{
+    return GeometricField<Type, PatchField, GeoMesh>::New
+    (
+        name,
+        regOpt,
+        mesh,
+        dt.value(),
+        dt.dimensions(),
+        patchFieldTypes,
+        actualPatchTypes
+    );
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh>>
 Foam::GeometricField<Type, PatchField, GeoMesh>::New
@@ -246,11 +410,8 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
     const word& patchFieldType
 )
 {
-    const bool caching = tgf().db().cacheTemporaryObject(name);
+    auto ptr = tmp<GeometricField<Type, PatchField, GeoMesh>>::New
-
-    return tmp<GeometricField<Type, PatchField, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
         IOobject
         (
             name,
@@ -259,11 +420,21 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
             tgf().db(),
             IOobjectOption::NO_READ,
             IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
+            IOobjectOption::NO_REGISTER
         ),
         tgf,
         patchFieldType
     );
+
+    if
+    (
+        ptr->db().is_cacheTemporaryObject(ptr.get())
+    )
+    {
+        ptr.protect(true);
+        ptr->checkIn();
+    }
+    return ptr;
 }
 
 
@@ -275,11 +446,8 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
     const tmp<GeometricField<Type, PatchField, GeoMesh>>& tgf
 )
 {
-    const bool caching = tgf().db().cacheTemporaryObject(name);
+    auto ptr = tmp<GeometricField<Type, PatchField, GeoMesh>>::New
-
-    return tmp<GeometricField<Type, PatchField, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
         IOobject
         (
             name,
@@ -288,10 +456,20 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
             tgf().db(),
             IOobjectOption::NO_READ,
             IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
+            IOobjectOption::NO_REGISTER
         ),
         tgf
     );
+
+    if
+    (
+        ptr->db().is_cacheTemporaryObject(ptr.get())
+    )
+    {
+        ptr.protect(true);
+        ptr->checkIn();
+    }
+    return ptr;
 }
 
 
@@ -305,11 +483,8 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
     const wordList& actualPatchTypes
 )
 {
-    const bool caching = tgf().db().cacheTemporaryObject(name);
+    auto ptr = tmp<GeometricField<Type, PatchField, GeoMesh>>::New
-
-    return tmp<GeometricField<Type, PatchField, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
         IOobject
         (
             name,
@@ -318,12 +493,22 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
             tgf().db(),
             IOobjectOption::NO_READ,
             IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
+            IOobjectOption::NO_REGISTER
         ),
         tgf,
         patchFieldTypes,
         actualPatchTypes
     );
+
+    if
+    (
+        ptr->db().is_cacheTemporaryObject(ptr.get())
+    )
+    {
+        ptr.protect(true);
+        ptr->checkIn();
+    }
+    return ptr;
 }
 
 
@@ -338,11 +523,8 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
     const word& patchFieldType
 )
 {
-    const bool caching = fld.db().cacheTemporaryObject(name);
+    auto ptr = tmp<GeometricField<Type, PatchField, GeoMesh>>::New
-
-    return tmp<GeometricField<Type, PatchField, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
         IOobject
         (
             name,
@@ -350,12 +532,22 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
             fld.db(),
             IOobjectOption::NO_READ,
             IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
+            IOobjectOption::NO_REGISTER
         ),
         fld.mesh(),
         dims,
         patchFieldType
     );
+
+    if
+    (
+        ptr->db().is_cacheTemporaryObject(ptr.get())
+    )
+    {
+        ptr.protect(true);
+        ptr->checkIn();
+    }
+    return ptr;
 }
 
 
@@ -370,11 +562,8 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
     const word& patchFieldType
 )
 {
-    const bool caching = fld.db().cacheTemporaryObject(name);
+    auto ptr = tmp<GeometricField<Type, PatchField, GeoMesh>>::New
-
-    return tmp<GeometricField<Type, PatchField, GeoMesh>>::NewImmovable
     (
-        caching,        // (true: immovable, false: movable)
         IOobject
         (
             name,
@@ -382,13 +571,23 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::New
             fld.db(),
             IOobjectOption::NO_READ,
             IOobjectOption::NO_WRITE,
-            caching     // (true: REGISTER, false: NO_REGISTER)
+            IOobjectOption::NO_REGISTER
         ),
         fld.mesh(),
         dt.value(),
         dt.dimensions(),
         patchFieldType
     );
+
+    if
+    (
+        ptr->db().is_cacheTemporaryObject(ptr.get())
+    )
+    {
+        ptr.protect(true);
+        ptr->checkIn();
+    }
+    return ptr;
 }
 
 

src/OpenFOAM/fields/cloud/cloud.C

@@ -37,6 +37,7 @@ namespace Foam
 }
 
 const Foam::word Foam::cloud::prefix("lagrangian");
+
 Foam::word Foam::cloud::defaultName("defaultCloud");
 
 const Foam::Enum<Foam::cloud::geometryType>
@@ -51,7 +52,7 @@ Foam::cloud::geometryTypeNames
 
 Foam::cloud::cloud(const objectRegistry& obr)
 :
-    cloud(obr, defaultName)
+    cloud(obr, cloud::defaultName)
 {}
 
 
@@ -63,10 +64,11 @@ Foam::cloud::cloud(const objectRegistry& obr, const word& cloudName)
         (
             cloudName,
             obr.time().timeName(),
-            prefix,
+            cloud::prefix,
             obr,
             IOobject::NO_READ,
-            IOobject::AUTO_WRITE
+            IOobject::AUTO_WRITE,
+            IOobject::REGISTER
         )
     )
 {}

src/OpenFOAM/fields/cloud/cloud.H

@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef cloud_H
+#ifndef Foam_cloud_H
-#define cloud_H
+#define Foam_cloud_H
 
 #include "objectRegistry.H"
 #include "Enum.H"
@@ -59,16 +59,6 @@ class cloud
 :
     public objectRegistry
 {
-
-    // Private Member Functions
-
-        //- No copy construct
-        cloud(const cloud&) = delete;
-
-        //- No copy assignment
-        void operator=(const cloud&) = delete;
-
-
 public:
 
     //- Cloud geometry type (internal or IO representations)
@@ -78,11 +68,11 @@ public:
         POSITIONS       //!< positions
     };
 
+    //- Named enumerations "coordinates", "positions"
     static const Enum<geometryType> geometryTypeNames;
 
 
-    //- Runtime type information
+    // Static Data Members
-    TypeName("cloud");
 
         //- The prefix to local: %lagrangian
         static const word prefix;
@@ -91,12 +81,25 @@ public:
         static word defaultName;
 
 
+    //- Runtime type information
+    TypeName("cloud");
+
+
+    // Generated Methods
+
+        //- No copy construct
+        cloud(const cloud&) = delete;
+
+        //- No copy assignment
+        void operator=(const cloud&) = delete;
+
+
     // Constructors
 
-        //- Construct for given objectRegistry and default cloud name
+        //- Construct for given objectRegistry and efault cloud name
         explicit cloud(const objectRegistry& obr);
 
-        //- Construct for given objectRegistry and named cloud instance
+        //- Construct for given objectRegistry and named cloud
         cloud(const objectRegistry& obr, const word& cloudName);
 
 
@@ -139,7 +142,7 @@ public:
             //- Locate an IOField within object registry
             //  \return nullptr if not found or wrong type
             template<class Type>
-            inline static const IOField<Type>* findIOField
+            static const IOField<Type>* findIOField
             (
                 const word& fieldName,
                 const objectRegistry& obr
@@ -150,7 +153,7 @@ public:
 
             //- Locate the "position" IOField within object registry
             //  \return nullptr if not found or wrong type
-            inline static const IOField<point>* findIOPosition
+            static const IOField<point>* findIOPosition
             (
                 const objectRegistry& obr
             )
@@ -161,7 +164,7 @@ public:
             //- Lookup an IOField within object registry
             //  Fatal if not found or wrong type
             template<class Type>
-            inline static const IOField<Type>& lookupIOField
+            static const IOField<Type>& lookupIOField
             (
                 const word& fieldName,
                 const objectRegistry& obr

src/OpenFOAM/global/fileOperations/collatedFileOperation/OFstreamCollator.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2017-2018 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -361,10 +361,13 @@ bool Foam::OFstreamCollator::write
     off_t totalSize = 0;
     label maxLocalSize = 0;
     {
-        for (const label recvSize : recvSizes)
+        if (UPstream::master(localComm_))
         {
-            totalSize += recvSize;
+            for (const label recvSize : recvSizes)
-            maxLocalSize = max(maxLocalSize, recvSize);
+            {
+                totalSize += recvSize;
+                maxLocalSize = max(maxLocalSize, recvSize);
+            }
         }
         Pstream::broadcasts(localComm_, totalSize, maxLocalSize);
     }

src/OpenFOAM/global/fileOperations/masterUncollatedFileOperation/masterUncollatedFileOperation.C

@@ -2007,6 +2007,8 @@ Foam::fileOperations::masterUncollatedFileOperation::readStream
             << " fName : " << fName << " readOnProc:" << readOnProc << endl;
     }
 
+    // Close old stream
+    io.close();
 
     autoPtr<ISstream> isPtr;
     bool isCollated = false;
@@ -2044,7 +2046,9 @@ Foam::fileOperations::masterUncollatedFileOperation::readStream
                             << " doing straight IFstream input from "
                             << fName << endl;
                     }
-                    io = headerIO;
+
+                    // Copy IOobject header information
+                    io.IOobject::operator=(headerIO);
                     return isPtr;
                 }
             }

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

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,6 +30,7 @@ License
 #include "GAMGInterface.H"
 #include "processorGAMGInterface.H"
 #include "cyclicLduInterface.H"
+#include "PrecisionAdaptor.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -409,12 +410,16 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
 
     if (debug & 2)
     {
+        const auto& coarseAddr = meshLevels_[fineLevelIndex].lduAddr();
+
         Pout<< "GAMGAgglomeration :"
             << " agglomerated level " << fineLevelIndex
             << " from nCells:" << fineMeshAddr.size()
             << " nFaces:" << upperAddr.size()
             << " to nCells:" << nCoarseCells
-            << " nFaces:" << nCoarseFaces
+            << " nFaces:" << nCoarseFaces << nl
+            << "    lower:" << flatOutput(coarseAddr.lowerAddr()) << nl
+            << "    upper:" << flatOutput(coarseAddr.upperAddr()) << nl
             << endl;
     }
 }
@@ -430,13 +435,29 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
     const label levelIndex
 )
 {
-    const lduMesh& myMesh = meshLevels_[levelIndex-1];
+    // - Assemble all the procIDs in meshComm onto a single master
+    //   (procIDs[0]). This constructs a new communicator ('comm') first.
+    // - The master communicates with neighbouring masters using
+    //   allMeshComm
 
+    const lduMesh& myMesh = meshLevels_[levelIndex-1];
+    const label nOldInterfaces = myMesh.interfaces().size();
 
     procAgglomMap_.set(levelIndex, new labelList(procAgglomMap));
     agglomProcIDs_.set(levelIndex, new labelList(procIDs));
     procCommunicator_[levelIndex] = allMeshComm;
 
+    procAgglomCommunicator_.set
+    (
+        levelIndex,
+        new UPstream::communicator
+        (
+            meshComm,
+            procIDs
+        )
+    );
+    const label comm = agglomCommunicator(levelIndex);
+
     // These could only be set on the master procs but it is
     // quite convenient to also have them on the slaves
     procCellOffsets_.set(levelIndex, new labelList(0));
@@ -447,7 +468,7 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
 
     // Collect meshes
     PtrList<lduPrimitiveMesh> otherMeshes;
-    lduPrimitiveMesh::gather(meshComm, myMesh, procIDs, otherMeshes);
+    lduPrimitiveMesh::gather(comm, myMesh, otherMeshes);
 
     if (Pstream::myProcNo(meshComm) == procIDs[0])
     {
@@ -476,6 +497,44 @@ void Foam::GAMGAgglomeration::procAgglomerateLduAddressing
     }
 
 
+    // Scatter the procBoundaryMap back to the originating processor
+    // so it knows which proc boundaries are to be kept. This is used
+    // so we only send over interfaceFields on kept processors (see
+    // GAMGSolver::procAgglomerateMatrix)
+    // TBD: using sub-communicator here (instead of explicit procIDs). Should
+    //      use sub-communicators more in other places.
+    {
+        const CompactListList<label> data
+        (
+            CompactListList<label>::pack<labelList>
+            (
+                procBoundaryMap_[levelIndex]
+            )
+        );
+        const labelList localSizes = data.localSizes();
+        const labelList& localStarts = data.offsets();
+
+        // Make space
+        procBoundaryMap_[levelIndex].setSize(procIDs.size());
+        labelList& bMap = procBoundaryMap_[levelIndex][Pstream::myProcNo(comm)];
+        bMap.setSize(nOldInterfaces);
+
+        // Scatter relevant section to originating processor
+        UPstream::scatter
+        (
+            data.values().cdata(),
+
+            // Pass as List<int> for MPI
+            ConstPrecisionAdaptor<int, label, List>(localSizes).cref(),
+            ConstPrecisionAdaptor<int, label, List>(localStarts).cref(),
+
+            bMap.data(),
+            bMap.size(),
+            comm
+        );
+    }
+
+
     // Combine restrict addressing
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

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

@@ -69,6 +69,7 @@ void Foam::GAMGAgglomeration::compactLevels
     {
         procAgglomMap_.setSize(nCreatedLevels);
         agglomProcIDs_.setSize(nCreatedLevels);
+        procAgglomCommunicator_.setSize(nCreatedLevels);
         procCellOffsets_.setSize(nCreatedLevels);
         procFaceMap_.setSize(nCreatedLevels);
         procBoundaryMap_.setSize(nCreatedLevels);
@@ -289,6 +290,7 @@ Foam::GAMGAgglomeration::GAMGAgglomeration
     {
         procAgglomMap_.setSize(maxLevels_);
         agglomProcIDs_.setSize(maxLevels_);
+        procAgglomCommunicator_.setSize(maxLevels_);
         procCellOffsets_.setSize(maxLevels_);
         procFaceMap_.setSize(maxLevels_);
         procBoundaryMap_.setSize(maxLevels_);
@@ -569,6 +571,12 @@ Foam::label Foam::GAMGAgglomeration::procCommunicator(const label leveli) const
 }
 
 
+Foam::label Foam::GAMGAgglomeration::agglomCommunicator(const label leveli) const
+{
+    return procAgglomCommunicator_[leveli];
+}
+
+
 const Foam::labelList& Foam::GAMGAgglomeration::cellOffsets
 (
     const label leveli

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

@@ -131,6 +131,11 @@ protected:
             //- Communicator for given level
             mutable labelList procCommunicator_;
 
+            //- Communicator for collecting contributions. Note self-contained
+            //  communicator for ease of cleanup - this one gets allocated
+            //  internally, not by the GAMGProcAgglomerations
+            mutable PtrList<UPstream::communicator> procAgglomCommunicator_;
+
             //- Mapping from processor to procMeshLevel cells
             mutable PtrList<labelList> procCellOffsets_;
 
@@ -423,6 +428,17 @@ public:
                 const bool procAgglom
             ) const;
 
+            //- Prolong (interpolate by injection) cell field. Return reference
+            //  to (potentially scattered) coarse field
+            template<class Type>
+            const Field<Type>& prolongField
+            (
+                Field<Type>& ff,        // fine-level field
+                Field<Type>& work,      // work storage
+                const Field<Type>& cf,  // coarse-level field
+                const label coarseLevelIndex
+            ) const;
+
 
         // Processor agglomeration. Note that the mesh and agglomeration is
         // stored per fineLevel (even though it is the coarse level mesh that
@@ -466,6 +482,9 @@ public:
             //- Communicator for current level or -1
             label procCommunicator(const label fineLeveli) const;
 
+            //- Communicator for collecting contributions
+            label agglomCommunicator(const label fineLeveli) const;
+
             //- Mapping from processor to procMesh cells
             const labelList& cellOffsets(const label fineLeveli) const;
 

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

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -176,4 +177,56 @@ void Foam::GAMGAgglomeration::prolongField
 }
 
 
+template<class Type>
+const Foam::Field<Type>& Foam::GAMGAgglomeration::prolongField
+(
+    Field<Type>& ff,
+    Field<Type>& allCf,      // work storage
+    const Field<Type>& cf,
+    const label levelIndex
+) const
+{
+    const labelList& fineToCoarse = restrictAddressing_[levelIndex];
+
+    const label coarseLevelIndex = levelIndex+1;
+
+    if (hasProcMesh(coarseLevelIndex))
+    {
+        const label coarseComm =
+            UPstream::parent(procCommunicator_[coarseLevelIndex]);
+
+        const List<label>& procIDs = agglomProcIDs(coarseLevelIndex);
+        const labelList& offsets = cellOffsets(coarseLevelIndex);
+
+        const label localSize = nCells_[levelIndex];
+        allCf.resize_nocopy(localSize);
+
+        globalIndex::scatter
+        (
+            offsets,
+            coarseComm,
+            procIDs,
+            cf,
+            allCf,
+            UPstream::msgType(),
+            Pstream::commsTypes::nonBlocking    //Pstream::commsTypes::scheduled
+        );
+
+        forAll(fineToCoarse, i)
+        {
+            ff[i] = allCf[fineToCoarse[i]];
+        }
+        return allCf;
+    }
+    else
+    {
+        forAll(fineToCoarse, i)
+        {
+            ff[i] = cf[fineToCoarse[i]];
+        }
+        return cf;
+    }
+}
+
+
 // ************************************************************************* //

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

@@ -125,12 +125,12 @@ Foam::GAMGSolver::GAMGSolver
                     {
                         if (fineMeshInterfaces.set(intI))
                         {
-                            OStringStream os;
+                            OStringStream os(IOstreamOption::BINARY);
                             refCast<const GAMGInterface>
                             (
                                 fineMeshInterfaces[intI]
                             ).write(os);
-                            IStringStream is(os.str());
+                            IStringStream is(os.str(), IOstreamOption::BINARY);
 
                             dummyPrimMeshInterfaces.set
                             (
@@ -238,6 +238,8 @@ Foam::GAMGSolver::GAMGSolver
                     {
                         Pout<< "        " << i
                             << "\ttype:" << interfaces[i].type()
+                            << "\tsize:"
+                            << interfaces[i].interface().faceCells().size()
                             << endl;
                     }
                 }

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,7 +45,8 @@ Description
       - Coarse matrix scaling: performed by correction scaling, using steepest
         descent optimisation.
       - Type of cycle: V-cycle with optional pre-smoothing.
-      - Coarsest-level matrix solved using PCG or PBiCGStab.
+      - Coarsest-level matrix solved using any lduSolver (PCG, PBiCGStab,
+        smoothSolver) or direct solver on master processor
 
 SourceFiles
     GAMGSolver.C
@@ -190,9 +191,8 @@ class GAMGSolver
         //- Collect matrices from other processors
         void gatherMatrices
         (
-            const labelList& procIDs,
+            const label destLevel,
-            const lduMesh& dummyMesh,
+            const label comm,
-            const label meshComm,
 
             const lduMatrix& mat,
             const FieldField<Field, scalar>& interfaceBouCoeffs,
@@ -202,8 +202,7 @@ class GAMGSolver
             PtrList<lduMatrix>& otherMats,
             PtrList<FieldField<Field, scalar>>& otherBouCoeffs,
             PtrList<FieldField<Field, scalar>>& otherIntCoeffs,
-            List<boolList>& otherTransforms,
+            PtrList<PtrList<lduInterfaceField>>& otherInterfaces
-            List<List<label>>& otherRanks
         ) const;
 
         //- Agglomerate processor matrices

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

@@ -283,74 +283,136 @@ void Foam::GAMGSolver::agglomerateInterfaceCoefficients
 
 void Foam::GAMGSolver::gatherMatrices
 (
-    const labelList& procIDs,
+    const label destLevel,
-    const lduMesh& dummyMesh,
+    const label comm,
-    const label meshComm,
 
+    // Local matrix
     const lduMatrix& mat,
     const FieldField<Field, scalar>& interfaceBouCoeffs,
     const FieldField<Field, scalar>& interfaceIntCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
 
+    // Remote matrices
     PtrList<lduMatrix>& otherMats,
     PtrList<FieldField<Field, scalar>>& otherBouCoeffs,
     PtrList<FieldField<Field, scalar>>& otherIntCoeffs,
-    List<boolList>& otherTransforms,
+    PtrList<PtrList<lduInterfaceField>>& otherInterfaces
-    List<List<label>>& otherRanks
 ) const
 {
     if (debug & 2)
     {
+        const auto& procIDs = UPstream::procID(comm);
+
         Pout<< "GAMGSolver::gatherMatrices :"
             << " collecting matrices from procs:" << procIDs
-            << " using comm:" << meshComm << endl;
+            << " using comm:" << comm << endl;
     }
 
-    if (Pstream::myProcNo(meshComm) == procIDs[0])
+    const auto& boundaryMap = agglomeration_.boundaryMap(destLevel);
+
+    // Use PstreamBuffers
+    PstreamBuffers pBufs
+    (
+        UPstream::commsTypes::nonBlocking,
+        UPstream::msgType(),
+        comm
+    );
+
+    // Send to master
+    if (!UPstream::master(comm))
     {
-        // Master.
+        // Mark valid interfaces
-        otherMats.setSize(procIDs.size()-1);
+        // -1   : not set
-        otherBouCoeffs.setSize(procIDs.size()-1);
+        // >= 0 : coupled interface (might also be unmerged processor boundary)
-        otherIntCoeffs.setSize(procIDs.size()-1);
+        //
-        otherTransforms.setSize(procIDs.size()-1);
+        // Note: most processor interfaces will disappear. Originally
-        otherRanks.setSize(procIDs.size()-1);
+        // we did not know which ones were kept but this is now stored
+        // on the boundaryMap (even on the slave processors). So we can
+        // already filter here and avoid sending across typeNames etc.
 
-        for (label proci = 1; proci < procIDs.size(); proci++)
+        const label proci = UPstream::myProcNo(comm);
+
+        labelList validInterface(interfaces.size(), -1);
+        forAll(interfaces, intI)
         {
-            label otherI = proci-1;
+            const label allIntI = boundaryMap[proci][intI];
+            if (interfaces.set(intI) && allIntI != -1)
+            {
+                validInterface[intI] = intI;
+            }
+        }
 
-            IPstream fromProc
+        UOPstream toMaster(UPstream::masterNo(), pBufs);
-            (
-                Pstream::commsTypes::scheduled,
-                procIDs[proci],
-                0,          // bufSize
-                Pstream::msgType(),
-                meshComm
-            );
 
-            otherMats.set(otherI, new lduMatrix(dummyMesh, fromProc));
+        toMaster<< mat << token::SPACE << validInterface;
+
+        forAll(validInterface, intI)
+        {
+            if (validInterface[intI] != -1)
+            {
+                const auto& interface = refCast<const GAMGInterfaceField>
+                (
+                    interfaces[intI]
+                );
+
+                toMaster
+                    << interfaceBouCoeffs[intI]
+                    << interfaceIntCoeffs[intI]
+                    << interface.type();
+                interface.write(toMaster);
+            }
+        }
+    }
+
+    // Wait for finish
+    pBufs.finishedGathers();
+
+    // Consume
+    if (UPstream::master(comm))
+    {
+        const label nProcs = UPstream::nProcs(comm);
+
+        const lduMesh& destMesh = agglomeration_.meshLevel(destLevel);
+        lduInterfacePtrsList destInterfaces = destMesh.interfaces();
+
+        // Master.
+        otherMats.setSize(nProcs-1);
+        otherBouCoeffs.setSize(nProcs-1);
+        otherIntCoeffs.setSize(nProcs-1);
+        otherInterfaces.setSize(nProcs-1);
+
+        for (const int proci : UPstream::subProcs(comm))
+        {
+            const label otherI = proci-1;
+
+            UIPstream fromProc(proci, pBufs);
+
+            otherMats.set(otherI, new lduMatrix(destMesh, fromProc));
 
             // Receive number of/valid interfaces
-            boolList& procTransforms = otherTransforms[otherI];
+            // >= 0 : remote interface index
-            List<label>& procRanks = otherRanks[otherI];
+            // -1   : invalid interface
+            const labelList validInterface(fromProc);
 
-            fromProc >> procTransforms;
-            fromProc >> procRanks;
-
-            // Size coefficients
             otherBouCoeffs.set
             (
                 otherI,
-                new FieldField<Field, scalar>(procRanks.size())
+                new FieldField<Field, scalar>(validInterface.size())
             );
             otherIntCoeffs.set
             (
                 otherI,
-                new FieldField<Field, scalar>(procRanks.size())
+                new FieldField<Field, scalar>(validInterface.size())
             );
-            forAll(procRanks, intI)
+            otherInterfaces.set
+            (
+                otherI,
+                new PtrList<lduInterfaceField>(validInterface.size())
+            );
+
+            forAll(validInterface, intI)
             {
-                if (procRanks[intI] != -1)
+                if (validInterface[intI] != -1)
                 {
                     otherBouCoeffs[otherI].set
                     (
@@ -362,49 +424,25 @@ void Foam::GAMGSolver::gatherMatrices
                         intI,
                         new scalarField(fromProc)
                     );
-                }
-            }
-        }
-    }
-    else
-    {
-        // Send to master
 
-        // Count valid interfaces
+                    const word coupleType(fromProc);
-        boolList procTransforms(interfaceBouCoeffs.size(), false);
+
-        List<label> procRanks(interfaceBouCoeffs.size(), -1);
+                    const label allIntI = boundaryMap[proci][intI];
-        forAll(interfaces, intI)
+
-        {
+                    otherInterfaces[otherI].set
-            if (interfaces.set(intI))
-            {
-                const processorLduInterfaceField& interface =
-                    refCast<const processorLduInterfaceField>
                     (
-                        interfaces[intI]
+                        intI,
+                        GAMGInterfaceField::New
+                        (
+                            coupleType,
+                            refCast<const GAMGInterface>
+                            (
+                                destInterfaces[allIntI]
+                            ),
+                            fromProc
+                        ).release()
                     );
-
+                }
-                procTransforms[intI] = interface.doTransform();
-                procRanks[intI] = interface.rank();
-            }
-        }
-
-        OPstream toMaster
-        (
-            Pstream::commsTypes::scheduled,
-            procIDs[0],
-            0,
-            Pstream::msgType(),
-            meshComm
-        );
-
-        toMaster << mat << procTransforms << procRanks;
-        forAll(procRanks, intI)
-        {
-            if (procRanks[intI] != -1)
-            {
-                toMaster
-                    << interfaceBouCoeffs[intI]
-                    << interfaceIntCoeffs[intI];
             }
         }
     }
@@ -434,9 +472,11 @@ void Foam::GAMGSolver::procAgglomerateMatrix
         interfaceLevelsBouCoeffs_[levelI];
     const FieldField<Field, scalar>& coarsestIntCoeffs =
         interfaceLevelsIntCoeffs_[levelI];
-    const lduMesh& coarsestMesh = coarsestMatrix.mesh();
 
-    label coarseComm = coarsestMesh.comm();
+    // Communicator containing all processors to combine (=agglomProcIDs).
+    // Result will be on master of communicator.
+    const label agglomComm = agglomeration_.agglomCommunicator(levelI+1);
+
 
     // Gather all matrix coefficients onto agglomProcIDs[0]
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -444,43 +484,29 @@ void Foam::GAMGSolver::procAgglomerateMatrix
     PtrList<lduMatrix> otherMats;
     PtrList<FieldField<Field, scalar>> otherBouCoeffs;
     PtrList<FieldField<Field, scalar>> otherIntCoeffs;
-    List<boolList> otherTransforms;
+    PtrList<PtrList<lduInterfaceField>> otherInterfaces;
-    List<List<label>> otherRanks;
     gatherMatrices
     (
-        agglomProcIDs,
+        levelI+1,                       // allMesh level (only on master)
-        coarsestMesh,
+        agglomComm,
-        coarseComm,
 
-        coarsestMatrix,
+        coarsestMatrix,                 // master before gathering
-        coarsestBouCoeffs,
+        coarsestBouCoeffs,              // master before gathering
-        coarsestIntCoeffs,
+        coarsestIntCoeffs,              // master before gathering
-        coarsestInterfaces,
+        coarsestInterfaces,             // master before gathering
 
-        otherMats,
+        otherMats,                      // slave info
-        otherBouCoeffs,
+        otherBouCoeffs,                 // slave info
-        otherIntCoeffs,
+        otherIntCoeffs,                 // slave info
-        otherTransforms,
+        otherInterfaces                 // slave info
-        otherRanks
     );
 
 
-    if (Pstream::myProcNo(coarseComm) == agglomProcIDs[0])
+    if (UPstream::master(agglomComm))
     {
         // Agglomerate all matrix
         // ~~~~~~~~~~~~~~~~~~~~~~
 
-        //Pout<< "Own matrix:" << coarsestMatrix.info() << endl;
-        //
-        //forAll(otherMats, i)
-        //{
-        //    Pout<< "** otherMats " << i << " "
-        //        << otherMats[i].info()
-        //        << endl;
-        //}
-        //Pout<< endl;
-
-
         const lduMesh& allMesh = agglomeration_.meshLevel(levelI+1);
         const labelList& cellOffsets = agglomeration_.cellOffsets(levelI+1);
         const labelListList& faceMap = agglomeration_.faceMap(levelI+1);
@@ -566,7 +592,8 @@ void Foam::GAMGSolver::procAgglomerateMatrix
             allInterfaceIntCoeffs.set(intI, new scalarField(size));
         }
 
-        labelList nBounFaces(allMeshInterfaces.size());
+        UPtrList<lduInterfaceField> otherFlds(0);
+
         forAll(boundaryMap, proci)
         {
             const FieldField<Field, scalar>& procBouCoeffs
@@ -582,6 +609,7 @@ void Foam::GAMGSolver::procAgglomerateMatrix
               : otherIntCoeffs[proci-1]
             );
 
+
             const labelList& bMap = boundaryMap[proci];
             forAll(bMap, procIntI)
             {
@@ -594,47 +622,76 @@ void Foam::GAMGSolver::procAgglomerateMatrix
 
                     if (!allInterfaces.set(allIntI))
                     {
-                        // Construct lduInterfaceField
+                        const GAMGInterface& intf = refCast<const GAMGInterface>
+                        (
+                            allMeshInterfaces[allIntI]
+                        );
 
-                        bool doTransform = false;
-                        int rank = -1;
                         if (proci == 0)
                         {
-                            const processorGAMGInterfaceField& procInt =
+                            // Clone my local interfaceField. Since it is from
-                                refCast
+                            // this processor it will already exist, even if it
-                                <
+                            // is a processor one.
-                                    const processorGAMGInterfaceField
+
-                                >
+                            const auto& ffld =
+                            refCast<const GAMGInterfaceField>
+                            (
+                                coarsestInterfaces[procIntI]
+                            );
+
+                            allPrimitiveInterfaces.set
+                            (
+                                allIntI,
+                                ffld.clone
                                 (
-                                    coarsestInterfaces[procIntI]
+                                    intf,
-                                );
+                                    otherFlds
-                            doTransform = procInt.doTransform();
+                                ).release()
-                            rank = procInt.rank();
+                            );
                         }
                         else
                         {
-                            doTransform =
+                            // Recreate a remote interfaceField
-                                otherTransforms[proci-1][procIntI];
+                            if (otherInterfaces[proci-1].set(procIntI))
-                            rank = otherRanks[proci-1][procIntI];
+                            {
-                        }
+                                const auto& ffld =
-
+                                refCast<const GAMGInterfaceField>
-                        allPrimitiveInterfaces.set
-                        (
-                            allIntI,
-                            GAMGInterfaceField::New
-                            (
-                                refCast<const GAMGInterface>
                                 (
-                                    allMeshInterfaces[allIntI]
+                                    otherInterfaces[proci-1][procIntI]
-                                ),
+                                );
-                                doTransform,
+
-                                rank
+                                allPrimitiveInterfaces.set
-                            ).ptr()
+                                (
-                        );
+                                    allIntI,
+                                    ffld.clone
+                                    (
+                                        intf,
+                                        otherFlds
+                                    ).release()
+                                );
+                            }
+                            else
+                            {
+                                // Recreate a default interfaceField with
+                                // sensible defaults.
+                                // Should not occur since all unmerged
+                                // processor interfaces get transferred.
+                                allPrimitiveInterfaces.set
+                                (
+                                    allIntI,
+                                    GAMGInterfaceField::New
+                                    (
+                                        intf,
+                                        false,  //doTransform,
+                                        0       //rank
+                                    ).ptr()
+                                );
+                            }
+                        }
                         allInterfaces.set
                         (
                             allIntI,
-                            &allPrimitiveInterfaces[allIntI]
+                            allPrimitiveInterfaces.get(allIntI)
                         );
                     }
 
@@ -669,7 +726,6 @@ void Foam::GAMGSolver::procAgglomerateMatrix
                     const scalarField& procBou = procBouCoeffs[procIntI];
                     const scalarField& procInt = procIntCoeffs[procIntI];
 
-
                     forAll(map, i)
                     {
                         if (map[i] >= 0)
@@ -703,6 +759,7 @@ void Foam::GAMGSolver::procAgglomerateMatrix
             }
         }
 
+
         //Pout<< "** Assembled allMatrix:" << allMatrix.info() << endl;
         //
         //forAll(allInterfaces, intI)
@@ -712,19 +769,18 @@ void Foam::GAMGSolver::procAgglomerateMatrix
         //        Pout<< "    patch:" << intI
         //            << " type:" << allInterfaces[intI].type()
         //            << " size:"
-        //            << allInterfaces[intI].interface().
+        //            << allInterfaces[intI].interface().faceCells().size()
-        //                faceCells().size()
         //            << endl;
         //
-        //        //const scalarField& bouCoeffs = allInterfaceBouCoeffs[intI];
+        //        const scalarField& bouCoeffs = allInterfaceBouCoeffs[intI];
-        //        //const scalarField& intCoeffs = allInterfaceIntCoeffs[intI];
+        //        const scalarField& intCoeffs = allInterfaceIntCoeffs[intI];
-        //        //forAll(bouCoeffs, facei)
+        //        forAll(bouCoeffs, facei)
-        //        //{
+        //        {
-        //        //    Pout<< "        " << facei
+        //            Pout<< "        " << facei
-        //        //        << "\tbou:" << bouCoeffs[facei]
+        //                << "\tbou:" << bouCoeffs[facei]
-        //        //        << "\tint:" << intCoeffs[facei]
+        //                << "\tint:" << intCoeffs[facei]
-        //        //        << endl;
+        //                << endl;
-        //        //}
+        //        }
         //    }
         //}
     }

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

@@ -8,6 +8,7 @@
     Copyright (C) 2011-2017 OpenFOAM Foundation
     Copyright (C) 2016-2021,2023 OpenCFD Ltd.
     Copyright (C) 2023 Huawei (Yu Ankun)
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -301,6 +302,9 @@ void Foam::GAMGSolver::Vcycle
 
     solveScalarField dummyField(0);
 
+    // Work storage for prolongation
+    solveScalarField work;
+
     for (label leveli = coarsestLevel - 1; leveli >= 0; leveli--)
     {
         if (coarseCorrFields.set(leveli))
@@ -321,16 +325,18 @@ void Foam::GAMGSolver::Vcycle
                 preSmoothedCoarseCorrField = coarseCorrFields[leveli];
             }
 
-            agglomeration_.prolongField
+
+            // Prolong correction to leveli
+            const auto& cf = agglomeration_.prolongField
             (
-                coarseCorrFields[leveli],
+                coarseCorrFields[leveli],   // current level
+                work,
                 (
                     coarseCorrFields.set(leveli + 1)
                   ? coarseCorrFields[leveli + 1]
                   : dummyField              // dummy value
                 ),
-                leveli + 1,
+                leveli + 1
-                true
             );
 
 
@@ -346,43 +352,21 @@ void Foam::GAMGSolver::Vcycle
                     ACf.operator const solveScalarField&()
                 );
 
-            if (interpolateCorrection_) //&& leveli < coarsestLevel - 2)
+            if (interpolateCorrection_)
             {
-                if
+                // Normal operation : have both coarse level and fine
+                // level. No processor agglomeration
+                interpolate
                 (
-                    coarseCorrFields.set(leveli+1)
+                    coarseCorrFields[leveli],
-                 && (
+                    ACfRef,
-                        matrixLevels_[leveli].mesh().comm()
+                    matrixLevels_[leveli],
-                     == matrixLevels_[leveli+1].mesh().comm()
+                    interfaceLevelsBouCoeffs_[leveli],
-                    )
+                    interfaceLevels_[leveli],
-                )
+                    agglomeration_.restrictAddressing(leveli + 1),
-                {
+                    cf,
-                    // Normal operation : have both coarse level and fine
+                    cmpt
-                    // level. No processor agglomeration
+                );
-                    interpolate
-                    (
-                        coarseCorrFields[leveli],
-                        ACfRef,
-                        matrixLevels_[leveli],
-                        interfaceLevelsBouCoeffs_[leveli],
-                        interfaceLevels_[leveli],
-                        agglomeration_.restrictAddressing(leveli + 1),
-                        coarseCorrFields[leveli + 1],
-                        cmpt
-                    );
-                }
-                else
-                {
-                    interpolate
-                    (
-                        coarseCorrFields[leveli],
-                        ACfRef,
-                        matrixLevels_[leveli],
-                        interfaceLevelsBouCoeffs_[leveli],
-                        interfaceLevels_[leveli],
-                        cmpt
-                    );
-                }
             }
 
             // Scale coarse-grid correction field

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/GAMGInterfaceField/GAMGInterfaceField.C

@@ -35,6 +35,7 @@ namespace Foam
     defineTypeNameAndDebug(GAMGInterfaceField, 0);
     defineRunTimeSelectionTable(GAMGInterfaceField, lduInterface);
     defineRunTimeSelectionTable(GAMGInterfaceField, lduInterfaceField);
+    defineRunTimeSelectionTable(GAMGInterfaceField, Istream);
 }
 
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/GAMGInterfaceField/GAMGInterfaceField.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -104,6 +104,18 @@ public:
             (GAMGCp, doTransform, rank)
         );
 
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            GAMGInterfaceField,
+            Istream,
+            (
+                const GAMGInterface& GAMGCp,
+                Istream& is
+            ),
+            (GAMGCp, is)
+        );
+
 
     // Selectors
 
@@ -124,6 +136,15 @@ public:
             const int rank
         );
 
+        //- Return a pointer to a new interface created on freestore given
+        //  the fine interface and stream
+        static autoPtr<GAMGInterfaceField> New
+        (
+            const word& patchFieldType,
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        );
+
 
     // Constructors
 
@@ -150,6 +171,28 @@ public:
             interface_(GAMGCp)
         {}
 
+        //- Construct from GAMG interface and fine level interface field
+        GAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        )
+        :
+            lduInterfaceField(GAMGCp),
+            interface_(GAMGCp)
+        {}
+
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterfaceField> clone
+        (
+            const GAMGInterface& GAMGCp,
+            const UPtrList<lduInterfaceField>& other // other
+        ) const = 0;
+
+
+    //- Destructor
+    virtual ~GAMGInterfaceField() = default;
+
 
     // Member Functions
 
@@ -161,6 +204,11 @@ public:
                 return interface_;
             }
 
+
+        // I/O
+
+            //- Write to stream
+            virtual void write(Ostream&) const = 0;
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/GAMGInterfaceField/GAMGInterfaceFieldNew.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2021,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -79,4 +79,27 @@ Foam::autoPtr<Foam::GAMGInterfaceField> Foam::GAMGInterfaceField::New
 }
 
 
+Foam::autoPtr<Foam::GAMGInterfaceField> Foam::GAMGInterfaceField::New
+(
+    const word& patchFieldType,
+    const GAMGInterface& GAMGCp,
+    Istream& is
+)
+{
+    auto* ctorPtr = IstreamConstructorTable(patchFieldType);
+
+    if (!ctorPtr)
+    {
+        FatalErrorInLookup
+        (
+            "GAMGInterfaceField",
+            patchFieldType,
+            *IstreamConstructorTablePtr_
+        ) << exit(FatalError);
+    }
+
+    return autoPtr<GAMGInterfaceField>(ctorPtr(GAMGCp, is));
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/cyclicGAMGInterfaceField/cyclicGAMGInterfaceField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,6 +47,12 @@ namespace Foam
         cyclicGAMGInterfaceField,
         lduInterfaceField
     );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        cyclicGAMGInterfaceField,
+        Istream
+    );
 
     // Add under name cyclicSlip
     addNamedToRunTimeSelectionTable
@@ -63,6 +69,13 @@ namespace Foam
         lduInterfaceField,
         cyclicSlip
     );
+    addNamedToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        cyclicGAMGInterfaceField,
+        Istream,
+        cyclicSlip
+    );
 }
 
 
@@ -101,6 +114,38 @@ Foam::cyclicGAMGInterfaceField::cyclicGAMGInterfaceField
 {}
 
 
+Foam::cyclicGAMGInterfaceField::cyclicGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    Istream& is
+)
+:
+    GAMGInterfaceField(GAMGCp, is),
+    cyclicInterface_(refCast<const cyclicGAMGInterface>(GAMGCp)),
+    doTransform_(readBool(is)),
+    rank_(readLabel(is))
+{}
+
+
+Foam::cyclicGAMGInterfaceField::cyclicGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    const lduInterfaceField& local,
+    const UPtrList<lduInterfaceField>& other
+)
+:
+    GAMGInterfaceField(GAMGCp, local),
+    cyclicInterface_(refCast<const cyclicGAMGInterface>(GAMGCp)),
+    doTransform_(false),
+    rank_(0)
+{
+    const auto& p = refCast<const cyclicLduInterfaceField>(local);
+
+    doTransform_ = p.doTransform();
+    rank_ = p.rank();
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::cyclicGAMGInterfaceField::updateInterfaceMatrix
@@ -133,4 +178,12 @@ void Foam::cyclicGAMGInterfaceField::updateInterfaceMatrix
 }
 
 
+void Foam::cyclicGAMGInterfaceField::write(Ostream& os) const
+{
+    //GAMGInterfaceField::write(os);
+    os  << token::SPACE << doTransform()
+        << token::SPACE << rank();
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/cyclicGAMGInterfaceField/cyclicGAMGInterfaceField.H

@@ -100,6 +100,39 @@ public:
             const int rank
         );
 
+        //- Construct from GAMG interface and Istream
+        cyclicGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        );
+
+        //- Construct from GAMG interface and local and remote fields
+        cyclicGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            const lduInterfaceField& local, // local
+            const UPtrList<lduInterfaceField>& other // other
+        );
+
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterfaceField> clone
+        (
+            const GAMGInterface& GAMGCp,
+            const UPtrList<lduInterfaceField>& other // other
+        ) const
+        {
+            return autoPtr<GAMGInterfaceField>
+            (
+                new cyclicGAMGInterfaceField
+                (
+                    GAMGCp,
+                    *this,      // local field
+                    other       // other fields
+                )
+            );
+        }
+
 
     //- Destructor
     virtual ~cyclicGAMGInterfaceField() = default;
@@ -157,6 +190,12 @@ public:
                 const direction cmpt,
                 const Pstream::commsTypes commsType
             ) const;
+
+
+        // I/O
+
+            //- Write to stream
+            virtual void write(Ostream&) const;
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorCyclicGAMGInterfaceField/processorCyclicGAMGInterfaceField.C

@@ -46,6 +46,12 @@ namespace Foam
         processorCyclicGAMGInterfaceField,
         lduInterfaceField
     );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        processorCyclicGAMGInterfaceField,
+        Istream
+    );
 }
 
 
@@ -72,6 +78,16 @@ Foam::processorCyclicGAMGInterfaceField::processorCyclicGAMGInterfaceField
 {}
 
 
+Foam::processorCyclicGAMGInterfaceField::processorCyclicGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    Istream& is
+)
+:
+    processorGAMGInterfaceField(GAMGCp, is)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::processorCyclicGAMGInterfaceField::~processorCyclicGAMGInterfaceField()

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorCyclicGAMGInterfaceField/processorCyclicGAMGInterfaceField.H

@@ -87,6 +87,24 @@ public:
             const int rank
         );
 
+        //- Construct from GAMG interface and Istream
+        processorCyclicGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        );
+
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterfaceField> clone
+        (
+            const GAMGInterface& GAMGCp,
+            const UPtrList<lduInterfaceField>& other
+        ) const
+        {
+            NotImplemented;
+            return autoPtr<GAMGInterfaceField>(nullptr);
+        }
+
 
     // Destructor
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.C

@@ -47,6 +47,12 @@ namespace Foam
         processorGAMGInterfaceField,
         lduInterfaceField
     );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        processorGAMGInterfaceField,
+        Istream
+    );
 }
 
 
@@ -88,6 +94,19 @@ Foam::processorGAMGInterfaceField::processorGAMGInterfaceField
 {}
 
 
+Foam::processorGAMGInterfaceField::processorGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    Istream& is
+)
+:
+    GAMGInterfaceField(GAMGCp, is),
+    procInterface_(refCast<const processorGAMGInterface>(GAMGCp)),
+    doTransform_(readBool(is)),
+    rank_(readLabel(is))
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
@@ -193,4 +212,12 @@ void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
 }
 
 
+void Foam::processorGAMGInterfaceField::write(Ostream& os) const
+{
+    //GAMGInterfaceField::write(os);
+    os  << token::SPACE << doTransform()
+        << token::SPACE << rank();
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.H

@@ -119,6 +119,31 @@ public:
             const int rank
         );
 
+        //- Construct from GAMG interface and Istream
+        processorGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            Istream& is
+        );
+
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterfaceField> clone
+        (
+            const GAMGInterface& patch,
+            const UPtrList<lduInterfaceField>& other
+        ) const
+        {
+            return autoPtr<GAMGInterfaceField>
+            (
+                new processorGAMGInterfaceField
+                (
+                    patch,
+                    doTransform_,
+                    rank_
+                )
+            );
+        }
+
 
     //- Destructor
     virtual ~processorGAMGInterfaceField() = default;
@@ -201,6 +226,12 @@ public:
             {
                 return rank_;
             }
+
+
+        // I/O
+
+            //- Write to stream
+            virtual void write(Ostream&) const;
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterface.H

@@ -197,6 +197,26 @@ public:
             Istream& is
         );
 
+        //- Construct by assembling and return a clone.
+        virtual autoPtr<GAMGInterface> clone
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            const labelList& interfaceMap,  // current to fine interface
+            const labelUList& faceCells,
+            const labelUList& faceRestrictAddresssing,
+            const labelUList& faceOffsets,
+            const lduInterfacePtrsList& allInterfaces,
+            const label coarseComm,
+            const label myProcNo,
+            const labelList& procAgglomMap
+        ) const
+        {
+            NotImplemented;
+            return nullptr;
+        }
+
+
 
     // Member Functions
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/cyclicGAMGInterface/cyclicGAMGInterface.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019,2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -162,6 +162,38 @@ Foam::cyclicGAMGInterface::cyclicGAMGInterface
 {}
 
 
+Foam::cyclicGAMGInterface::cyclicGAMGInterface
+(
+    const label index,
+    const lduInterfacePtrsList& coarseInterfaces,
+    const lduInterface& fineInterface,
+    const labelList& interfaceMap,
+    const labelUList& faceCells,
+    const labelUList& faceRestrictAddresssing,
+    const labelUList& faceOffsets,
+    const lduInterfacePtrsList& allInterfaces
+)
+:
+    GAMGInterface
+    (
+        index,
+        coarseInterfaces,
+        faceCells,
+        faceRestrictAddresssing
+    ),
+    neighbPatchID_
+    (
+        interfaceMap.find
+        (
+            refCast<const cyclicLduInterface>(fineInterface).neighbPatchID()
+        )
+    ),
+    owner_(refCast<const cyclicLduInterface>(fineInterface).owner()),
+    forwardT_(refCast<const cyclicLduInterface>(fineInterface).forwardT()),
+    reverseT_(refCast<const cyclicLduInterface>(fineInterface).reverseT())
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::labelField> Foam::cyclicGAMGInterface::internalFieldTransfer

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/cyclicGAMGInterface/cyclicGAMGInterface.H

@@ -108,6 +108,51 @@ public:
             Istream& is
         );
 
+        //- Construct from fine level interface,
+        //- local and neighbour restrict addressing
+        cyclicGAMGInterface
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            const lduInterface& fineInterface,
+            const labelList& interfaceMap,
+            const labelUList& faceCells,
+            const labelUList& faceRestrictAddresssing,
+            const labelUList& faceOffsets,
+            const lduInterfacePtrsList& allInterfaces
+        );
+
+        //- Construct by assembling and returning a clone.
+        virtual autoPtr<GAMGInterface> clone
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            const labelList& interfaceMap,
+            const labelUList& faceCells,
+            const labelUList& faceRestrictAddresssing,
+            const labelUList& faceOffsets,
+            const lduInterfacePtrsList& allInterfaces,
+            const label coarseComm,
+            const label myProcNo,
+            const labelList& procAgglomMap
+        ) const
+        {
+            return autoPtr<GAMGInterface>
+            (
+                new cyclicGAMGInterface
+                (
+                    index,
+                    coarseInterfaces,
+                    *this,
+                    interfaceMap,
+                    faceCells,
+                    faceRestrictAddresssing,
+                    faceOffsets,
+                    allInterfaces
+                )
+            );
+        }
+
 
     //- Destructor
     virtual ~cyclicGAMGInterface() = default;

src/OpenFOAM/memory/tmp/tmp.H

@@ -120,10 +120,6 @@ public:
         //- Construct, taking ownership of the pointer.
         inline explicit tmp(T* p);
 
-        //- Construct, taking ownership of the pointer,
-        //- with specified protection against moving.
-        inline tmp(bool immovable, T* p);
-
         //- Implicit construct for a const reference to an object.
         inline constexpr tmp(const T& obj) noexcept;
 
@@ -173,19 +169,6 @@ public:
             return tmp<T>(new T(std::forward<Args>(args)...));
         }
 
-        //- Construct movable/immovable tmp with forwarding arguments
-        //  \param immovable create as non-movable pointer
-        //  \param args list of arguments with which an instance of T
-        //      will be constructed.
-        //
-        //  \note Similar to std::make_shared, but the overload for
-        //      array types is not disabled.
-        template<class... Args>
-        static tmp<T> NewImmovable(bool immovable, Args&&... args)
-        {
-            return tmp<T>(immovable, new T(std::forward<Args>(args)...));
-        }
-
         //- Construct tmp from derived type with forwarding arguments
         //  \param args list of arguments with which an instance of U
         //      will be constructed.
@@ -256,16 +239,16 @@ public:
         //- delete object and set pointer to nullptr
         inline void clear() const noexcept;
 
+        //- Use specified protection against moving for the managed pointer.
+        //- No-op for references.
+        inline void protect(bool on) noexcept;
+
         //- Clear existing and transfer ownership.
         inline void reset(tmp<T>&& other) noexcept;
 
         //- Delete managed temporary object and set to new given pointer.
         inline void reset(T* p = nullptr) noexcept;
 
-        //- Delete managed temporary object and set to new given pointer,
-        //- with specified protection against moving.
-        inline void reset(bool immovable, T* p) noexcept;
-
         //- Avoid inadvertent casting (to object or pointer)
         void reset(const autoPtr<T>&) = delete;
 

src/OpenFOAM/memory/tmp/tmpI.H

@@ -90,15 +90,6 @@ inline Foam::tmp<T>::tmp(T* p)
 }
 
 
-template<class T>
-inline Foam::tmp<T>::tmp(bool immovable, T* p)
-:
-    tmp<T>(p)
-{
-    if (ptr_ && immovable) type_ = CACHE_PTR;
-}
-
-
 template<class T>
 inline constexpr Foam::tmp<T>::tmp(const T& obj) noexcept
 :
@@ -319,6 +310,28 @@ inline void Foam::tmp<T>::clear() const noexcept
 }
 
 
+template<class T>
+inline void Foam::tmp<T>::protect(bool on) noexcept
+{
+    if (on)
+    {
+        // ON: from PTR -> CACHE_PTR, but not nullptr
+        if (ptr_ && type_ == PTR)
+        {
+            type_ = CACHE_PTR;
+        }
+    }
+    else
+    {
+        // OFF: from CACHE_PTR -> PTR
+        if (type_ == CACHE_PTR)
+        {
+            type_ = PTR;
+        }
+    }
+}
+
+
 template<class T>
 inline void Foam::tmp<T>::reset(T* p) noexcept
 {
@@ -328,14 +341,6 @@ inline void Foam::tmp<T>::reset(T* p) noexcept
 }
 
 
-template<class T>
-inline void Foam::tmp<T>::reset(bool immovable, T* p) noexcept
-{
-    reset(p);
-    if (ptr_ && immovable) type_ = CACHE_PTR;
-}
-
-
 template<class T>
 inline void Foam::tmp<T>::reset(tmp<T>&& other) noexcept
 {

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C

@@ -28,6 +28,7 @@ License
 
 #include "lduPrimitiveMesh.H"
 #include "processorLduInterface.H"
+#include "processorCyclicGAMGInterface.H"
 #include "edgeHashes.H"
 #include "labelPair.H"
 #include "processorGAMGInterface.H"
@@ -308,6 +309,56 @@ Foam::labelList Foam::lduPrimitiveMesh::upperTriOrder
 }
 
 
+Foam::label Foam::lduPrimitiveMesh::findConnectedInterface
+(
+    const lduMesh& myMesh,
+    const PtrList<lduPrimitiveMesh>& otherMeshes,
+    const labelPairList& procAndInterfaces,
+
+    const label nbrProci,
+    const label myRank
+) const
+{
+    // Find mesh, interfacei in procAndInterfaces
+    label nbrInti = -1;
+    for (const auto& procAndInterface : procAndInterfaces)
+    {
+        const label proci = procAndInterface[0];
+
+        if (proci == nbrProci)
+        {
+            const label interfacei = procAndInterface[1];
+            const lduInterfacePtrsList interfaces =
+                mesh
+                (
+                    myMesh,
+                    otherMeshes,
+                    proci
+                ).interfaces();
+            const processorLduInterface& pldui =
+                refCast<const processorLduInterface>
+                (
+                    interfaces[interfacei]
+                );
+
+            if (pldui.neighbProcNo() == myRank)
+            {
+                nbrInti = procAndInterface[1];
+                break;
+            }
+        }
+    }
+
+
+    if (nbrInti == -1)
+    {
+        FatalErrorInFunction
+            << "procAndInterfaces:" << procAndInterfaces << abort(FatalError);
+    }
+    return nbrInti;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::lduPrimitiveMesh::lduPrimitiveMesh
@@ -472,7 +523,6 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     const label nMeshes = otherMeshes.size()+1;
 
     const label myAgglom = procAgglomMap[UPstream::myProcNo(currentComm)];
-
     if (lduPrimitiveMesh::debug)
     {
         Pout<< "I am " << UPstream::myProcNo(currentComm)
@@ -481,6 +531,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
             << endl;
     }
 
+    const lduInterfacePtrsList myInterfaces = myMesh.interfaces();
 
     forAll(procIDs, i)
     {
@@ -522,24 +573,28 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
           + procMesh.lduAddr().lowerAddr().size();
     }
 
-    // Count how faces get added. Interfaces inbetween get merged.
+    // Count how faces get added. Proc interfaces inbetween get merged.
 
-    // Merged interfaces: map from two coarse processors back to
+    // Merged proc interfaces: map from two coarse processors back to
     // - procMeshes
     // - interface in procMesh
     // (estimate size from number of patches of mesh0)
-    EdgeMap<labelPairList> mergedMap(2*myMesh.interfaces().size());
+    EdgeMap<labelPairList> mergedMap(2*myInterfaces.size());
 
-    // Unmerged interfaces: map from two coarse processors back to
+    // Unmerged proc interfaces: map from two coarse processors back to
     // - procMeshes
     // - interface in procMesh
     EdgeMap<labelPairList> unmergedMap(mergedMap.size());
 
+    // (unmerged) global interfaces. These are present on all processors
+    // in the same order (and keep the order in the merged mesh)
+    List<DynamicList<label>> procToGlobal(nMeshes);
+
+
     boundaryMap.setSize(nMeshes);
     boundaryFaceMap.setSize(nMeshes);
 
 
-    label nOtherInterfaces = 0;
     labelList nCoupledFaces(nMeshes, Zero);
 
     for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
@@ -560,6 +615,16 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
 
                 if (isA<processorLduInterface>(ldui))
                 {
+                    if (isA<processorCyclicGAMGInterface>(ldui))
+                    {
+                        FatalErrorInFunction
+                            << "At mesh from processor " << procIDs[procMeshI]
+                            << " have interface " << intI
+                            << " of unhandled type " << ldui.type()
+                            << ". Adapt decomposition to avoid these"
+                            << exit(FatalError);
+                    }
+
                     const processorLduInterface& pldui =
                         refCast<const processorLduInterface>(ldui);
 
@@ -582,7 +647,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                         // Merged interface
                         if (debug)
                         {
-                            Pout<< "merged interface: myProcNo:"
+                            Pout<< "merged proc interface: myProcNo:"
                                 << pldui.myProcNo()
                                 << " nbr:" << pldui.neighbProcNo()
                                 << " size:" << ldui.faceCells().size()
@@ -607,7 +672,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                     {
                         if (debug)
                         {
-                            Pout<< "external interface: myProcNo:"
+                            Pout<< "external proc interface: myProcNo:"
                                 << pldui.myProcNo()
                                 << " nbr:" << pldui.neighbProcNo()
                                 << " size:" << ldui.faceCells().size()
@@ -623,53 +688,80 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                 else
                 {
                     // Still external (non proc) interface
-                    FatalErrorInFunction
+                    procToGlobal[procMeshI].append(intI);
-                        << "At mesh from processor " << procIDs[procMeshI]
-                        << " have interface " << intI
-                        << " of unhandled type " << interfaces[intI].type()
-                        << exit(FatalError);
-
-                    ++nOtherInterfaces;
                 }
             }
         }
     }
 
 
-
+    // Check that all processors have any global patches in the same order
     if (debug)
     {
-        Pout<< "Remaining interfaces:" << endl;
+        const auto& global0 = procToGlobal[0];
-        forAllConstIters(unmergedMap, iter)
+        for (label procMeshI = 1; procMeshI < nMeshes; ++procMeshI)
         {
-            Pout<< "    agglom procEdge:" << iter.key() << endl;
+            const auto& global = procToGlobal[procMeshI];
-            const labelPairList& elems = iter.val();
+            if (global != global0)
-            forAll(elems, i)
             {
-                label procMeshI = elems[i][0];
+                FatalErrorInFunction
-                label interfacei = elems[i][1];
+                    << "At mesh from processor " << procIDs[procMeshI]
-                const lduInterfacePtrsList interfaces =
+                    << " have global interfaces " << global
-                    mesh(myMesh, otherMeshes, procMeshI).interfaces();
+                    << " which differ from those on processor "
-
+                    << procIDs[procMeshI]
-                const processorLduInterface& pldui =
+                    << " : " << global0
-                    refCast<const processorLduInterface>
+                    << exit(FatalError);
-                    (
-                        interfaces[interfacei]
-                    );
-
-                Pout<< "        proc:" << procIDs[procMeshI]
-                    << " interfacei:" << interfacei
-                    << " between:" << pldui.myProcNo()
-                    << " and:" << pldui.neighbProcNo()
-                    << endl;
             }
-            Pout<< endl;
         }
     }
+
+
     if (debug)
     {
+        Pout<< "Global interfaces:" << endl;
+        const auto& global0 = procToGlobal[0];
+        for (const label intI : global0)
+        {
+            Pout<< "        interfacei:" << intI
+                << " type:" << myInterfaces[intI].type()
+                << endl;
+
+        }
+        Pout<< endl;
+
+        Pout<< "Remaining interfaces:" << endl;
+        for (const auto& iter : unmergedMap.csorted())
+        {
+            Pout<< "    agglom procEdge:" << iter.key() << endl;
+            const labelPairList& elems = iter.val();
+            forAll(elems, i)
+            {
+                label procMeshI = elems[i][0];
+                label interfacei = elems[i][1];
+                const lduInterfacePtrsList interfaces =
+                    mesh(myMesh, otherMeshes, procMeshI).interfaces();
+
+                const processorLduInterface& pldui =
+                    refCast<const processorLduInterface>
+                    (
+                        interfaces[interfacei]
+                    );
+
+                Pout<< "        proc:" << procIDs[procMeshI]
+                    << " interfacei:" << interfacei
+                    << " between:" << pldui.myProcNo()
+                    << " and:" << pldui.neighbProcNo()
+                    << " localsize:"
+                    << interfaces[interfacei].faceCells().size()
+                    << endl;
+            }
+            Pout<< endl;
+        }
+        Pout<< endl;
+
         Pout<< "Merged interfaces:" << endl;
-        forAllConstIters(mergedMap, iter)
+        //forAllConstIters(mergedMap, iter)
+        for (const auto& iter : mergedMap.csorted())
         {
             Pout<< "    agglom procEdge:" << iter.key() << endl;
             const labelPairList& elems = iter.val();
@@ -690,10 +782,13 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                     << " interfacei:" << interfacei
                     << " between:" << pldui.myProcNo()
                     << " and:" << pldui.neighbProcNo()
+                    << " localsize:"
+                    << interfaces[interfacei].faceCells().size()
                     << endl;
             }
             Pout<< endl;
         }
+        Pout<< endl;
     }
 
 
@@ -780,43 +875,14 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
                         {
                             const labelPairList& elems = fnd();
 
-                            // Find nbrP in elems
+                            const label nbrIntI = findConnectedInterface
-                            label nbrIntI = -1;
+                            (
-                            forAll(elems, i)
+                                myMesh,
-                            {
+                                otherMeshes,
-                                label proci = elems[i][0];
+                                elems,
-                                label interfacei = elems[i][1];
+                                nbrProcMeshI,
-                                const lduInterfacePtrsList interfaces =
+                                procIDs[procMeshI]
-                                    mesh
+                            );
-                                    (
-                                        myMesh,
-                                        otherMeshes,
-                                        proci
-                                    ).interfaces();
-                                const processorLduInterface& pldui =
-                                    refCast<const processorLduInterface>
-                                    (
-                                        interfaces[interfacei]
-                                    );
-
-                                if
-                                (
-                                    elems[i][0] == nbrProcMeshI
-                                 && pldui.neighbProcNo() == procIDs[procMeshI]
-                                )
-                                {
-                                    nbrIntI = elems[i][1];
-                                    break;
-                                }
-                            }
-
-
-                            if (nbrIntI == -1)
-                            {
-                                FatalErrorInFunction
-                                    << "elems:" << elems << abort(FatalError);
-                            }
-
 
                             const lduInterfacePtrsList nbrInterfaces = mesh
                             (
@@ -930,12 +996,118 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     // Kept interfaces
     // ~~~~~~~~~~~~~~~
 
-    interfaces_.setSize(unmergedMap.size() + nOtherInterfaces);
+    interfaces_.setSize(unmergedMap.size() + procToGlobal[0].size());
     primitiveInterfaces_.setSize(interfaces_.size());
 
     label allInterfacei = 0;
 
-    forAllConstIters(unmergedMap, iter)
+
+    // Global interfaces
+    // ~~~~~~~~~~~~~~~~~
+    // (e.g. cyclicAMI)
+
+    {
+        const auto& global0 = procToGlobal[0];
+
+        for (const label interfacei : global0)
+        {
+            //Pout<< "        interfacei:" << interfacei
+            //    << " type:" << interfaces[interfacei].type()
+            //    << endl;
+
+            // Just add all individual face-cells in processor order
+
+            // Count
+            label n = 0;
+            for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
+            {
+                const auto& procMesh = mesh(myMesh, otherMeshes, procMeshI);
+                n += procMesh.interfaces()[interfacei].faceCells().size();
+            }
+
+            // Size
+            labelField allFaceCells(n);
+            labelField allFaceRestrictAddressing(n);
+            labelList faceOffsets(nMeshes+1, 0);
+            lduInterfacePtrsList allProcInterfaces(nMeshes);
+            n = 0;
+
+            // Fill
+            for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
+            {
+                faceOffsets[procMeshI] = n;
+
+                const auto& procMesh = mesh(myMesh, otherMeshes, procMeshI);
+                const lduInterfacePtrsList interfaces = procMesh.interfaces();
+
+                allProcInterfaces.set(procMeshI, &interfaces[interfacei]);
+                boundaryMap[procMeshI][interfacei] = allInterfacei;
+                labelList& bfMap = boundaryFaceMap[procMeshI][interfacei];
+
+                const labelUList& l = interfaces[interfacei].faceCells();
+                bfMap.setSize(l.size());
+
+                forAll(l, facei)
+                {
+                    allFaceCells[n] = cellOffsets[procMeshI]+l[facei];
+                    allFaceRestrictAddressing[n] = n;
+                    bfMap[facei] = n;
+                    n++;
+                }
+            }
+
+            // For convenience populate last element with size
+            faceOffsets.last() = n;
+
+
+            const auto& myFineInterface =
+                refCast<const GAMGInterface>(myInterfaces[interfacei]);
+
+            autoPtr<GAMGInterface> ppPtr
+            (
+                myFineInterface.clone
+                (
+                    allInterfacei,
+                    interfaces_,
+                    global0,
+                    allFaceCells,
+                    allFaceRestrictAddressing,
+                    faceOffsets,
+                    allProcInterfaces,
+                    comm_,
+                    myAgglom,
+                    procAgglomMap
+                )
+            );
+
+            primitiveInterfaces_.set
+            (
+                allInterfacei,
+                ppPtr.ptr()
+            );
+
+            interfaces_.set(allInterfacei, &primitiveInterfaces_[allInterfacei]);
+
+            if (debug)
+            {
+                Pout<< "Created " << interfaces_[allInterfacei].type()
+                    << " interface at " << allInterfacei
+                    << " comm:" << comm_
+                    << " myProcNo:" << myAgglom
+                    << " nFaces:" << allFaceCells.size()
+                    << endl;
+            }
+
+            allInterfacei++;
+        }
+    }
+
+
+    // Unmerged processor interfaces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    //forAllConstIters(unmergedMap, iter)
+    for (const auto& iter : unmergedMap.csorted())
     {
         const labelPairList& elems = iter.val();
 
@@ -1085,6 +1257,58 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     }
 
 
+    if (allInterfacei != interfaces_.size())
+    {
+        FatalErrorInFunction << "allInterfacei:" << allInterfacei
+            << " interfaces_:" << interfaces_.size() << exit(FatalError);
+    }
+
+
+    if (debug)
+    {
+        Pout<< endl
+            << "Created new lduPrimitiveMesh:" << nl
+            << "    cells:" << this->lduAddr().size() << nl
+            << "    internal face lower:"
+            << this->lduAddr().lowerAddr().size() << nl
+            << "    internal faces upper:"
+            << this->lduAddr().upperAddr().size()
+            << endl;
+        forAll(interfaces_, i)
+        {
+            if (interfaces_.set(i))
+            {
+                Pout<< "    interface:" << i << " type:" << interfaces_[i].type()
+                    << nl
+                    << "    faceCells:" << flatOutput(interfaces_[i].faceCells())
+                    << endl;
+            }
+        }
+
+
+        Pout<< "Original input meshes:" << endl;
+        forAll(boundaryMap, procMeshI)
+        {
+            const auto& procMesh = mesh(myMesh, otherMeshes, procMeshI);
+            const lduInterfacePtrsList interfaces = procMesh.interfaces();
+
+            Pout<< "    proc:" << procMeshI
+                << " interfaces:" << interfaces.size() << endl;
+            forAll(interfaces, inti)
+            {
+                if (interfaces.set(inti))
+                {
+                    Pout<< "        int:" << inti
+                        << " type:" << interfaces[inti].type()
+                        << " size:" << interfaces[inti].faceCells().size()
+                        << " maps to:" << boundaryMap[procMeshI][inti]
+                        << endl;
+                }
+            }
+        }
+    }
+
+
     patchSchedule_ = nonBlockingSchedule<processorGAMGInterface>(interfaces_);
 
     if (debug)
@@ -1111,89 +1335,25 @@ void Foam::lduPrimitiveMesh::gather
 (
     const label comm,
     const lduMesh& mesh,
-    const labelList& procIDs,
     PtrList<lduPrimitiveMesh>& otherMeshes
 )
 {
     // Force calculation of schedule (since does parallel comms)
     (void)mesh.lduAddr().patchSchedule();
 
-    if (Pstream::myProcNo(comm) == procIDs[0])
+    // Use PstreamBuffers
+    PstreamBuffers pBufs
+    (
+        Pstream::commsTypes::nonBlocking,
+        UPstream::msgType(),
+        comm
+    );
+
+    // Send to master
+    if (!Pstream::master(comm))
     {
-        // Master.
-        otherMeshes.setSize(procIDs.size()-1);
-
-        for (label i = 1; i < procIDs.size(); ++i)
-        {
-            //Pout<< "on master :"
-            //    << " receiving from proc " << procIDs[i] << endl;
-
-            IPstream fromProc
-            (
-                Pstream::commsTypes::scheduled,
-                procIDs[i],
-                0,          // bufSize
-                Pstream::msgType(),
-                comm
-            );
-
-            label nCells = readLabel(fromProc);
-            labelList lowerAddr(fromProc);
-            labelList upperAddr(fromProc);
-            boolList validInterface(fromProc);
-
-
-            // Construct mesh without interfaces
-            otherMeshes.set
-            (
-                i-1,
-                new lduPrimitiveMesh
-                (
-                    nCells,
-                    lowerAddr,
-                    upperAddr,
-                    comm,
-                    true    // reuse
-                )
-            );
-
-            // Construct GAMGInterfaces
-            lduInterfacePtrsList newInterfaces(validInterface.size());
-            forAll(validInterface, intI)
-            {
-                if (validInterface[intI])
-                {
-                    word coupleType(fromProc);
-
-                    newInterfaces.set
-                    (
-                        intI,
-                        GAMGInterface::New
-                        (
-                            coupleType,
-                            intI,
-                            otherMeshes[i-1].rawInterfaces(),
-                            fromProc
-                        ).ptr()
-                    );
-                }
-            }
-
-            otherMeshes[i-1].addInterfaces
-            (
-                newInterfaces,
-                nonBlockingSchedule<processorGAMGInterface>
-                (
-                    newInterfaces
-                )
-            );
-       }
-    }
-    else if (procIDs.found(Pstream::myProcNo(comm)))
-    {
-        // Send to master
-
         const lduAddressing& addressing = mesh.lduAddr();
+
         lduInterfacePtrsList interfaces(mesh.interfaces());
         boolList validInterface(interfaces.size());
         forAll(interfaces, intI)
@@ -1201,14 +1361,7 @@ void Foam::lduPrimitiveMesh::gather
             validInterface[intI] = interfaces.set(intI);
         }
 
-        OPstream toMaster
+        UOPstream toMaster(UPstream::masterNo(), pBufs);
-        (
-            Pstream::commsTypes::scheduled,
-            procIDs[0],
-            0,
-            Pstream::msgType(),
-            comm
-        );
 
         toMaster
             << addressing.size()
@@ -1230,6 +1383,74 @@ void Foam::lduPrimitiveMesh::gather
             }
         }
     }
+
+    // Wait for finish
+    pBufs.finishedGathers();
+
+    // Consume
+    if (Pstream::master(comm))
+    {
+        const label nProcs = UPstream::nProcs(comm);
+
+        // Master.
+        otherMeshes.setSize(nProcs-1);
+
+        for (const int proci : UPstream::subProcs(comm))
+        {
+            UIPstream fromProc(proci, pBufs);
+
+            const label nCells = readLabel(fromProc);
+            labelList lowerAddr(fromProc);
+            labelList upperAddr(fromProc);
+            const boolList validInterface(fromProc);
+
+
+            // Construct mesh without interfaces
+            otherMeshes.set
+            (
+                proci-1,
+                new lduPrimitiveMesh
+                (
+                    nCells,
+                    lowerAddr,
+                    upperAddr,
+                    mesh.comm(),
+                    true    // reuse
+                )
+            );
+
+            // Construct GAMGInterfaces
+            lduInterfacePtrsList newInterfaces(validInterface.size());
+            forAll(validInterface, intI)
+            {
+                if (validInterface[intI])
+                {
+                    word coupleType(fromProc);
+
+                    newInterfaces.set
+                    (
+                        intI,
+                        GAMGInterface::New
+                        (
+                            coupleType,
+                            intI,
+                            otherMeshes[proci-1].rawInterfaces(),
+                            fromProc
+                        ).ptr()
+                    );
+                }
+            }
+
+            otherMeshes[proci-1].addInterfaces
+            (
+                newInterfaces,
+                nonBlockingSchedule<processorGAMGInterface>
+                (
+                    newInterfaces
+                )
+            );
+        }
+    }
 }
 
 

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H

@@ -82,6 +82,17 @@ class lduPrimitiveMesh
         //- Get size of all meshes
         static label totalSize(const PtrList<lduPrimitiveMesh>&);
 
+        label findConnectedInterface
+        (
+            const lduMesh& myMesh,
+            const PtrList<lduPrimitiveMesh>& otherMeshes,
+            const labelPairList& procAndInterfaces,
+
+            const label nbrProci,
+            const label myRank
+        ) const;
+
+
         //- No copy construct
         lduPrimitiveMesh(const lduPrimitiveMesh&) = delete;
 
@@ -277,13 +288,12 @@ public:
                 const label meshI
             );
 
-            //- Gather meshes from other processors onto procIDs[0].
+            //- Gather meshes from other processors using agglomComm.
-            //  Received meshes get GAMGInterface and communicator comm
+            //  Received meshes get GAMGInterface.
             static void gather
             (
-                const label comm,
+                const label agglomComm,
                 const lduMesh& mesh,
-                const labelList& procIDs,
                 PtrList<lduPrimitiveMesh>& otherMeshes
             );
 

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMeshTemplates.C

@@ -84,6 +84,8 @@ Foam::lduSchedule Foam::lduPrimitiveMesh::nonBlockingSchedule
     // (where interfaces are defined) but must retain the full list length
     // for later (external) bookkeeping
 
+    schedule.setSize(patchEvali);
+
     return schedule;
 }
 

src/OpenFOAM/meshes/meshState/meshState.C

@@ -31,8 +31,12 @@ License
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-const Foam::word Foam::meshState::solverPerformanceDictName = "solver";
+const Foam::word Foam::meshState::controlsDictName = "controls";
 const Foam::word Foam::meshState::meshDictName = "mesh";
+const Foam::word Foam::meshState::solverPerformanceDictName = "solver";
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
 bool Foam::meshState::getBoolEntry
 (
@@ -87,6 +91,7 @@ Foam::meshState::meshState
     ),
     prevTimeIndex_(-1)
 {
+    (void)subDictOrAdd(controlsDictName);
     (void)subDictOrAdd(meshDictName);
     (void)subDictOrAdd(solverPerformanceDictName);
 }
@@ -100,6 +105,18 @@ void Foam::meshState::reset(const meshState& ms)
 }
 
 
+Foam::dictionary& Foam::meshState::controlsDict()
+{
+    return subDict(controlsDictName);
+}
+
+
+const Foam::dictionary& Foam::meshState::controlsDict() const
+{
+    return subDict(controlsDictName);
+}
+
+
 Foam::dictionary& Foam::meshState::meshDict()
 {
     return subDict(meshDictName);
@@ -126,25 +143,25 @@ const Foam::dictionary& Foam::meshState::solverPerformanceDict() const
 
 bool Foam::meshState::isFirstIteration() const
 {
-    return getBoolEntry(solverPerformanceDict(), "firstIteration");
+    return getBoolEntry(controlsDict(), "firstIteration");
 }
 
 
 bool Foam::meshState::isFinalIteration() const
 {
-    return getBoolEntry(solverPerformanceDict(), "finalIteration");
+    return getBoolEntry(controlsDict(), "finalIteration");
 }
 
 
 void Foam::meshState::setFirstIteration(bool on)
 {
-    return setBoolEntry(solverPerformanceDict(), "firstIteration", on);
+    return setBoolEntry(controlsDict(), "firstIteration", on);
 }
 
 
 void Foam::meshState::setFinalIteration(bool on)
 {
-    return setBoolEntry(solverPerformanceDict(), "finalIteration", on);
+    return setBoolEntry(controlsDict(), "finalIteration", on);
 }
 
 

src/OpenFOAM/meshes/meshState/meshState.H

@@ -31,7 +31,8 @@ Description
     Database for mesh data, solution data, solver performance and other
     reduced data.
 
-    polyMesh has a meshState data member so that every GeometricField has access.
+    polyMesh has a meshState data member so that
+    every GeometricField has access.
 
 SourceFiles
     meshState.C
@@ -88,9 +89,15 @@ public:
 
     // Static Data Members
 
-        static const word solverPerformanceDictName;
+        //- Name for "controls" sub-dictionary
+        static const word controlsDictName;
+
+        //- Name for "mesh" sub-dictionary
         static const word meshDictName;
 
+        //- Name for "solver" (solver performance) sub-dictionary
+        static const word solverPerformanceDictName;
+
 
     // Constructors
 
@@ -134,6 +141,10 @@ public:
         //- Reset the dictionary
         void reset(const meshState& ms);
 
+        //- Return the dictionary of controls
+        dictionary& controlsDict();
+        const dictionary& controlsDict() const;
+
         //- Return the dictionary of mesh data, typically populated by
         //- the polyMesh::checkXXX functions, checkMesh etc.
         dictionary& meshDict();
@@ -144,16 +155,16 @@ public:
         dictionary& solverPerformanceDict();
         const dictionary& solverPerformanceDict() const;
 
-        //- True if "firstIteration" entry exists and is set
+        //- True if "firstIteration" entry exists (in controls) and is set
         bool isFirstIteration() const;
 
-        //- True if "finalIteration" entry exists and is set
+        //- True if "finalIteration" entry exists (in controls) and is set
         bool isFinalIteration() const;
 
-        //- Add/remove "firstIteration" entry
+        //- Add/remove "firstIteration" entry (from controls)
         void setFirstIteration(bool on);
 
-        //- Add/remove "finalIteration" entry
+        //- Add/remove "finalIteration" entry (from controls)
         void setFinalIteration(bool on);
 
         //- Add/set the solverPerformance entry for the named field

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C

@@ -87,7 +87,12 @@ void Foam::globalMeshData::initProcAddr()
 
     if (Pstream::parRun())
     {
-        PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+        PstreamBuffers pBufs
+        (
+            Pstream::commsTypes::nonBlocking,
+            UPstream::msgType(),
+            mesh_.comm()
+        );
 
         // Send indices of my processor patches to my neighbours
         for (const label patchi : processorPatches_)
@@ -1738,7 +1743,7 @@ Foam::globalMeshData::globalMeshData(const polyMesh& mesh)
         processorTopology::New<processorPolyPatch>
         (
             mesh.boundaryMesh(),
-            UPstream::worldComm
+            mesh_.comm()
         )
     ),
     processorPatches_(),

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -146,6 +146,21 @@ void Foam::mapDistribute::printLayout(Ostream& os) const
 }
 
 
+Foam::UPtrList<const Foam::mapDistributeBase> Foam::mapDistribute::extractBase
+(
+    const UPtrList<const mapDistribute>& maps
+)
+{
+    UPtrList<const mapDistributeBase> baseMaps(maps.size());
+    forAll(maps, i)
+    {
+        const mapDistributeBase& map = maps[i];
+        baseMaps.set(i, &map);
+    }
+    return baseMaps;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::mapDistribute::mapDistribute() noexcept
@@ -437,6 +452,32 @@ Foam::mapDistribute::mapDistribute
 }
 
 
+Foam::mapDistribute::mapDistribute
+(
+    const UPtrList<const mapDistribute>& maps,
+    const labelList& localRanks,
+    const label newComm,
+    const labelListList& newToOldRanks, // from rank in newComm to
+                                        // ranks in (old)comm
+    labelList& startOfLocal,
+    List<Map<label>>& compactMaps
+)
+:
+    mapDistributeBase
+    (
+        extractBase(maps),
+        localRanks,
+        newComm,
+        newToOldRanks,
+        startOfLocal,
+        compactMaps
+    )
+{
+    // TBD. -have mapDistributeBase::add or something
+    //      -set transforms from individual maps
+}
+
+
 Foam::autoPtr<Foam::mapDistribute> Foam::mapDistribute::clone() const
 {
     return autoPtr<mapDistribute>::New(*this);

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H

@@ -198,6 +198,12 @@ class mapDistribute
             const TransformOp& top
         ) const;
 
+        //- Helper: convert mapDistribute to mapDistributeBase
+        static UPtrList<const mapDistributeBase> extractBase
+        (
+            const UPtrList<const mapDistribute>& maps
+        );
+
 
 public:
 
@@ -393,6 +399,27 @@ public:
             const label comm = UPstream::worldComm
         );
 
+        //- Construct from multiple maps and processor collation
+        //  Assumes all local data first. Sorts contributions of maps
+        //  in processor order i.e.  constructed map has all local data first.
+        //  Returns
+        //  - startOfLocal : per input map the start of the local data. Extends
+        //                   one beyond number of maps so overall local size
+        //                   is startOfLocal.last()
+        //  - compactMaps  : per input map from slot position in the input map
+        //                   to new slot position. (note there is no information
+        //                   returned about which processor it is from)
+        mapDistribute
+        (
+            const UPtrList<const mapDistribute>& maps,
+            const labelList& localRanks,
+            const label newComm,
+            const labelListList& newToOldRanks, // from rank in newComm to
+                                                // ranks in (old)comm
+            labelList& startOfLocal,            // per map start of local data
+            List<Map<label>>& compactMaps       // per map old slot to new slot
+        );
+
         //- Construct from Istream
         explicit mapDistribute(Istream& is);
 

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBase.C

@@ -399,9 +399,9 @@ void Foam::mapDistributeBase::calcCompactAddressing
 
     for (const label globalIdx : elements)
     {
-        if (globalIdx != -1 && !globalNumbering.isLocal(globalIdx))
+        if (globalIdx != -1 && !globalNumbering.isLocal(myRank, globalIdx))
         {
-            label proci = globalNumbering.whichProcID(globalIdx);
+            label proci = globalNumbering.whichProcID(myRank, globalIdx);
             nNonLocal[proci]++;
         }
     }
@@ -421,9 +421,9 @@ void Foam::mapDistributeBase::calcCompactAddressing
     // Collect all (non-local) elements needed.
     for (const label globalIdx : elements)
     {
-        if (globalIdx != -1 && !globalNumbering.isLocal(globalIdx))
+        if (globalIdx != -1 && !globalNumbering.isLocal(myRank, globalIdx))
         {
-            label proci = globalNumbering.whichProcID(globalIdx);
+            label proci = globalNumbering.whichProcID(myRank, globalIdx);
             label index = globalNumbering.toLocal(proci, globalIdx);
             label nCompact = compactMap[proci].size();
             compactMap[proci].insert(index, nCompact);
@@ -449,9 +449,9 @@ void Foam::mapDistributeBase::calcCompactAddressing
     {
         for (const label globalIdx : cCells)
         {
-            if (globalIdx != -1 && !globalNumbering.isLocal(globalIdx))
+            if (globalIdx != -1 && !globalNumbering.isLocal(myRank, globalIdx))
             {
-                label proci = globalNumbering.whichProcID(globalIdx);
+                label proci = globalNumbering.whichProcID(myRank, globalIdx);
                 nNonLocal[proci]++;
             }
         }
@@ -474,9 +474,9 @@ void Foam::mapDistributeBase::calcCompactAddressing
     {
         for (const label globalIdx : cCells)
         {
-            if (globalIdx != -1 && !globalNumbering.isLocal(globalIdx))
+            if (globalIdx != -1 && !globalNumbering.isLocal(myRank, globalIdx))
             {
-                label proci = globalNumbering.whichProcID(globalIdx);
+                label proci = globalNumbering.whichProcID(myRank, globalIdx);
                 label index = globalNumbering.toLocal(proci, globalIdx);
                 label nCompact = compactMap[proci].size();
                 compactMap[proci].insert(index, nCompact);
@@ -504,7 +504,7 @@ void Foam::mapDistributeBase::exchangeAddressing
 
     compactStart.setSize(nProcs);
     compactStart[myRank] = 0;
-    constructSize_ = globalNumbering.localSize();
+    constructSize_ = globalNumbering.localSize(myRank);
     forAll(compactStart, proci)
     {
         if (proci != myRank)
@@ -526,7 +526,7 @@ void Foam::mapDistributeBase::exchangeAddressing
         if (proci == myRank)
         {
             // All my own elements are used
-            label nLocal = globalNumbering.localSize();
+            label nLocal = globalNumbering.localSize(myRank);
             wantedRemoteElements[proci] = identity(nLocal);
             constructMap_[proci] = identity(nLocal);
         }
@@ -561,7 +561,7 @@ void Foam::mapDistributeBase::exchangeAddressing
     // Renumber elements
     for (label& elem : elements)
     {
-        elem = renumber(globalNumbering, compactMap, elem);
+        elem = renumber(globalNumbering, comm_, compactMap, elem);
     }
 }
 
@@ -584,7 +584,7 @@ void Foam::mapDistributeBase::exchangeAddressing
 
     compactStart.setSize(nProcs);
     compactStart[myRank] = 0;
-    constructSize_ = globalNumbering.localSize();
+    constructSize_ = globalNumbering.localSize(myRank);
     forAll(compactStart, proci)
     {
         if (proci != myRank)
@@ -606,7 +606,7 @@ void Foam::mapDistributeBase::exchangeAddressing
         if (proci == myRank)
         {
             // All my own elements are used
-            label nLocal = globalNumbering.localSize();
+            label nLocal = globalNumbering.localSize(myRank);
             wantedRemoteElements[proci] = identity(nLocal);
             constructMap_[proci] = identity(nLocal);
         }
@@ -643,7 +643,7 @@ void Foam::mapDistributeBase::exchangeAddressing
     {
         for (label& celli : cCells)
         {
-            celli = renumber(globalNumbering, compactMap, celli);
+            celli = renumber(globalNumbering, comm_, compactMap, celli);
         }
     }
 }
@@ -1005,6 +1005,452 @@ Foam::mapDistributeBase::mapDistributeBase
 {}
 
 
+Foam::mapDistributeBase::mapDistributeBase
+(
+    const UPtrList<const mapDistributeBase>& maps,
+    const labelList& localRanks,
+    const label newComm,
+    const labelListList& newToOldRanks,// from newComm to comm_
+    labelList& startOfLocal,
+    List<Map<label>>& compactMaps
+)
+:
+    constructSize_(0),
+    subHasFlip_(false),
+    constructHasFlip_(false),
+    comm_(-1),
+    schedulePtr_(nullptr)
+{
+    if (maps.empty())
+    {
+        return;
+    }
+
+    comm_ = newComm;
+    subHasFlip_ = maps[0].subHasFlip();
+    constructHasFlip_ = maps[0].constructHasFlip();
+
+    const label nNewRanks = newToOldRanks.size();
+    const label myNewRank = UPstream::myProcNo(newComm);
+    if (nNewRanks != UPstream::nProcs(newComm))
+    {
+        FatalErrorInFunction<< "nNewRanks:" << nNewRanks
+            << " nProcs:" << UPstream::nProcs(newComm)
+            << exit(FatalError);
+    }
+
+    if (localRanks.size() != maps.size())
+    {
+        FatalErrorInFunction
+            << "Number of maps:" << maps.size()
+            << " number of localRanks:" << localRanks.size()
+            << exit(FatalError);
+    }
+
+    // Sanity checks
+    const auto& map0 = maps[0];
+    forAll(maps, mapi)
+    {
+        const auto& map = maps[mapi];
+
+        if
+        (
+            (map.comm() != map0.comm())
+         || (map.subHasFlip() != map0.subHasFlip())
+         || (map.constructHasFlip() != map0.constructHasFlip())
+        )
+        {
+            FatalErrorInFunction
+                << "Maps should all be the same form"
+                << " Map " << mapi
+                << " has comm:" << map.comm()
+                << " subHasFlip:" << map.subHasFlip()
+                << " constructHasFlip:" << map.constructHasFlip()
+                << " which is different from map 0"
+                << exit(FatalError);
+        }
+
+        const label localRank = localRanks[mapi];
+        const auto& constructOwn = maps[mapi].constructMap()[localRank];
+        forAll(constructOwn, i)
+        {
+            if (constructOwn[i] != i)
+            {
+                FatalErrorInFunction
+                    << "Maps constructMap not identity."
+                    << " Map " << mapi
+                    << " constructMap:" << flatOutput(constructOwn)
+                    << exit(FatalError);
+            }
+        }
+    }
+
+
+    constructMap_.resize_nocopy(nNewRanks);
+    subMap_.resize_nocopy(nNewRanks);
+
+
+    // Store starts
+    startOfLocal.setSize(maps.size()+1);
+    compactMaps.resize_nocopy(maps.size());
+
+    label constructi = 0;
+    forAll(maps, mapi)
+    {
+        startOfLocal[mapi] = constructi;
+        const label localRank = localRanks[mapi];
+        const auto& map = maps[mapi].constructMap()[localRank];
+
+        // Presize compaction array
+        const label nRemote = maps[mapi].constructSize()-map.size();
+        compactMaps[mapi].resize(2*nRemote);
+
+        constructi += map.size();
+    }
+    startOfLocal.last() = constructi;
+
+
+    // Determine start of constructed remote data. This is used to get the
+    // local offset which can then be used to get the relative subMap location.
+    labelListList startOfRemote(maps.size());
+    forAll(maps, mapi)
+    {
+        const label nOldProcs = maps[mapi].constructMap().size();
+        labelList& starts = startOfRemote[mapi];
+
+        starts.setSize(nOldProcs, labelMax);
+        forAll(maps[mapi].constructMap(), oldProci)
+        {
+            const labelList& map = maps[mapi].constructMap()[oldProci];
+            forAll(map, i)
+            {
+                const label index
+                (
+                    constructHasFlip_
+                  ? mag(map[i])-1
+                  : map[i]
+                );
+                starts[oldProci] = min(starts[oldProci], index);
+            }
+        }
+    }
+
+
+    // Construct map
+    // ~~~~~~~~~~~~~
+    // - all localRanks:
+    //      - data gets appended in map order
+    //      - map is just an offset (startOfLocal)
+    // - all previously remote ranks:
+    //      - data is already present according to startOfLocal
+    //      - map is old-to-new index
+    // - all still remote ranks:
+    //      - data gets appended in map order after the startOfLocal
+    //      - map is old-to-new index
+
+
+    // Append local (= myRank) data. TBD: assumes subMap and constructMap
+    // are identity maps.
+    {
+        labelList& myConstruct = constructMap_[myNewRank];
+        myConstruct.resize_nocopy(constructi);
+        constructi = 0;
+        forAll(maps, mapi)
+        {
+            const label localRank = localRanks[mapi];
+            const auto& map = maps[mapi].constructMap()[localRank];
+            const label offset = startOfLocal[mapi];
+
+            forAll(map, i)
+            {
+                if (constructHasFlip_)
+                {
+                    forAll(map, i)
+                    {
+                        if (map[i] < 0)
+                        {
+                            myConstruct[constructi++] = map[i]-offset;
+                        }
+                        else
+                        {
+                            myConstruct[constructi++] = map[i]+offset;
+                        }
+                    }
+                }
+                else
+                {
+                    myConstruct[constructi++] = map[i]+offset;
+                }
+            }
+        }
+    }
+
+    // Filter remote construct data
+    {
+        // Remote ranks that are now local
+        //  - store new index for mapping stencils
+        //  - no need to construct since already
+        const auto& oldProcs = newToOldRanks[myNewRank];
+
+        forAll(maps, mapi)
+        {
+            for (const label oldProci : oldProcs)
+            {
+                if (oldProci != localRanks[mapi])
+                {
+                    const auto& map = maps[mapi].constructMap()[oldProci];
+
+                    if (!map.size())
+                    {
+                        continue;
+                    }
+
+
+                    // The slots come from a local map so we can look up the
+                    // new location
+                    const label sourceMapi = localRanks.find(oldProci);
+                    const auto& subMap =
+                        maps[sourceMapi].subMap()[localRanks[mapi]];
+
+                    //Pout<< "From oldRank:" << oldProci
+                    //    << " sending to masterRank:" << localRanks[mapi]
+                    //    << " elements:" << flatOutput(subMap)
+                    //    << nl
+                    //    << "   received as elements:" << flatOutput(map)
+                    //    << endl;
+
+                    if (map.size() != subMap.size())
+                    {
+                        FatalErrorInFunction << "Problem:"
+                            << "oldProci:" << oldProci
+                            << " mapi:" << mapi
+                            << " constructMap:" << map.size()
+                            << " sourceMapi:" << sourceMapi
+                            << " subMap:" << subMap.size()
+                            << exit(FatalError);
+                    }
+
+                    const label offset = startOfLocal[sourceMapi];
+                    // Construct map starts after the local data
+                    const label nMapLocal = startOfRemote[mapi][oldProci];
+
+                    auto& cptMap = compactMaps[mapi];
+                    forAll(map, i)
+                    {
+                        // old slot position to new slot position
+                        const label index
+                        (
+                            constructHasFlip_
+                          ? mag(map[i])-1
+                          : map[i]
+                        );
+                        const label newIndex = subMap[index-nMapLocal]+offset;
+
+                        // Note: should always warn for duplicates? Or only if
+                        // different?
+                        if
+                        (
+                           !cptMap.insert(index, newIndex)
+                         && cptMap[index] != newIndex
+                        )
+                        {
+                            FatalErrorInFunction<< "Duplicate insertion"
+                                << "From oldProc:" << oldProci
+                                << " on map:" << mapi
+                                << " at index:" << i
+                                << " have construct slot:" << index
+                                << " new index:" << newIndex
+                                << " but already have entry:" << cptMap[index]
+                                << " on for that slot"
+                                << exit(FatalError);
+                        }
+                    }
+                }
+            }
+        }
+
+
+        // Remote ranks that are still remote
+        //  - store new index for mapping stencils
+        //  - append to construction
+
+        // Either loop over all old ranks and filter out ones already handled
+        // or loop over all new ranks and avoid myNewRank
+
+        forAll(newToOldRanks, newProci)
+        {
+            if (newProci != myNewRank)
+            {
+                const auto& oldProcs = newToOldRanks[newProci];
+
+                label allSize = 0;
+                forAll(maps, mapi)
+                {
+                    for (const label oldProci : oldProcs)
+                    {
+                        allSize += maps[mapi].constructMap()[oldProci].size();
+                    }
+                }
+
+                labelList& myConstruct = constructMap_[newProci];
+                myConstruct.resize_nocopy(allSize);
+
+                allSize = 0;
+                forAll(maps, mapi)
+                {
+                    for (const label oldProci : oldProcs)
+                    {
+                        const auto& map = maps[mapi].constructMap()[oldProci];
+                        // Construct map starts after the local data
+                        const label nMapLocal = startOfRemote[mapi][oldProci];
+                        SubList<label> slice(myConstruct, map.size(), allSize);
+
+                        if (constructHasFlip_)
+                        {
+                            forAll(map, i)
+                            {
+                                if (map[i] < 0)
+                                {
+                                    slice[i] = map[i]+nMapLocal-constructi;
+                                }
+                                else
+                                {
+                                    slice[i] = map[i]-nMapLocal+constructi;
+                                }
+                            }
+
+                            auto& cptMap = compactMaps[mapi];
+                            forAll(map, i)
+                            {
+                                cptMap.insert(mag(map[i])-1,mag(slice[i])-1);
+                            }
+                        }
+                        else
+                        {
+                            forAll(map, i)
+                            {
+                                slice[i] = map[i]-nMapLocal+constructi;
+                                compactMaps[mapi].insert(map[i], slice[i]);
+                            }
+                        }
+                        allSize += map.size();
+                        constructi += map.size();
+                    }
+                }
+            }
+        }
+    }
+
+
+    // Sub (=send) map
+    // ~~~~~~~~~~~~~~~
+    // - all localRanks:
+    //      - get appended in map order
+    // - all previously remote ranks:
+    //      - not needed. Stay empty
+    // - all still remote ranks:
+    //      - convert to new local index
+
+    // Append local (= myRank) data
+    {
+        label allSize = 0;
+        forAll(maps, mapi)
+        {
+            const label localRank = localRanks[mapi];
+            allSize += maps[mapi].subMap()[localRank].size();
+        }
+
+        labelList& mySub = subMap_[myNewRank];
+        mySub.resize_nocopy(allSize);
+        allSize = 0;
+        forAll(maps, mapi)
+        {
+            const label localRank = localRanks[mapi];
+            const auto& map = maps[mapi].subMap()[localRank];
+            SubList<label> slice(mySub, map.size(), allSize);
+
+            if (subHasFlip_)
+            {
+                forAll(slice, i)
+                {
+                    if (map[i] < 0)
+                    {
+                        slice[i] = map[i]-startOfLocal[mapi];
+                    }
+                    else
+                    {
+                        slice[i] = map[i]+startOfLocal[mapi];
+                    }
+                }
+            }
+            else
+            {
+                forAll(slice, i)
+                {
+                    slice[i] = map[i]+startOfLocal[mapi];
+                }
+            }
+            allSize += map.size();
+        }
+    }
+    // Filter remote sub data
+    forAll(newToOldRanks, newProci)
+    {
+        if (newProci != myNewRank)
+        {
+            const auto& oldProcs = newToOldRanks[newProci];
+
+            label allSize = 0;
+            forAll(maps, mapi)
+            {
+                for (const label oldProci : oldProcs)
+                {
+                    allSize += maps[mapi].subMap()[oldProci].size();
+                }
+            }
+
+            labelList& mySub = subMap_[newProci];
+            mySub.resize_nocopy(allSize);
+
+            allSize = 0;
+            for (const label oldProci : oldProcs)
+            {
+                forAll(maps, mapi)
+                {
+                    const auto& map = maps[mapi].subMap()[oldProci];
+                    SubList<label> slice(mySub, map.size(), allSize);
+                    if (subHasFlip_)
+                    {
+                        forAll(map, i)
+                        {
+                            if (map[i] < 0)
+                            {
+                                slice[i] = map[i]-startOfLocal[mapi];
+                            }
+                            else
+                            {
+                                slice[i] = map[i]+startOfLocal[mapi];
+                            }
+                        }
+                    }
+                    else
+                    {
+                        forAll(map, i)
+                        {
+                            slice[i] = map[i]+startOfLocal[mapi];
+                        }
+                    }
+                    allSize += map.size();
+                }
+            }
+        }
+    }
+
+
+    constructSize_ = constructi;
+}
+
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::labelList Foam::mapDistributeBase::subMapSizes() const
@@ -1088,21 +1534,24 @@ void Foam::mapDistributeBase::transfer(mapDistributeBase& rhs)
 Foam::label Foam::mapDistributeBase::renumber
 (
     const globalIndex& globalNumbering,
+    const label comm,
     const List<Map<label>>& compactMap,
     const label globalI
 )
 {
+    const label myRank = Pstream::myProcNo(comm);
+
     if (globalI == -1)
     {
         return globalI;
     }
-    if (globalNumbering.isLocal(globalI))
+    if (globalNumbering.isLocal(myRank, globalI))
     {
-        return globalNumbering.toLocal(globalI);
+        return globalNumbering.toLocal(myRank, globalI);
     }
     else
     {
-        label proci = globalNumbering.whichProcID(globalI);
+        label proci = globalNumbering.whichProcID(myRank, globalI);
         label index = globalNumbering.toLocal(proci, globalI);
         return compactMap[proci][index];
     }

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBase.H

@@ -484,6 +484,27 @@ public:
             const label comm = UPstream::worldComm
         );
 
+        //- Construct from multiple maps and processor collation
+        //  Assumes all local data first. Sorts contributions of maps
+        //  in processor order i.e.  constructed map has all local data first.
+        //  Returns
+        //  - startOfLocal : per input map the start of the local data. Extends
+        //                   one beyond number of maps so overall local size
+        //                   is startOfLocal.last()
+        //  - compactMaps  : per input map from slot position in the input map
+        //                   to new slot position. (note there is no information
+        //                   returned about which processor it is from)
+        mapDistributeBase
+        (
+            const UPtrList<const mapDistributeBase>& maps,
+            const labelList& localRanks,
+            const label newComm,
+            const labelListList& newToOldRanks, // from rank in newComm to
+                                                // ranks in (old)comm
+            labelList& startOfLocal,            // per map start of local data
+            List<Map<label>>& compactMaps       // per map old slot to new slot
+        );
+
         //- Construct from Istream
         explicit mapDistributeBase(Istream& is);
 
@@ -645,6 +666,7 @@ public:
         static label renumber
         (
             const globalIndex&,
+            const label comm,
             const List<Map<label>>& compactMap,
             const label globalElement
         );
@@ -665,6 +687,36 @@ public:
             const bool hasFlip
         );
 
+        //- Helper for renumbering the map elements. Assumes local elements
+        //- are first, followed by any remote. Local elements get offset,
+        //- remote elements are mapped.
+        //
+        //  \param[in,out] map  The map to be renumbered
+        //  \param localSize    elements < localSize are offset
+        //  \param offset       offset
+        //  \param cMap         map for non-local elements
+        //  \param hasFlip      True if map has flip addressing
+        //
+        //  \return max-size needed for new addressing (eg, constructSize)
+        static label renumberMap
+        (
+            labelList& map,
+            const label localSize,
+            const label offset,
+            const Map<label>& cMap,
+            const bool hasFlip
+        );
+
+        //- Helper for a list of maps. Calls above renumberMap for all elements
+        static label renumberMap
+        (
+            labelListList& mapElements,
+            const label localSize,
+            const label offset,
+            const Map<label>& cMap,
+            const bool hasFlip
+        );
+
 
     // Compaction
 

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBaseSubset.C

@@ -329,6 +329,104 @@ Foam::label Foam::mapDistributeBase::renumberMap
 }
 
 
+Foam::label Foam::mapDistributeBase::renumberMap
+(
+    labelList& map,
+    const label localSize,
+    const label offset,
+    const Map<label>& cMap,
+    const bool hasFlip
+)
+{
+    label maxIndex = -1;
+
+    // Transcribe the map
+    if (hasFlip)
+    {
+        for (label& val : map)
+        {
+            // Unflip indexed value
+            const label index = mag(val)-1;
+            if (index < localSize)
+            {
+                // Local element
+                if (val < 0)
+                {
+                    val -= offset;
+                }
+                else
+                {
+                    val += offset;
+                }
+            }
+            else
+            {
+                // Remote element
+                if (val < 0)
+                {
+                    val = -cMap[index]-1;
+                }
+                else
+                {
+                    val = cMap[index]+1;
+                }
+            }
+            maxIndex = max(maxIndex, mag(val)-1);
+        }
+    }
+    else
+    {
+        for (label& val : map)
+        {
+            // Get indexed value (no flipping)
+            if (val < localSize)
+            {
+                val += offset;
+            }
+            else
+            {
+                val = cMap[val];
+            }
+            maxIndex = max(maxIndex, val);
+        }
+    }
+
+    return (maxIndex+1);
+}
+
+
+Foam::label Foam::mapDistributeBase::renumberMap
+(
+    labelListList& mapElements,
+    const label localSize,
+    const label offset,
+    const Map<label>& cMap,
+    const bool hasFlip
+)
+{
+    label maxIndex = -1;
+
+    // Transcribe the map
+    for (labelList& map : mapElements)
+    {
+        maxIndex = max
+        (
+            maxIndex,
+            renumberMap
+            (
+                map,
+                localSize,
+                offset,
+                cMap,
+                hasFlip
+            )
+        );
+    }
+
+    return (maxIndex+1);
+}
+
+
 void Foam::mapDistributeBase::renumberVisitOrder
 (
     const labelUList& origElements,

src/Pstream/Allwmake

@@ -3,9 +3,10 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/wmake/scripts/AllwmakeParseArguments
 
 #------------------------------------------------------------------------------
+# Never need/want openmp for MPI interfaces, or for static linkage
 
-echo "wmake $targetType dummy (mpi=$WM_MPLIB)"
+echo "wmake -no-openmp $targetType${targetType:+ }dummy (mpi=$WM_MPLIB)"
-wmake $targetType dummy
+wmake -no-openmp $targetType dummy
 
 ./Allwmake-mpi $targetType $*
 

src/Pstream/Allwmake-mpi

@@ -7,6 +7,7 @@ cd "${0%/*}" || exit                                # Run from this directory
 # Environment
 #   - FOAM_MPI_LIBBIN (optional: defaults to FOAM_LIBBIN/FOAM_MPI)
 
-wmakeLibMpi mpi
+# Never need/want openmp for MPI interfaces
+wmakeLibMpi -no-openmp mpi
 
 #------------------------------------------------------------------------------

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2019 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -189,105 +189,150 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate
 {
     const label patchi = patch.index();
 
-    const tmp<scalarField> tnutw = turbModel.nut(patchi);
+    const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
-    const scalarField& nutw = tnutw();
+    const scalar Cmu75 = pow(wallCoeffs_.Cmu(), 0.75);
+    const scalar kappa = wallCoeffs_.kappa();
+    const scalar yPlusLam = wallCoeffs_.yPlusLam();
 
     const scalarField& y = turbModel.y()[patchi];
 
+    const labelUList& faceCells = patch.faceCells();
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
     const tmp<volScalarField> tk = turbModel.k();
     const volScalarField& k = tk();
 
-    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
+    // Calculate y-plus
-
+    const auto yPlus = [&](const label facei) -> scalar
-    const scalarField magGradUw(mag(Uw.snGrad()));
-
-    const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
-    const scalar Cmu75 = pow(wallCoeffs_.Cmu(), 0.75);
-    const scalar kappa = wallCoeffs_.kappa();
-    const scalar yPlusLam = wallCoeffs_.yPlusLam();
-
-    // Set epsilon and G
-    forAll(nutw, facei)
     {
-        const label celli = patch.faceCells()[facei];
+        return
+        (
+            Cmu25*y[facei]*sqrt(k[faceCells[facei]])/nuw[facei]
+        );
+    };
 
-        const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei];
+    // Contribution from the viscous sublayer
+    const auto epsilonVis = [&](const label facei) -> scalar
+    {
+        return
+        (
+            cornerWeights[facei]*2.0*k[faceCells[facei]]*nuw[facei]
+          / sqr(y[facei])
+        );
+    };
 
-        const scalar w = cornerWeights[facei];
+    // Contribution from the inertial sublayer
+    const auto epsilonLog = [&](const label facei) -> scalar
+    {
+        return
+        (
+            cornerWeights[facei]*Cmu75*pow(k[faceCells[facei]], 1.5)
+          / (kappa*y[facei])
+        );
+    };
 
-        // Contribution from the viscous sublayer
+    switch (blender_)
-        const scalar epsilonVis = w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
+    {
-
+        case blenderType::STEPWISE:
-        // Contribution from the inertial sublayer
-        const scalar epsilonLog = w*Cmu75*pow(k[celli], 1.5)/(kappa*y[facei]);
-
-        switch (blender_)
         {
-            case blenderType::STEPWISE:
+            forAll(faceCells, facei)
             {
-                if (lowReCorrection_ && yPlus < yPlusLam)
+                if (lowReCorrection_ && yPlus(facei) < yPlusLam)
                 {
-                    epsilon0[celli] += epsilonVis;
+                    epsilon0[faceCells[facei]] += epsilonVis(facei);
                 }
                 else
                 {
-                    epsilon0[celli] += epsilonLog;
+                    epsilon0[faceCells[facei]] += epsilonLog(facei);
                 }
-                break;
-            }
-
-            case blenderType::BINOMIAL:
-            {
-                // (ME:Eqs. 15-16)
-                epsilon0[celli] +=
-                    pow
-                    (
-                        pow(epsilonVis, n_) + pow(epsilonLog, n_),
-                        scalar(1)/n_
-                    );
-                break;
-            }
-
-            case blenderType::MAX:
-            {
-                // (PH:Eq. 27)
-                epsilon0[celli] += max(epsilonVis, epsilonLog);
-                break;
-            }
-
-            case blenderType::EXPONENTIAL:
-            {
-                // (PH:p. 193)
-                const scalar Gamma = 0.001*pow4(yPlus)/(scalar(1) + yPlus);
-                const scalar invGamma = scalar(1)/(Gamma + ROOTVSMALL);
-                epsilon0[celli] +=
-                    epsilonVis*exp(-Gamma) + epsilonLog*exp(-invGamma);
-                break;
-            }
-
-            case blenderType::TANH:
-            {
-                // (KAS:Eqs. 33-34)
-                const scalar phiTanh = tanh(pow4(0.1*yPlus));
-                const scalar b1 = epsilonVis + epsilonLog;
-                const scalar b2 =
-                    pow(pow(epsilonVis, 1.2) + pow(epsilonLog, 1.2), 1.0/1.2);
-
-                epsilon0[celli] += phiTanh*b1 + (1 - phiTanh)*b2;
-                break;
             }
+            break;
         }
 
-        if (!lowReCorrection_ || (yPlus > yPlusLam))
+        case blenderType::BINOMIAL:
         {
-            G0[celli] +=
+            forAll(faceCells, facei)
-                w
+            {
+                // (ME:Eqs. 15-16)
+                epsilon0[faceCells[facei]] +=
+                    pow
+                    (
+                        pow(epsilonVis(facei), n_) + pow(epsilonLog(facei), n_),
+                        scalar(1)/n_
+                    );
+            }
+            break;
+        }
+
+        case blenderType::MAX:
+        {
+            forAll(faceCells, facei)
+            {
+                // (PH:Eq. 27)
+                epsilon0[faceCells[facei]] +=
+                    max(epsilonVis(facei), epsilonLog(facei));
+            }
+            break;
+        }
+
+        case blenderType::EXPONENTIAL:
+        {
+            forAll(faceCells, facei)
+            {
+                // (PH:p. 193)
+                const scalar yPlusFace = yPlus(facei);
+                const scalar Gamma =
+                    0.001*pow4(yPlusFace)/(scalar(1) + yPlusFace);
+                const scalar invGamma = scalar(1)/(Gamma + ROOTVSMALL);
+
+                epsilon0[faceCells[facei]] +=
+                (
+                    epsilonVis(facei)*exp(-Gamma)
+                  + epsilonLog(facei)*exp(-invGamma)
+                );
+            }
+            break;
+        }
+
+        case blenderType::TANH:
+        {
+            forAll(faceCells, facei)
+            {
+                // (KAS:Eqs. 33-34)
+                const scalar epsilonVisFace = epsilonVis(facei);
+                const scalar epsilonLogFace = epsilonLog(facei);
+                const scalar b1 = epsilonVisFace + epsilonLogFace;
+                const scalar b2 =
+                    pow
+                    (
+                        pow(epsilonVisFace, 1.2) + pow(epsilonLogFace, 1.2),
+                        1.0/1.2
+                    );
+                const scalar phiTanh = tanh(pow4(0.1*yPlus(facei)));
+
+                epsilon0[faceCells[facei]] += phiTanh*b1 + (1 - phiTanh)*b2;
+            }
+            break;
+        }
+    }
+
+    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
+    const scalarField magGradUw(mag(Uw.snGrad()));
+
+    const tmp<scalarField> tnutw = turbModel.nut(patchi);
+    const scalarField& nutw = tnutw();
+
+    forAll(faceCells, facei)
+    {
+        if (!lowReCorrection_ || (yPlus(facei) > yPlusLam))
+        {
+            G0[faceCells[facei]] +=
+                cornerWeights[facei]
                *(nutw[facei] + nuw[facei])
                *magGradUw[facei]
-               *Cmu25*sqrt(k[celli])
+               *Cmu25*sqrt(k[faceCells[facei]])
                /(kappa*y[facei]);
         }
     }

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,7 +32,6 @@ License
 #include "volFields.H"
 #include "addToRunTimeSelectionTable.H"
 
-
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField>
@@ -40,6 +39,10 @@ Foam::nutUWallFunctionFvPatchScalarField::calcNut() const
 {
     const label patchi = patch().index();
 
+    const scalar kappa = wallCoeffs_.kappa();
+    const scalar E = wallCoeffs_.E();
+    const scalar yPlusLam = wallCoeffs_.yPlusLam();
+
     const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
@@ -52,88 +55,107 @@ Foam::nutUWallFunctionFvPatchScalarField::calcNut() const
     const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
     const scalarField magUp(mag(Uw.patchInternalField() - Uw));
 
-    const tmp<scalarField> tnuw = turbModel.nu(patchi);
-    const scalarField& nuw = tnuw();
-
-    const scalar kappa = wallCoeffs_.kappa();
-    const scalar E = wallCoeffs_.E();
-    const scalar yPlusLam = wallCoeffs_.yPlusLam();
-
     tmp<scalarField> tyPlus = calcYPlus(magUp);
     const scalarField& yPlus = tyPlus();
 
+    // Viscous sublayer contribution
+    const tmp<scalarField> tnutVis = turbModel.nu(patchi);
+    const scalarField& nutVis = tnutVis();
+
+    // Inertial sublayer contribution
+    const auto nutLog = [&](const label facei) -> scalar
+    {
+        const scalar yPlusFace = yPlus[facei];
+        return
+        (
+            nutVis[facei]*yPlusFace*kappa/log(max(E*yPlusFace, 1 + 1e-4))
+        );
+    };
+
     auto tnutw = tmp<scalarField>::New(patch().size(), Zero);
     auto& nutw = tnutw.ref();
 
-    forAll(yPlus, facei)
+    switch (blender_)
     {
-        // Viscous sublayer contribution
+        case blenderType::STEPWISE:
-        const scalar nutVis = 0;
-
-        // Inertial sublayer contribution
-        const scalar nutLog =
-            nuw[facei]
-           *(yPlus[facei]*kappa/log(max(E*yPlus[facei], 1 + 1e-4)) - 1.0);
-
-        switch (blender_)
         {
-            case blenderType::STEPWISE:
+            forAll(nutw, facei)
             {
                 if (yPlus[facei] > yPlusLam)
                 {
-                    nutw[facei] = nutLog;
+                    nutw[facei] = nutLog(facei);
                 }
                 else
                 {
-                    nutw[facei] = nutVis;
+                    nutw[facei] = nutVis[facei];
                 }
-                break;
             }
+            break;
+        }
 
-            case blenderType::MAX:
+        case blenderType::MAX:
+        {
+            forAll(nutw, facei)
             {
                 // (PH:Eq. 27)
-                nutw[facei] = max(nutVis, nutLog);
+                nutw[facei] = max(nutVis[facei], nutLog(facei));
-                break;
             }
+            break;
+        }
 
-            case blenderType::BINOMIAL:
+        case blenderType::BINOMIAL:
+        {
+            forAll(nutw, facei)
             {
                 // (ME:Eqs. 15-16)
                 nutw[facei] =
                     pow
                     (
-                        pow(nutVis, n_) + pow(nutLog, n_),
+                        pow(nutVis[facei], n_) + pow(nutLog(facei), n_),
                         scalar(1)/n_
                     );
-                break;
             }
+            break;
+        }
 
-            case blenderType::EXPONENTIAL:
+        case blenderType::EXPONENTIAL:
+        {
+            forAll(nutw, facei)
             {
                 // (PH:Eq. 31)
                 const scalar Gamma =
                     0.01*pow4(yPlus[facei])/(1 + 5*yPlus[facei]);
                 const scalar invGamma = scalar(1)/(Gamma + ROOTVSMALL);
 
-                nutw[facei] = nutVis*exp(-Gamma) + nutLog*exp(-invGamma);
+                nutw[facei] =
-                break;
+                    nutVis[facei]*exp(-Gamma) + nutLog(facei)*exp(-invGamma);
             }
+            break;
+        }
 
-            case blenderType::TANH:
+        case blenderType::TANH:
+        {
+            forAll(nutw, facei)
             {
                 // (KAS:Eqs. 33-34)
-                const scalar phiTanh = tanh(pow4(0.1*yPlus[facei]));
+                const scalar nutLogFace = nutLog(facei);
-                const scalar b1 = nutVis + nutLog;
+                const scalar b1 = nutVis[facei] + nutLogFace;
                 const scalar b2 =
-                    pow(pow(nutVis, 1.2) + pow(nutLog, 1.2), 1.0/1.2);
+                    pow
+                    (
+                        pow(nutVis[facei], 1.2) + pow(nutLogFace, 1.2),
+                        1.0/1.2
+                    );
+                const scalar phiTanh = tanh(pow4(0.1*yPlus[facei]));
 
                 nutw[facei] = phiTanh*b1 + (1 - phiTanh)*b2;
-                break;
             }
+            break;
         }
     }
 
+    nutw -= nutVis;
+
     return tnutw;
 }
 

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016, 2019 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,6 +40,11 @@ calcNut() const
 {
     const label patchi = patch().index();
 
+    const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
+    const scalar kappa = wallCoeffs_.kappa();
+    const scalar E = wallCoeffs_.E();
+    const scalar yPlusLam = wallCoeffs_.yPlusLam();
+
     const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
@@ -49,93 +54,118 @@ calcNut() const
         )
     );
 
+    const labelUList& faceCells = patch().faceCells();
+
     const scalarField& y = turbModel.y()[patchi];
 
     const tmp<volScalarField> tk = turbModel.k();
     const volScalarField& k = tk();
 
-    const tmp<scalarField> tnuw = turbModel.nu(patchi);
+    // Viscous sublayer contribution
-    const scalarField& nuw = tnuw();
+    const tmp<scalarField> tnutVis = turbModel.nu(patchi);
+    const scalarField& nutVis = tnutVis();
 
-    const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
+    // Calculate y-plus
-    const scalar kappa = wallCoeffs_.kappa();
+    const auto yPlus = [&](const label facei) -> scalar
-    const scalar E = wallCoeffs_.E();
+    {
-    const scalar yPlusLam = wallCoeffs_.yPlusLam();
+        return (Cmu25*y[facei]*sqrt(k[faceCells[facei]])/nutVis[facei]);
+    };
+
+    // Inertial sublayer contribution
+    const auto nutLog = [&](const label facei) -> scalar
+    {
+        const scalar yPlusFace = yPlus(facei);
+        return
+        (
+            nutVis[facei]*yPlusFace*kappa
+          / log(max(E*yPlusFace, 1 + 1e-4))
+        );
+    };
 
     auto tnutw = tmp<scalarField>::New(patch().size(), Zero);
     auto& nutw = tnutw.ref();
 
-    forAll(nutw, facei)
+    switch (blender_)
     {
-        const label celli = patch().faceCells()[facei];
+        case blenderType::STEPWISE:
-
-        const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei];
-
-        // Viscous sublayer contribution
-        const scalar nutVis = 0;
-
-        // Inertial sublayer contribution
-        const scalar nutLog =
-            nuw[facei]*(yPlus*kappa/log(max(E*yPlus, 1 + 1e-4)) - 1.0);
-
-        switch (blender_)
         {
-            case blenderType::STEPWISE:
+            forAll(nutw, facei)
             {
-                if (yPlus > yPlusLam)
+                if (yPlus(facei) > yPlusLam)
                 {
-                    nutw[facei] = nutLog;
+                    nutw[facei] = nutLog(facei);
                 }
                 else
                 {
-                    nutw[facei] = nutVis;
+                    nutw[facei] = nutVis[facei];
                 }
-                break;
             }
+            break;
+        }
 
-            case blenderType::MAX:
+        case blenderType::MAX:
+        {
+            forAll(nutw, facei)
             {
                 // (PH:Eq. 27)
-                nutw[facei] = max(nutVis, nutLog);
+                nutw[facei] = max(nutVis[facei], nutLog(facei));
-                break;
             }
+            break;
+        }
 
-            case blenderType::BINOMIAL:
+        case blenderType::BINOMIAL:
+        {
+            forAll(nutw, facei)
             {
                 // (ME:Eqs. 15-16)
                 nutw[facei] =
                     pow
                     (
-                        pow(nutVis, n_) + pow(nutLog, n_),
+                        pow(nutVis[facei], n_) + pow(nutLog(facei), n_),
                         scalar(1)/n_
                     );
-                break;
             }
+            break;
+        }
 
-            case blenderType::EXPONENTIAL:
+        case blenderType::EXPONENTIAL:
+        {
+            forAll(nutw, facei)
             {
                 // (PH:Eq. 31)
-                const scalar Gamma = 0.01*pow4(yPlus)/(1 + 5*yPlus);
+                const scalar yPlusFace = yPlus(facei);
+                const scalar Gamma = 0.01*pow4(yPlusFace)/(1 + 5*yPlusFace);
                 const scalar invGamma = scalar(1)/(Gamma + ROOTVSMALL);
 
-                nutw[facei] = nutVis*exp(-Gamma) + nutLog*exp(-invGamma);
+                nutw[facei] =
-                break;
+                    nutVis[facei]*exp(-Gamma) + nutLog(facei)*exp(-invGamma);
             }
+            break;
+        }
 
-            case blenderType::TANH:
+        case blenderType::TANH:
+        {
+            forAll(nutw, facei)
             {
                 // (KAS:Eqs. 33-34)
-                const scalar phiTanh = tanh(pow4(0.1*yPlus));
+                const scalar nutLogFace = nutLog(facei);
-                const scalar b1 = nutVis + nutLog;
+                const scalar b1 = nutVis[facei] + nutLogFace;
                 const scalar b2 =
-                    pow(pow(nutVis, 1.2) + pow(nutLog, 1.2), 1.0/1.2);
+                    pow
+                    (
+                        pow(nutVis[facei], 1.2) + pow(nutLogFace, 1.2),
+                        1.0/1.2
+                    );
+                const scalar phiTanh = tanh(pow4(0.1*yPlus(facei)));
 
                 nutw[facei] = phiTanh*b1 + (1 - phiTanh)*b2;
-                break;
             }
+            break;
         }
     }
 
+    nutw -= nutVis;
+
     return tnutw;
 }
 
@@ -231,8 +261,8 @@ yPlus() const
     tmp<scalarField> tkwc = k.boundaryField()[patchi].patchInternalField();
     const scalarField& kwc = tkwc();
 
-    tmp<scalarField> tnuw = turbModel.nu(patchi);
+    tmp<scalarField> tnutVis = turbModel.nu(patchi);
-    const scalarField& nuw = tnuw();
+    const scalarField& nutVis = tnutVis();
 
     tmp<scalarField> tnuEff = turbModel.nuEff(patchi);
     const scalarField& nuEff = tnuEff();
@@ -249,13 +279,13 @@ yPlus() const
     forAll(yPlus, facei)
     {
         // inertial sublayer
-        yPlus[facei] = Cmu25*y[facei]*sqrt(kwc[facei])/nuw[facei];
+        yPlus[facei] = Cmu25*y[facei]*sqrt(kwc[facei])/nutVis[facei];
 
         if (yPlusLam > yPlus[facei])
         {
             // viscous sublayer
             yPlus[facei] =
-                y[facei]*sqrt(nuEff[facei]*magGradUw[facei])/nuw[facei];
+                y[facei]*sqrt(nuEff[facei]*magGradUw[facei])/nutVis[facei];
         }
     }