ENH: functionObjects: add optional flag to suppress warnings for missing patches

ENH: solver function objects: add outer-loop convergence checks

See merge request Development/openfoam!699

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/setRegionFluidFields.H

@@ -8,11 +8,13 @@
 
     volVectorField& U1 = phase1.URef();
     surfaceScalarField& phi1 = phase1.phiRef();
-    const surfaceScalarField& alphaPhi1 = phase1.alphaPhi();
+    const tmp<surfaceScalarField> talphaPhi1 = phase1.alphaPhi();
+    const auto& alphaPhi1 = talphaPhi1();
 
     volVectorField& U2 = phase2.URef();
     surfaceScalarField& phi2 = phase2.phiRef();
-    const surfaceScalarField& alphaPhi2 = phase2.alphaPhi();
+    const tmp<surfaceScalarField> talphaPhi2 = phase2.alphaPhi();
+    const auto& alphaPhi2 = talphaPhi2();
 
     surfaceScalarField& phi = fluid.phi();
 

applications/solvers/multiphase/compressibleInterFoam/VoFphaseCompressibleTurbulenceModels/compressibleInterPhaseTransportModel.C

@@ -74,8 +74,11 @@ Foam::compressibleInterPhaseTransportModel::compressibleInterPhaseTransportModel
         const volScalarField& alpha1(mixture_.alpha1());
         const volScalarField& alpha2(mixture_.alpha2());
 
-        const volScalarField& rho1 = mixture_.thermo1().rho();
-        const volScalarField& rho2 = mixture_.thermo2().rho();
+        const tmp<volScalarField> trho1 = mixture_.thermo1().rho();
+        const tmp<volScalarField> trho2 = mixture_.thermo2().rho();
+
+        const auto& rho1 = trho1();
+        const auto& rho2 = trho2();
 
         alphaRhoPhi1_ =
         (
@@ -185,8 +188,8 @@ void Foam::compressibleInterPhaseTransportModel::correctPhasePhi()
 {
     if (twoPhaseTransport_)
     {
-        const volScalarField& rho1 = mixture_.thermo1().rho();
-        const volScalarField& rho2 = mixture_.thermo2().rho();
+        const tmp<volScalarField> rho1 = mixture_.thermo1().rho();
+        const tmp<volScalarField> rho2 = mixture_.thermo2().rho();
 
         alphaRhoPhi1_.ref() = fvc::interpolate(rho1)*alphaPhi10_;
         alphaRhoPhi2_.ref() = fvc::interpolate(rho2)*(phi_ - alphaPhi10_);

applications/solvers/multiphase/reactingMultiphaseEulerFoam/EEqns.H

@@ -12,8 +12,11 @@ for (int Ecorr=0; Ecorr<nEnergyCorrectors; Ecorr++)
         phaseModel& phase = fluid.anisothermalPhases()[anisothermalPhasei];
 
         const volScalarField& alpha = phase;
-        const volScalarField& rho = phase.rho();
-        const volVectorField& U = phase.U();
+        const tmp<volScalarField> trho = phase.rho();
+        const tmp<volVectorField> tU = phase.U();
+
+        const auto& rho = trho();
+        const auto& U = tU();
 
         fvScalarMatrix EEqn
         (

applications/solvers/multiphase/reactingMultiphaseEulerFoam/YEqns.H

@@ -11,7 +11,9 @@
 
         UPtrList<volScalarField>& Y = phase.YActiveRef();
         const volScalarField& alpha = phase;
-        const volScalarField& rho = phase.rho();
+        const tmp<volScalarField> trho = phase.rho();
+
+        const auto& rho = trho();
 
         forAll(Y, i)
         {

applications/solvers/multiphase/reactingMultiphaseEulerFoam/pU/UEqns.H

@@ -14,9 +14,11 @@ PtrList<fvVectorMatrix> UEqns(phases.size());
         phaseModel& phase = fluid.movingPhases()[movingPhasei];
 
         const volScalarField& alpha = phase;
-        const volScalarField& rho = phase.rho();
+        const tmp<volScalarField> trho = phase.rho();
         volVectorField& U = phase.URef();
 
+        const auto& rho = trho();
+
         UEqns.set
         (
             phase.index(),

applications/solvers/multiphase/reactingMultiphaseEulerFoam/pUf/UEqns.H

@@ -17,9 +17,11 @@ PtrList<fvVectorMatrix> UEqns(phases.size());
         phaseModel& phase = fluid.movingPhases()[movingPhasei];
 
         const volScalarField& alpha = phase;
-        const volScalarField& rho = phase.rho();
+        const tmp<volScalarField> trho = phase.rho();
         volVectorField& U = phase.URef();
 
+        const auto& rho = trho();
+
         UEqns.set
         (
             phase.index(),

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/createFieldRefs.H

@@ -6,11 +6,13 @@ const volScalarField& alpha2 = phase2;
 
 volVectorField& U1 = phase1.URef();
 surfaceScalarField& phi1 = phase1.phiRef();
-const surfaceScalarField& alphaPhi1 = phase1.alphaPhi();
+const tmp<surfaceScalarField> talphaPhi1 = phase1.alphaPhi();
+const auto& alphaPhi1 = talphaPhi1();
 
 volVectorField& U2 = phase2.URef();
 surfaceScalarField& phi2 = phase2.phiRef();
-const surfaceScalarField& alphaPhi2 = phase2.alphaPhi();
+const tmp<surfaceScalarField> talphaPhi2 = phase2.alphaPhi();
+const auto& alphaPhi2 = talphaPhi2();
 
 surfaceScalarField& phi = fluid.phi();
 

applications/utilities/mesh/conversion/vtkUnstructuredToFoam/Make/options

@@ -1,5 +1,7 @@
 EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude
 
 EXE_LIBS = \
+    -lfiniteVolume \
     -lfileFormats

applications/utilities/mesh/conversion/vtkUnstructuredToFoam/vtkUnstructuredToFoam.C

@@ -34,23 +34,84 @@ Description
     Convert legacy VTK file (ascii) containing an unstructured grid
     to an OpenFOAM mesh without boundary information.
 
+Usage
+    \b vtkUnstructuredToFoam \<XXX.vtk\>
+
+    Options:
+      - \par -no-fields
+        Do not attempt to recreate volFields
+
 Note
     The .vtk format does not contain any boundary information.
-    It is purely a description of the internal mesh.
+    It is purely a description of the internal mesh. This also limits the
+    usefulness of reconstructing the volFields.
+
     Not extensively tested.
 
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
 #include "Time.H"
-#include "polyMesh.H"
+#include "fvMesh.H"
 #include "IFstream.H"
 #include "vtkUnstructuredReader.H"
 
+#include "columnFvMesh.H"
+#include "scalarIOField.H"
+#include "vectorIOField.H"
+#include "volFields.H"
+
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+template<class Type>
+void constructVolFields(fvMesh& mesh, const vtkUnstructuredReader& reader)
+{
+    const auto fields(reader.cellData().csorted<IOField<Type>>());
+    for (const auto& field : fields)
+    {
+        Info<< "Constructing volField " << field.name() << endl;
+
+        // field is
+        //  - cell data followed by
+        //  - boundary face data
+
+
+        auto tfld = GeometricField<Type, fvPatchField, volMesh>::New
+        (
+            field.name(),
+            mesh,
+            dimless
+        );
+        auto& fld = tfld.ref();
+        fld.instance() = mesh.time().timeName();
+        fld.writeOpt() = IOobject::AUTO_WRITE;
+
+        // Fill cell values
+        fld.internalFieldRef().field() =
+            UIndirectList<Type>(field, reader.cellMap());
+
+        // Fill boundary values
+        const auto& map = reader.faceMap();
+        if (map.size())
+        {
+            for (auto& pfld : fld.boundaryFieldRef())
+            {
+                const auto& pp = pfld.patch();
+
+                forAll(pfld, i)
+                {
+                    const label bFacei = pp.patch().offset()+i;
+                    pfld[i] = field[map[bFacei]];
+                }
+            }
+        }
+
+        regIOobject::store(std::move(tfld));
+    }
+}
+
 
 int main(int argc, char *argv[])
 {
@@ -61,11 +122,14 @@ int main(int argc, char *argv[])
     );
 
     argList::noParallel();
+    argList::addOptionCompat("no-fields", {"noFields", 2106});
     argList::addArgument("vtk-file", "The input legacy ascii vtk file");
 
     #include "setRootCase.H"
     #include "createTime.H"
 
+    const bool doFields = !args.found("no-fields");
+
     IFstream mshStream(args.get<fileName>(1));
 
     vtkUnstructuredReader reader(runTime, mshStream);
@@ -96,6 +160,24 @@ int main(int argc, char *argv[])
     mesh.removeFiles();
     mesh.write();
 
+
+    if (doFields)
+    {
+        // Re-read mesh as fvMesh so we can have fields
+        Info<< "Re-reading mesh ..." << endl;
+        #include "createMesh.H"
+
+        constructVolFields<scalar>(mesh, reader);
+        constructVolFields<vector>(mesh, reader);
+        constructVolFields<sphericalTensor>(mesh, reader);
+        constructVolFields<symmTensor>(mesh, reader);
+        constructVolFields<tensor>(mesh, reader);
+
+        // No need to write the mesh, only fields
+        mesh.thisDb().write();
+    }
+
+
     Info<< "End\n" << endl;
 
     return 0;

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C

@@ -743,12 +743,12 @@ Foam::label Foam::cellShapeControlMesh::estimateCellCount
                 cit->vertex(3)->point()
             );
 
-            pointFromPoint centre = topoint(CGAL::centroid(tet));
+            const auto tetCentre = CGAL::centroid(tet);
 
             if
             (
-                Pstream::parRun()
-             && !decomposition().positionOnThisProcessor(centre)
+                UPstream::parRun()
+             && !decomposition().positionOnThisProcessor(topoint(tetCentre))
             )
             {
                 continue;

applications/utilities/mesh/manipulation/transformPoints/transformPoints.C

@@ -83,6 +83,10 @@ Note
     pitch (rotation about y)
     yaw (rotation about z)
 
+    - with -rotate and two exactly opposing vectors it will actually mirror
+    the geometry. Use any of the other rotation options instead or use
+    two steps, each 90 degrees.
+
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -176,7 +176,7 @@ autoPtr<labelIOList> procAddressing
     const objectRegistry& procRegistry,
     const word& name,
     const word& instance,
-    const word& local = fvMesh::meshSubDir
+    const word& local = polyMesh::meshSubDir
 )
 {
     return autoPtr<labelIOList>::New
@@ -218,17 +218,22 @@ const labelIOList& procAddressing
     PtrList<labelIOList>& procAddressingList
 )
 {
-    const fvMesh& procMesh = procMeshList[proci];
+    const auto& procMesh = procMeshList[proci];
 
-    if (!procAddressingList.set(proci))
-    {
-        procAddressingList.set
+    return procAddressingList.try_emplace
+    (
+        proci,
+        IOobject
         (
-            proci,
-            procAddressing(procMesh, name, procMesh.facesInstance())
-        );
-    }
-    return procAddressingList[proci];
+            name,
+            procMesh.facesInstance(),
+            polyMesh::meshSubDir,
+            procMesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            IOobject::NO_REGISTER
+        )
+    );
 }
 
 

src/OpenFOAM/interpolations/interpolateSplineXY/interpolateSplineXY.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -60,11 +61,11 @@ Type Foam::interpolateSplineXY
     label n = xOld.size();
 
     // early exit if out of bounds or only one value
-    if (n == 1 || x < xOld[0])
+    if (n == 1 || x <= xOld[0])
     {
         return yOld[0];
     }
-    if (x > xOld[n - 1])
+    if (x >= xOld[n - 1])
     {
         return yOld[n - 1];
     }

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C

@@ -2443,26 +2443,45 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
     // Seed all boundary faces with owner value. This is to make sure that
     // they are visited (probably only important for coupled faces since
     // these need to be visited from both sides)
+    List<refinementData> nbrCellInfo;
+    syncTools::swapBoundaryCellList(mesh_, allCellInfo, nbrCellInfo);
+
     for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
     {
         // Check if face already handled in loop above
         if (!allFaceInfo[facei].valid(dummyTrackData))
         {
-            label own = faceOwner[facei];
+            const label own = faceOwner[facei];
+            const auto& nbrInfo = nbrCellInfo[facei-mesh_.nInternalFaces()];
 
-            // Seed face with transported data from owner.
-            refinementData faceData;
-            faceData.updateFace
-            (
-                mesh_,
-                facei,
-                own,
-                allCellInfo[own],
-                FaceCellWave<refinementData, int>::propagationTol(),
-                dummyTrackData
-            );
-            seedFaces.append(facei);
-            seedFacesInfo.append(faceData);
+            if (allCellInfo[own].count() > nbrInfo.count())
+            {
+                allFaceInfo[facei].updateFace
+                (
+                    mesh_,
+                    facei,
+                    own,
+                    allCellInfo[own],
+                    FaceCellWave<refinementData, int>::propagationTol(),
+                    dummyTrackData
+                );
+                seedFaces.append(facei);
+                seedFacesInfo.append(allFaceInfo[facei]);
+            }
+            else if (allCellInfo[own].count() < nbrInfo.count())
+            {
+                allFaceInfo[facei].updateFace
+                (
+                    mesh_,
+                    facei,
+                    -1,         // Lucky! neighbCelli not used!
+                    nbrInfo,
+                    FaceCellWave<refinementData, int>::propagationTol(),
+                    dummyTrackData
+                );
+                seedFaces.append(facei);
+                seedFacesInfo.append(allFaceInfo[facei]);
+            }
         }
     }
 

src/dynamicMesh/polyTopoChange/polyTopoChange/refinementData.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2022,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +28,13 @@ License
 
 #include "refinementData.H"
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::refinementData Foam::transform(const tensor&, const refinementData val)
+{
+    return val;
+}
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 Foam::Ostream& Foam::operator<<

src/dynamicMesh/polyTopoChange/polyTopoChange/refinementData.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2020,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -216,6 +216,9 @@ public:
 
 // * * * * * * * * * * * * * * * * * Traits  * * * * * * * * * * * * * * * * //
 
+//- No-op rotational transform for base types
+refinementData transform(const tensor&, const refinementData val);
+
 //- Contiguous data for refinementData
 template<> struct is_contiguous<refinementData> : std::true_type {};
 

src/finiteVolume/fvMesh/wallDist/patchDistMethods/advectionDiffusion/advectionDiffusionPatchDistMethod.C

@@ -79,7 +79,7 @@ Foam::patchDistMethods::advectionDiffusion::advectionDiffusion
     tolerance_(coeffs_.getOrDefault<scalar>("tolerance", 1e-3)),
     maxIter_(coeffs_.getOrDefault<int>("maxIter", 10)),
     predicted_(false),
-    checkAndWriteMesh_(coeffs_.getOrDefault("checkAndWriteMesh", true))
+    checkAndWriteMesh_(coeffs_.getOrDefault("checkAndWriteMesh", false))
 {}
 
 

src/functionObjects/field/Curle/Curle.C

@@ -67,7 +67,8 @@ Foam::functionObjects::Curle::Curle
     c0_(0),
     rawFilePtrs_(),
     inputSurface_(),
-    surfaceWriterPtr_(nullptr)
+    surfaceWriterPtr_(nullptr),
+    warnOnNoPatch_(true)
 {
     read(dict);
 }
@@ -85,8 +86,10 @@ bool Foam::functionObjects::Curle::read(const dictionary& dict)
     }
 
     dict.readIfPresent("p", pName_);
+    dict.readIfPresent("warnOnNoPatch", warnOnNoPatch_);
 
-    patchIDs_ = pbm.patchSet(dict.get<wordRes>("patches")).sortedToc();
+    patchIDs_ =
+        pbm.patchSet(dict.get<wordRes>("patches"), warnOnNoPatch_).sortedToc();
 
     if (patchIDs_.empty())
     {

src/functionObjects/field/Curle/Curle.H

@@ -181,6 +181,9 @@ class Curle
         //- Ouput surface when sampling a surface
         autoPtr<surfaceWriter> surfaceWriterPtr_;
 
+        //- Flag to suppress warnings for missing patches
+        bool warnOnNoPatch_;
+
 
 public:
 

src/functionObjects/field/binField/binModels/binModel/binModel.C

@@ -122,6 +122,7 @@ Foam::binModel::binModel
     mesh_(mesh),
     decomposePatchValues_(false),
     cumulative_(false),
+    warnOnNoPatch_(true),
     coordSysPtr_(nullptr),
     nBin_(1)
 {}
@@ -138,8 +139,11 @@ bool Foam::binModel::read(const dictionary& dict)
         return false;
     }
 
+    dict.readIfPresent("warnOnNoPatch", warnOnNoPatch_);
+
     // Can also use pbm.indices(), but no warnings...
-    patchIDs_ = pbm.patchSet(dict.get<wordRes>("patches")).sortedToc();
+    patchIDs_ =
+        pbm.patchSet(dict.get<wordRes>("patches"), warnOnNoPatch_).sortedToc();
     fieldNames_ = dict.get<wordHashSet>("fields").sortedToc();
 
     wordRes zoneNames;

src/functionObjects/field/binField/binModels/binModel/binModel.H

@@ -77,6 +77,9 @@ protected:
         //- increasing distance in binning direction
         bool cumulative_;
 
+        //- Flag to suppress warnings for missing patches
+        bool warnOnNoPatch_;
+
         //- Local coordinate system of bins
         autoPtr<coordinateSystem> coordSysPtr_;
 

src/functionObjects/field/columnAverage/columnAverage.C

@@ -128,7 +128,8 @@ Foam::functionObjects::columnAverage::columnAverage
 :
     fvMeshFunctionObject(name, runTime, dict),
     patchIDs_(),
-    fieldSet_(mesh_)
+    fieldSet_(mesh_),
+    warnOnNoPatch_(true)
 {
     read(dict);
 }
@@ -140,10 +141,13 @@ bool Foam::functionObjects::columnAverage::read(const dictionary& dict)
 {
     fvMeshFunctionObject::read(dict);
 
+    dict.readIfPresent("warnOnNoPatch", warnOnNoPatch_);
+
     patchIDs_ =
         mesh_.boundaryMesh().patchSet
         (
-            dict.get<wordRes>("patches")
+            dict.get<wordRes>("patches"),
+            warnOnNoPatch_
         ).sortedToc();
 
     fieldSet_.read(dict);

src/functionObjects/field/columnAverage/columnAverage.H

@@ -132,6 +132,9 @@ class columnAverage
         mutable autoPtr<globalIndex> globalPoints_;
         mutable autoPtr<meshStructure> meshStructurePtr_;
 
+        //- Flag to suppress warnings for missing patches
+        bool warnOnNoPatch_;
+
 
     // Private Member Functions
 

src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.C

@@ -51,7 +51,8 @@ Foam::heatTransferCoeffModel::heatTransferCoeffModel
 :
     mesh_(mesh),
     TName_(TName),
-    qrName_("qr")
+    qrName_("qr"),
+    warnOnNoPatch_(true)
 {}
 
 
@@ -150,8 +151,10 @@ bool Foam::heatTransferCoeffModel::read(const dictionary& dict)
     const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
 
     dict.readIfPresent("qr", qrName_);
+    dict.readIfPresent("warnOnNoPatch", warnOnNoPatch_);
 
-    patchIDs_ = pbm.patchSet(dict.get<wordRes>("patches")).sortedToc();
+    patchIDs_ =
+        pbm.patchSet(dict.get<wordRes>("patches"), warnOnNoPatch_).sortedToc();
 
     return true;
 }

src/functionObjects/field/heatTransferCoeff/heatTransferCoeffModels/heatTransferCoeffModel/heatTransferCoeffModel.H

@@ -76,6 +76,9 @@ protected:
         //- Name of radiative heat flux field
         word qrName_;
 
+        //- Flag to suppress warnings for missing patches
+        bool warnOnNoPatch_;
+
 
     // Protected Member Functions
 

src/functionObjects/field/nearWallFields/nearWallFields.C

@@ -245,7 +245,8 @@ Foam::functionObjects::nearWallFields::nearWallFields
 )
 :
     fvMeshFunctionObject(name, runTime, dict),
-    fieldSet_()
+    fieldSet_(),
+    warnOnNoPatch_(true)
 {
     read(dict);
 }
@@ -261,9 +262,11 @@ bool Foam::functionObjects::nearWallFields::read(const dictionary& dict)
 
     dict.readEntry("fields", fieldSet_);
     dict.readEntry("distance", distance_);
+    dict.readIfPresent("warnOnNoPatch", warnOnNoPatch_);
 
     // Can also use pbm.indices(), but no warnings...
-    patchIDs_ = pbm.patchSet(dict.get<wordRes>("patches")).sortedToc();
+    patchIDs_ =
+        pbm.patchSet(dict.get<wordRes>("patches"), warnOnNoPatch_).sortedToc();
 
 
     // Clear out any previously loaded fields

src/functionObjects/field/nearWallFields/nearWallFields.H

@@ -134,6 +134,9 @@ protected:
             //- From resulting back to original field
             HashTable<word> reverseFieldMap_;
 
+            //- Flag to suppress warnings for missing patches
+            bool warnOnNoPatch_;
+
 
         // Calculated addressing
 

src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C

@@ -206,7 +206,10 @@ Foam::functionObjects::regionSizeDistribution::findPatchRegions
 
     labelHashSet patchRegions(2*regions.nRegions());
 
-    labelHashSet patchSet(mesh_.boundaryMesh().patchSet(patchNames_));
+    labelHashSet patchSet
+    (
+        mesh_.boundaryMesh().patchSet(patchNames_, warnOnNoPatch_)
+    );
 
     for (const label patchi : patchSet)
     {
@@ -369,7 +372,8 @@ Foam::functionObjects::regionSizeDistribution::regionSizeDistribution
     writeFile(obr_, name),
     alphaName_(dict.get<word>("field")),
     patchNames_(dict.get<wordRes>("patches")),
-    isoPlanes_(dict.getOrDefault("isoPlanes", false))
+    isoPlanes_(dict.getOrDefault("isoPlanes", false)),
+    warnOnNoPatch_(true)
 {
     read(dict);
 }
@@ -416,6 +420,8 @@ bool Foam::functionObjects::regionSizeDistribution::read(const dictionary& dict)
          direction_.normalise();
     }
 
+    dict.readIfPresent("warnOnNoPatch", warnOnNoPatch_);
+
     return true;
 }
 

src/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H

@@ -249,6 +249,9 @@ class regionSizeDistribution
             //- Switch to enable iso-planes sampling
             bool isoPlanes_;
 
+        //- Flag to suppress warnings for missing patches
+        bool warnOnNoPatch_;
+
 
     // Private Member Functions
 

src/functionObjects/field/wallHeatFlux/wallHeatFlux.C

@@ -106,7 +106,8 @@ Foam::functionObjects::wallHeatFlux::wallHeatFlux
 :
     fvMeshFunctionObject(name, runTime, dict),
     writeFile(obr_, name, typeName, dict),
-    qrName_("qr")
+    qrName_("qr"),
+    warnOnNoPatch_(true)
 {
     read(dict);
 
@@ -144,12 +145,13 @@ bool Foam::functionObjects::wallHeatFlux::read(const dictionary& dict)
     writeFile::read(dict);
 
     dict.readIfPresent("qr", qrName_);
+    dict.readIfPresent("warnOnNoPatch", warnOnNoPatch_);
 
     wordRes patchNames;
     labelHashSet patchSet;
     if (dict.readIfPresent("patches", patchNames) && !patchNames.empty())
     {
-        patchSet = pbm.patchSet(patchNames);
+        patchSet = pbm.patchSet(patchNames, warnOnNoPatch_);
     }
 
     labelHashSet allWalls(pbm.findPatchIDs<wallPolyPatch>());

src/functionObjects/field/wallHeatFlux/wallHeatFlux.H

@@ -123,6 +123,9 @@ protected:
         //- Name of radiative heat flux name
         word qrName_;
 
+        //- Flag to suppress warnings for missing patches
+        bool warnOnNoPatch_;
+
 
     // Protected Member Functions
 

src/functionObjects/field/wallShearStress/wallShearStress.C

@@ -90,7 +90,8 @@ Foam::functionObjects::wallShearStress::wallShearStress
 :
     fvMeshFunctionObject(name, runTime, dict),
     writeFile(mesh_, name, typeName, dict),
-    writeFields_(true)  // May change in the future
+    writeFields_(true),  // May change in the future
+    warnOnNoPatch_(true)
 {
     read(dict);
 
@@ -127,6 +128,7 @@ bool Foam::functionObjects::wallShearStress::read(const dictionary& dict)
 
     writeFields_ = true;   // May change in the future
     dict.readIfPresent("writeFields", writeFields_);
+    dict.readIfPresent("warnOnNoPatch", warnOnNoPatch_);
 
     const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
 
@@ -134,7 +136,7 @@ bool Foam::functionObjects::wallShearStress::read(const dictionary& dict)
     labelHashSet patchSet;
     if (dict.readIfPresent("patches", patchNames) && !patchNames.empty())
     {
-        patchSet = pbm.patchSet(patchNames);
+        patchSet = pbm.patchSet(patchNames, warnOnNoPatch_);
     }
 
     labelHashSet allWalls(pbm.findPatchIDs<wallPolyPatch>());

src/functionObjects/field/wallShearStress/wallShearStress.H

@@ -133,6 +133,9 @@ class wallShearStress
         //- Write the shear stress field ?
         bool writeFields_;
 
+        //- Flag to suppress warnings for missing patches
+        bool warnOnNoPatch_;
+
 
 protected:
 

src/functionObjects/forces/forces/forces.C

@@ -550,7 +550,8 @@ Foam::functionObjects::forces::forces
     directForceDensity_(false),
     porosity_(false),
     writeFields_(false),
-    initialised_(false)
+    initialised_(false),
+    warnOnNoPatch_(true)
 {
     if (readFields)
     {
@@ -589,7 +590,8 @@ Foam::functionObjects::forces::forces
     directForceDensity_(false),
     porosity_(false),
     writeFields_(false),
-    initialised_(false)
+    initialised_(false),
+    warnOnNoPatch_(true)
 {
     if (readFields)
     {
@@ -613,10 +615,13 @@ bool Foam::functionObjects::forces::read(const dictionary& dict)
 
     initialised_ = false;
 
+    dict.readIfPresent("warnOnNoPatch", warnOnNoPatch_);
+
     Info<< type() << ' ' << name() << ':' << endl;
 
     // Can also use pbm.indices(), but no warnings...
-    patchIDs_ = pbm.patchSet(dict.get<wordRes>("patches")).sortedToc();
+    patchIDs_ =
+        pbm.patchSet(dict.get<wordRes>("patches"), warnOnNoPatch_).sortedToc();
 
     dict.readIfPresent("directForceDensity", directForceDensity_);
     if (directForceDensity_)

src/functionObjects/forces/forces/forces.H

@@ -264,6 +264,9 @@ protected:
             //- Flag of initialisation (internal)
             bool initialised_;
 
+            //- Flag to suppress warnings for missing patches
+            bool warnOnNoPatch_;
+
 
     // Protected Member Functions
 

src/functionObjects/solvers/electricPotential/electricPotential.C

@@ -208,6 +208,7 @@ Foam::functionObjects::electricPotential::electricPotential
         )
     ),
     fvOptions_(mesh_),
+    tol_(1),
     nCorr_(1),
     writeDerivedFields_(false),
     electricField_(false)
@@ -258,6 +259,7 @@ bool Foam::functionObjects::electricPotential::read(const dictionary& dict)
     dict.readIfPresent("sigma", sigma_);
     dict.readIfPresent("epsilonr", epsilonr_);
     dict.readIfPresent("nCorr", nCorr_);
+    dict.readIfPresent("tolerance", tol_);
     dict.readIfPresent("writeDerivedFields", writeDerivedFields_);
     dict.readIfPresent("electricField", electricField_);
 
@@ -346,6 +348,10 @@ bool Foam::functionObjects::electricPotential::execute()
 {
     Log << type() << " execute: " << name() << endl;
 
+    // Convergence monitor parameters
+    bool converged = false;
+    label iter = 0;
+
     tmp<volScalarField> tsigma = this->sigma();
     const auto& sigma = tsigma();
 
@@ -362,7 +368,9 @@ bool Foam::functionObjects::electricPotential::execute()
 
         fvOptions_.constrain(eVEqn);
 
-        eVEqn.solve();
+        ++iter;
+        converged = (eVEqn.solve().initialResidual() < tol_);
+        if (converged) break;
     }
 
     if (electricField_)
@@ -371,6 +379,14 @@ bool Foam::functionObjects::electricPotential::execute()
         E == -fvc::grad(eV);
     }
 
+    if (converged)
+    {
+        Log << type() << ": " << name() << ": "
+            << eV.name() << " is converged." << nl
+            << tab << "initial-residual tolerance: " << tol_ << nl
+            << tab << "outer iteration: " << iter << nl;
+    }
+
     Log << endl;
 
     return true;

src/functionObjects/solvers/electricPotential/electricPotential.H

@@ -124,6 +124,7 @@ Usage
         electricField         <bool>;
         E                     <word>;
         fvOptions             <dict>;
+        tolerance             <scalar>;
 
         // Inherited entries
         ...
@@ -143,6 +144,7 @@ Usage
       electricField | Flag to calculate electric field  | bool | no | false
       E    | Name of electric field           | word | no | electricPotential:E
       fvOptions | List of finite-volume options       | dict | no   | -
+      tolerance | Outer-loop initial-residual tolerance | scalar | no | 1
     \endtable
 
     The inherited entries are elaborated in:
@@ -218,6 +220,9 @@ class electricPotential
         //- Run-time selectable finite volume options
         fv::optionList fvOptions_;
 
+        //- Outer-loop initial-residual tolerance
+        scalar tol_;
+
         //- Number of corrector iterations
         int nCorr_;
 

src/functionObjects/solvers/energyTransport/energyTransport.C

@@ -26,11 +26,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "energyTransport.H"
-#include "surfaceFields.H"
-#include "fvmDdt.H"
-#include "fvmDiv.H"
-#include "fvmLaplacian.H"
-#include "fvmSup.H"
 #include "turbulentTransportModel.H"
 #include "turbulentFluidThermoModel.H"
 #include "addToRunTimeSelectionTable.H"
@@ -188,23 +183,6 @@ Foam::functionObjects::energyTransport::energyTransport
 )
 :
     fvMeshFunctionObject(name, runTime, dict),
-    fieldName_(dict.getOrDefault<word>("field", "T")),
-    phiName_(dict.getOrDefault<word>("phi", "phi")),
-    rhoName_(dict.getOrDefault<word>("rho", "rho")),
-    nCorr_(0),
-    schemesField_("unknown-schemesField"),
-    fvOptions_(mesh_),
-    multiphaseThermo_(dict.subOrEmptyDict("phaseThermos")),
-    Cp_("Cp", dimEnergy/dimMass/dimTemperature, 0, dict),
-    kappa_
-    (
-        "kappa",
-        dimEnergy/dimTime/dimLength/dimTemperature,
-        0,
-        dict
-    ),
-    rho_("rhoInf", dimDensity, 0, dict),
-    Prt_("Prt", dimless, 1, dict),
     rhoCp_
     (
         IOobject
@@ -218,7 +196,25 @@ Foam::functionObjects::energyTransport::energyTransport
         ),
         mesh_,
         dimensionedScalar(dimEnergy/dimTemperature/dimVolume, Zero)
-    )
+    ),
+    fvOptions_(mesh_),
+    multiphaseThermo_(dict.subOrEmptyDict("phaseThermos")),
+    Cp_("Cp", dimEnergy/dimMass/dimTemperature, 0, dict),
+    kappa_
+    (
+        "kappa",
+        dimEnergy/dimTime/dimLength/dimTemperature,
+        0,
+        dict
+    ),
+    rho_("rhoInf", dimDensity, 0, dict),
+    Prt_("Prt", dimless, 1, dict),
+    fieldName_(dict.getOrDefault<word>("field", "T")),
+    schemesField_("unknown-schemesField"),
+    phiName_(dict.getOrDefault<word>("phi", "phi")),
+    rhoName_(dict.getOrDefault<word>("rho", "rho")),
+    tol_(1),
+    nCorr_(0)
 {
     read(dict);
 
@@ -305,17 +301,14 @@ Foam::functionObjects::energyTransport::energyTransport
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::functionObjects::energyTransport::~energyTransport()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::functionObjects::energyTransport::read(const dictionary& dict)
 {
-    fvMeshFunctionObject::read(dict);
+    if (!fvMeshFunctionObject::read(dict))
+    {
+        return false;
+    }
 
     dict.readIfPresent("phi", phiName_);
     dict.readIfPresent("rho", rhoName_);
@@ -323,6 +316,7 @@ bool Foam::functionObjects::energyTransport::read(const dictionary& dict)
     schemesField_ = dict.getOrDefault("schemesField", fieldName_);
 
     dict.readIfPresent("nCorr", nCorr_);
+    dict.readIfPresent("tolerance", tol_);
 
     if (dict.found("fvOptions"))
     {
@@ -355,12 +349,16 @@ bool Foam::functionObjects::energyTransport::execute()
     scalar relaxCoeff = 0;
     mesh_.relaxEquation(schemesField_, relaxCoeff);
 
+    // Convergence monitor parameters
+    bool converged = false;
+    int iter = 0;
+
     if (phi.dimensions() == dimMass/dimTime)
     {
         rhoCp_ = rho()*Cp();
         const surfaceScalarField rhoCpPhi(fvc::interpolate(Cp())*phi);
 
-        for (label i = 0; i <= nCorr_; i++)
+        for (int i = 0; i <= nCorr_; ++i)
         {
             fvScalarMatrix sEqn
             (
@@ -376,7 +374,9 @@ bool Foam::functionObjects::energyTransport::execute()
 
             fvOptions_.constrain(sEqn);
 
-            sEqn.solve(schemesField_);
+            ++iter;
+            converged = (sEqn.solve(schemesField_).initialResidual() < tol_);
+            if (converged) break;
         }
     }
     else if (phi.dimensions() == dimVolume/dimTime)
@@ -393,7 +393,7 @@ bool Foam::functionObjects::energyTransport::execute()
             rhoCp
         );
 
-        for (label i = 0; i <= nCorr_; i++)
+        for (int i = 0; i <= nCorr_; ++i)
         {
             fvScalarMatrix sEqn
             (
@@ -408,7 +408,9 @@ bool Foam::functionObjects::energyTransport::execute()
 
             fvOptions_.constrain(sEqn);
 
-            sEqn.solve(schemesField_);
+            ++iter;
+            converged = (sEqn.solve(schemesField_).initialResidual() < tol_);
+            if (converged) break;
         }
     }
     else
@@ -419,6 +421,14 @@ bool Foam::functionObjects::energyTransport::execute()
             << dimVolume/dimTime << exit(FatalError);
     }
 
+    if (converged)
+    {
+        Log << type() << ": " << name() << ": "
+            << s.name() << " is converged." << nl
+            << tab << "initial-residual tolerance: " << tol_ << nl
+            << tab << "outer iteration: " << iter << nl;
+    }
+
     Log << endl;
 
     return true;

src/functionObjects/solvers/energyTransport/energyTransport.H

@@ -30,151 +30,180 @@ Group
     grpSolversFunctionObjects
 
 Description
-    Evolves a simplified energy transport equation for incompressible flows.
-    It takes into account the inertia, conduction and convection terms plus
-    a source.
+    Computes the simplified energy transport equation in single-phase or
+    two-phase flow, considering incompressible cases:
 
-    - The field name must be temperature and its BC's specified in the time
-      directory.
-    - The turbulence model should be incompressible
-    - In order to use in a incompressible multi phase a list of thermal
-      properties are needed. See below
+    \f[
+        \frac{\partial \rho \, C_p \, T}{\partial t}
+      + \nabla \cdot \left(\rho \, C_p \, \phi \, T \right)
+      - \nabla \cdot \left(\rho \, C_p \, \phi \right) \, T
+      - \nabla \cdot \left(\kappa_{eff} \, \nabla T \right)
+      = S_T
+    \f]
 
+    where:
+    \vartable
+      T       | Scalar field
+      \rho    | (Generic) Fluid density which is unity when not specified
+      C_p     | Specific heat capacity at constant pressure
+      \phi    | (Generic) Flux field
+      \kappa_{eff}  | Effective thermal conductivity
+      S_T     | Scalar field source term
+    \endvartable
 
 Usage
-    Example of function object specification to solve a energy transport
-    equation for a single phase flow plus a source term
+    Minimal example in \c system/controlDict.functions:
     \verbatim
-    functions
+    energyTransport1
     {
-        energy
-        {
-            type            energyTransport;
-            libs            (energyTransportFunctionObjects);
+        // Mandatory entries
+        type            energyTransport;
+        libs            (solverFunctionObjects);
 
-            enabled         true;
-            writeControl    writeTime;
-            writeInterval   1;
+        // Optional entries
+        field           <word>;
+        phi             <word>;
+        rho             <word>;
+        Cp              <scalar>;
+        kappa           <scalar>;
+        rhoInf          <scalar>;
+        Prt             <scalar>;
+        schemesField    <word>;
+        tolerance       <scalar>;
+        nCorr           <int>;
+        fvOptions       <dict>;
+        phaseThermos    <dict>;
 
-            field           T;
-
-            // volumetric Flux
-            phi             phi;
-
-            // Thermal properties
-            Cp              Cp    [J/kg/K]  1e3;
-            kappa           kappa [W/m/K]   0.0257;
-            rhoInf          rho   [kg/m^3]  1.2;
-
-            write           true;
-
-            fvOptions
-            {
-                viscousDissipation
-                {
-                    type            viscousDissipation;
-                    enabled         true;
-
-                    viscousDissipationCoeffs
-                    {
-                        fields          (T);
-                        rhoInf          $....rhoInf;
-                    }
-                }
-            }
-        }
+        // Inherited entries
+        ...
     }
     \endverbatim
 
-    Example of function object specification to solve a energy transport
-    equation for a multiphase phase flow plus a source term
-
-    equation:
-    \verbatim
-    functions
-    {
-        energy
-        {
-            type            energyTransport;
-            libs            (energyTransportFunctionObjects);
-
-            enabled         true;
-            writeControl    writeTime;
-            writeInterval   1;
-
-            field           T;
-
-            // rho field name
-            rho             rho;
-            // mass flux for multiphase
-            phi             rhoPhi;
-
-            write           true;
-
-            // Thermal properties of the phases
-            phaseThermos
-            {
-                alpha.air
-                {
-                    Cp          1e3;
-                    kappa       0.0243;
-                }
-                alpha.mercury
-                {
-                    Cp          140;
-                    kappa       8.2;
-                }
-                alpha.oil
-                {
-                    Cp          2e3;
-                    kappa       0.2;
-                }
-                alpha.water
-                {
-                    Cp          4e3;
-                    kappa       0.6;
-                }
-            }
-
-
-            fvOptions
-            {
-                viscousDissipation
-                {
-                    type            viscousDissipation;
-                    enabled         true;
-
-                    viscousDissipationCoeffs
-                    {
-                        fields          (T);
-                        rho             rho; //rho Field
-                    }
-                }
-            }
-        }
-    }
-    \endverbatim
-
-    Where the entries comprise:
+    where:
     \table
-        Property     | Description             | Required    | Default value
-        type         | Type name: energyTransport | yes      |
-        field        | Name of the scalar field | no         | T
-        phi          | Name of flux field      | no          | phi
-        rho          | Name of density field   | no          | rho
-        nCorr        | Number of correctors    | no          | 0
-        schemesField | Name of field to specify schemes | no | field name
-        fvOptions    | List of scalar sources  | no          |
-        Cp           | Heat capacity for single phase | no   | 0
-        rhoInf       | Density for single phase | no         | 0
-        kappa        | Thermal conductivity for single phase | no   | 0
-        Prt          | Turbulent Prandlt number | no         | 1.0
-        phaseThermos | Dictionary for multi-phase thermo |no | null
-        fvOptions    | Opotional extra sources | no          | null
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: energyTransport          | word | yes  | -
+      libs      | Library name: solverFunctionObjects | word | yes  | -
+      field     | Name of the passive-scalar field    | word | no   | s
+      phi       | Name of flux field                  | word | no   | phi
+      rho       | Name of density field               | word | no   | rho
+      Cp        | Specific heat capacity at constant pressure | scalar | no | 0
+      kappa     | Thermal conductivity                | scalar | no | 0
+      rhoInf    | Fluid density                       | scalar | no | 0
+      Prt       | Turbulent Prandtl number            | scalar | no | 1
+      schemesField | Name of field to specify schemes | word   | no | field
+      tolerance | Outer-loop initial-residual tolerance | scalar | no | 1
+      nCorr     | Number of outer-loop correctors     | int    | no | 0
+      fvOptions | List of finite-volume options       | dict   | no | -
+      phaseThermos | Dictionary for multi-phase thermo | dict  | no | null
     \endtable
 
-See also
-    Foam::functionObjects::fvMeshFunctionObject
+    The inherited entries are elaborated in:
+      - \link fvMeshFunctionObject.H \endlink
+      - \link fvOption.H \endlink
+
+    An example of function object specification to solve a energy transport
+    equation for a single phase flow plus a source term:
+    \verbatim
+    energyTransport1
+    {
+        // Mandatory entries
+        type            energyTransport;
+        libs            (solverFunctionObjects);
+
+        // Optional entries
+        field           T;
+        phi             phi;
+        Cp              Cp    [J/kg/K]  1e3;
+        kappa           kappa [W/m/K]   0.0257;
+        rhoInf          rho   [kg/m^3]  1.2;
+        fvOptions
+        {
+            viscousDissipation
+            {
+                type            viscousDissipation;
+                enabled         true;
+
+                viscousDissipationCoeffs
+                {
+                    fields          (T);
+                    rhoInf          $....rhoInf;
+                }
+            }
+        }
+
+        // Inherited entries
+        enabled         true;
+        writeControl    writeTime;
+        writeInterval   1;
+    }
+    \endverbatim
+
+    An example of function object specification to solve a energy transport
+    equation for a multiphase phase flow plus a source term:
+    \verbatim
+    energyTransport1
+    {
+        // Mandatory entries
+        type            energyTransport;
+        libs            (solverFunctionObjects);
+
+        // Optional entries
+        field           T;
+        rho             rho;
+        phi             rhoPhi;
+
+        // Thermal properties of the phases
+        phaseThermos
+        {
+            alpha.air
+            {
+                Cp          1e3;
+                kappa       0.0243;
+            }
+            alpha.mercury
+            {
+                Cp          140;
+                kappa       8.2;
+            }
+            alpha.oil
+            {
+                Cp          2e3;
+                kappa       0.2;
+            }
+            alpha.water
+            {
+                Cp          4e3;
+                kappa       0.6;
+            }
+        }
+
+        fvOptions
+        {
+            viscousDissipation
+            {
+                type            viscousDissipation;
+                enabled         true;
+
+                viscousDissipationCoeffs
+                {
+                    fields          (T);
+                    rho             rho;
+                }
+            }
+        }
+
+        // Inherited entries
+        enabled         true;
+        writeControl    writeTime;
+        writeInterval   1;
+    }
+    \endverbatim
+
+Note
+  - The field name must be temperature and its boundary conditions
+    specified in the time directory.
+  - The turbulence model should be incompressible.
 
 SourceFiles
     energyTransport.C
@@ -203,22 +232,10 @@ class energyTransport
 :
     public fvMeshFunctionObject
 {
-    // Private data
+    // Private Data
 
-        //- Name of the transport field.
-        word fieldName_;
-
-        //- Name of flux field
-        word phiName_;
-
-        //- Name of density field
-        word rhoName_;
-
-        //- Number of corrector iterations (optional)
-        label nCorr_;
-
-        //- Name of field whose schemes are used (optional)
-        word schemesField_;
+        //- Volumetric heat capacity field [J/m^3/K]
+        volScalarField rhoCp_;
 
         //- Run-time selectable finite volume options, e.g. sources, constraints
         fv::optionList fvOptions_;
@@ -229,7 +246,7 @@ class energyTransport
         //- List of phase names
         wordList phaseNames_;
 
-        //- List of phase heat capacities
+        //- List of phase specific heat capacities at constant pressure
         PtrList<dimensionedScalar> Cps_;
 
         //- List of phase thermal diffusivity for temperature [J/m/s/K]
@@ -238,7 +255,7 @@ class energyTransport
         //- Unallocated phase list
         UPtrList<volScalarField> phases_;
 
-        //- Heat capacity for single phase flows
+        //- Specific heat capacity at constant pressure for single phase flows
         dimensionedScalar Cp_;
 
         //- Thermal diffusivity for temperature for single phase flows
@@ -250,8 +267,23 @@ class energyTransport
         //- Turbulent Prandt number
         dimensionedScalar Prt_;
 
-        //- rhoCp
-        volScalarField rhoCp_;
+        //- Name of the transport field
+        word fieldName_;
+
+        //- Name of field whose schemes are used
+        word schemesField_;
+
+        //- Name of flux field
+        word phiName_;
+
+        //- Name of density field
+        word rhoName_;
+
+        //- Outer-loop initial-residual tolerance
+        scalar tol_;
+
+        //- Number of corrector iterations
+        int nCorr_;
 
 
     // Private Member Functions
@@ -262,21 +294,15 @@ class energyTransport
         //- Return the diffusivity field
         tmp<volScalarField> kappaEff() const;
 
-        //- Return rho field
+        //- Return the density field, rho
         tmp<volScalarField> rho() const;
 
-        //- Return Cp
+        //- Return the specific heat capacity at constant pressure field, Cp
         tmp<volScalarField> Cp() const;
 
-        //- Return kappa
+        //- Return the thermal diffusivity field
         tmp<volScalarField> kappa() const;
 
-        //- No copy construct
-        energyTransport(const energyTransport&) = delete;
-
-        //- No copy assignment
-        void operator=(const energyTransport&) = delete;
-
 
 public:
 
@@ -296,7 +322,7 @@ public:
 
 
     //- Destructor
-    virtual ~energyTransport();
+    virtual ~energyTransport() = default;
 
 
     // Member Functions

src/functionObjects/solvers/scalarTransport/scalarTransport.C

@@ -27,11 +27,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "scalarTransport.H"
-#include "surfaceFields.H"
-#include "fvmDdt.H"
-#include "fvmDiv.H"
-#include "fvmLaplacian.H"
-#include "fvmSup.H"
 #include "CMULES.H"
 #include "turbulentTransportModel.H"
 #include "turbulentFluidThermoModel.H"
@@ -172,7 +167,9 @@ Foam::functionObjects::scalarTransport::scalarTransport
 )
 :
     fvMeshFunctionObject(name, runTime, dict),
+    fvOptions_(mesh_),
     fieldName_(dict.getOrDefault<word>("field", "s")),
+    schemesField_("unknown-schemesField"),
     phiName_(dict.getOrDefault<word>("phi", "phi")),
     rhoName_(dict.getOrDefault<word>("rho", "rho")),
     nutName_(dict.getOrDefault<word>("nut", "none")),
@@ -182,11 +179,12 @@ Foam::functionObjects::scalarTransport::scalarTransport
         dict.getOrDefault<word>("phasePhiCompressed", "alphaPhiUn")
     ),
     D_(0),
-    constantD_(false),
+    alphaD_(1),
+    alphaDt_(1),
+    tol_(1),
     nCorr_(0),
     resetOnStartUp_(false),
-    schemesField_("unknown-schemesField"),
-    fvOptions_(mesh_),
+    constantD_(false),
     bounded01_(dict.getOrDefault("bounded01", true))
 {
     read(dict);
@@ -202,31 +200,30 @@ Foam::functionObjects::scalarTransport::scalarTransport
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::functionObjects::scalarTransport::~scalarTransport()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::functionObjects::scalarTransport::read(const dictionary& dict)
 {
-    fvMeshFunctionObject::read(dict);
+    if (!fvMeshFunctionObject::read(dict))
+    {
+        return false;
+    }
 
     dict.readIfPresent("phi", phiName_);
     dict.readIfPresent("rho", rhoName_);
     dict.readIfPresent("nut", nutName_);
     dict.readIfPresent("phase", phaseName_);
-    dict.readIfPresent("bounded01", bounded01_);
+    dict.readIfPresent("phasePhiCompressed", phasePhiCompressedName_);
 
     schemesField_ = dict.getOrDefault("schemesField", fieldName_);
-    constantD_ = dict.readIfPresent("D", D_);
-    alphaD_ = dict.getOrDefault<scalar>("alphaD", 1);
-    alphaDt_ = dict.getOrDefault<scalar>("alphaDt", 1);
 
+    dict.readIfPresent("alphaD", alphaD_);
+    dict.readIfPresent("alphaDt", alphaDt_);
+    dict.readIfPresent("tolerance", tol_);
     dict.readIfPresent("nCorr", nCorr_);
     dict.readIfPresent("resetOnStartUp", resetOnStartUp_);
+    constantD_ = dict.readIfPresent("D", D_);
+    dict.readIfPresent("bounded01", bounded01_);
 
     if (dict.found("fvOptions"))
     {
@@ -256,6 +253,10 @@ bool Foam::functionObjects::scalarTransport::execute()
     scalar relaxCoeff = 0;
     mesh_.relaxEquation(schemesField_, relaxCoeff);
 
+    // Convergence monitor parameters
+    bool converged = false;
+    int iter = 0;
+
     // Two phase scalar transport
     if (phaseName_ != "none")
     {
@@ -272,7 +273,7 @@ bool Foam::functionObjects::scalarTransport::execute()
 
         // Solve
         tmp<surfaceScalarField> tTPhiUD;
-        for (label i = 0; i <= nCorr_; i++)
+        for (int i = 0; i <= nCorr_; ++i)
         {
             fvScalarMatrix sEqn
             (
@@ -285,9 +286,13 @@ bool Foam::functionObjects::scalarTransport::execute()
 
             sEqn.relax(relaxCoeff);
             fvOptions_.constrain(sEqn);
-            sEqn.solve(schemesField_);
+
+            ++iter;
+            converged = (sEqn.solve(schemesField_).initialResidual() < tol_);
 
             tTPhiUD = sEqn.flux();
+
+            if (converged) break;
         }
 
         if (bounded01_)
@@ -307,9 +312,8 @@ bool Foam::functionObjects::scalarTransport::execute()
     {
         const volScalarField& rho = lookupObject<volScalarField>(rhoName_);
 
-        for (label i = 0; i <= nCorr_; i++)
+        for (int i = 0; i <= nCorr_; ++i)
         {
-
             fvScalarMatrix sEqn
             (
                 fvm::ddt(rho, s)
@@ -323,12 +327,14 @@ bool Foam::functionObjects::scalarTransport::execute()
 
             fvOptions_.constrain(sEqn);
 
-            sEqn.solve(schemesField_);
+            ++iter;
+            converged = (sEqn.solve(schemesField_).initialResidual() < tol_);
+            if (converged) break;
         }
     }
     else if (phi.dimensions() == dimVolume/dimTime)
     {
-        for (label i = 0; i <= nCorr_; i++)
+        for (int i = 0; i <= nCorr_; ++i)
         {
             fvScalarMatrix sEqn
             (
@@ -343,7 +349,9 @@ bool Foam::functionObjects::scalarTransport::execute()
 
             fvOptions_.constrain(sEqn);
 
-            sEqn.solve(schemesField_);
+            ++iter;
+            converged = (sEqn.solve(schemesField_).initialResidual() < tol_);
+            if (converged) break;
         }
     }
     else
@@ -354,6 +362,14 @@ bool Foam::functionObjects::scalarTransport::execute()
             << dimVolume/dimTime << exit(FatalError);
     }
 
+    if (converged)
+    {
+        Log << type() << ": " << name() << ": "
+            << s.name() << " is converged." << nl
+            << tab << "initial-residual tolerance: " << tol_ << nl
+            << tab << "outer iteration: " << iter << nl;
+    }
+
     Log << endl;
 
     return true;

src/functionObjects/solvers/scalarTransport/scalarTransport.H

@@ -31,110 +31,134 @@ Group
     grpSolversFunctionObjects
 
 Description
-    Evolves a passive scalar transport equation.
+    Computes the transport equation for a passive scalar in single-phase or
+    two-phase flow, considering both incompressible and compressible cases:
 
-    - To specify the field name set the \c field entry
-    - To employ the same numerical schemes as another field set
-      the \c schemesField entry,
-    - The diffusivity can be set manually using the 'D' entry, retrieved
-      from the turbulence model or specified nut
-    - Alternatively if a turbulence model is available a turbulent diffusivity
-      may be constructed from the laminar and turbulent viscosities using the
-      optional diffusivity coefficients \c alphaD and \c alphaDt (which default
-      to 1):
-      \verbatim
-          D = alphaD*nu + alphaDt*nut
-      \endverbatim
-    - To specify a transport quantity within a phase enter phase.
-    - bounded01 bounds the transported scalar within 0 and 1.
+    \f[
+        \frac{\partial \rho \, T}{\partial t}
+      + \nabla \cdot \left( \phi_\alpha \, T \right)
+      - \nabla \cdot (D_T \, \nabla T)
+      = \alpha \, S_T
+    \f]
+
+    where:
+    \vartable
+      T       | Passive scalar field
+      \rho    | (Generic) Fluid density which is unity when not specified
+      \phi_\alpha | (Generic) Flux field
+      \alpha  | Phase fraction which is unity for single-phase flows
+      D_T     | Diffusivity representing the diffusive transport of T
+      S_T     | Passive-scalar field source term
+    \endvartable
 
 Usage
-    Example of function object specification to solve a scalar transport
-    equation:
+    Minimal example in \c system/controlDict.functions:
     \verbatim
-    functions
+    scalarTransport1
     {
-        scalar1
-        {
-            type            scalarTransport;
-            libs            (solverFunctionObjects);
+        // Mandatory entries
+        type            scalarTransport;
+        libs            (solverFunctionObjects);
 
-            resetOnStartUp  no;
-            region          cabin;
-            field           H2O;
+        // Optional entries
+        field           <word>;
+        phi             <word>;
+        rho             <word>;
+        nut             <word>;
+        phase           <word>;
+        phasePhiCompressed  <word>;
+        schemesField    <word>;
+        bounded01       <bool>;
+        D               <scalar>;
+        alphaD          <scalar>;
+        alphaDt         <scalar>;
+        tolerance       <scalar>;
+        nCorr           <int>;
+        resetOnStartUp  <bool>;
+        fvOptions       <dict>;
 
-
-            fvOptions
-            {
-                ...
-            }
-        }
+        // Inherited entries
+        ...
     }
-    \endverbatim
 
-    Example of function object specification to solve a residence time
-    in a two phase flow:
-    equation:
-    \verbatim
-    functions
-    {
-        sTransport
-        {
-            type            scalarTransport;
-            libs            (solverFunctionObjects);
-
-            enabled         true;
-            writeControl    writeTime;
-            writeInterval   1;
-
-            field           s;
-            bounded01       false;
-            phase           alpha.water;
-
-            write           true;
-
-            fvOptions
-            {
-                unitySource
-                {
-                    type            scalarSemiImplicitSource;
-                    enabled         true;
-
-                    selectionMode   all;
-                    volumeMode      specific;
-
-                    sources
-                    {
-                        s           (1 0);
-                    }
-                }
-            }
-
-            resetOnStartUp  false;
-        }
-    }
-    \endverbatim
-
-    Where the entries comprise:
+    where:
     \table
-        Property     | Description             | Required    | Default value
-        type         | Type name: scalarTransport | yes      |
-        field        | Name of the scalar field | no         | s
-        phi          | Name of flux field      | no          | phi
-        rho          | Name of density field   | no          | rho
-        phase        | Name of the phase       | no          | none
-        nut          | Name of the turbulence viscosity | no | none
-        D            | Diffusion coefficient   | no          | auto generated
-        nCorr        | Number of correctors    | no          | 0
-        resetOnStartUp | Reset scalar to zero on start-up | no | no
-        schemesField | Name of field to specify schemes | no | field name
-        fvOptions    | List of scalar sources  | no          |
-        bounded01    | Bounds scalar between 0-1 for multiphase | no | true
-        phasePhiCompressed | Compressed flux for VOF | no | alphaPhiUn
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: scalarTransport          | word | yes  | -
+      libs      | Library name: solverFunctionObjects | word | yes  | -
+      field     | Name of the passive-scalar field    | word | no   | s
+      phi       | Name of flux field                  | word | no   | phi
+      rho       | Name of density field               | word | no   | rho
+      nut       | Name of the turbulence viscosity    | word | no   | none
+      phase     | Name of the phase                   | word | no   | none
+      phasePhiCompressed | Name of compressed VOF flux | word | no | alphaPhiUn
+      schemesField | Name of field to specify schemes | word  | no | field
+      bounded01 | Bounds scalar between 0-1 for multiphase  | bool | no | true
+      D         | Diffusion coefficient               | scalar | no | -
+      alphaD    | Laminar diffusivity coefficient     | scalar | no | 1
+      alphaDt   | Turbulent diffusivity coefficient   | scalar | no | 1
+      tolerance | Outer-loop initial-residual tolerance | scalar | no | 1
+      nCorr     | Number of outer-loop correctors     | int    | no | 0
+      resetOnStartUp | Flag to reset field to zero on start-up | bool | no | no
+      fvOptions | List of finite-volume options       | dict   | no | -
     \endtable
 
-See also
-    Foam::functionObjects::fvMeshFunctionObject
+    The inherited entries are elaborated in:
+      - \link fvMeshFunctionObject.H \endlink
+      - \link fvOption.H \endlink
+
+    An example of function object specification to solve a residence time
+    in a two-phase flow:
+    \verbatim
+    scalarTransport1
+    {
+        // Mandatory entries
+        type            scalarTransport;
+        libs            (solverFunctionObjects);
+
+        // Optional entries
+        field           s;
+        bounded01       false;
+        phase           alpha.water;
+        tolerance       1e-5;
+        resetOnStartUp  false;
+        fvOptions
+        {
+            unitySource
+            {
+                type            scalarSemiImplicitSource;
+                enabled         true;
+
+                selectionMode   all;
+                volumeMode      specific;
+
+                sources
+                {
+                    s           (1 0);
+                }
+            }
+        }
+
+        // Inherited entries
+        enabled         true;
+        writeControl    writeTime;
+        writeInterval   1;
+    }
+    \endverbatim
+
+Note
+  - To use the same numerical schemes as another field,
+  set the \c schemesField entry.
+  - The diffusivity can be set manually using the \c D entry, obtained
+  from the turbulence model or specified as `nut`.
+  - Alternatively, if a turbulence model is available, turbulent diffusivity
+  can be constructed from the laminar and turbulent viscosities using the
+  optional diffusivity coefficients \c alphaD and \c alphaDt
+  (which default to 1):
+
+    \f[
+        D = \alpha_D \, \nu + \alpha_{Dt} \, \nu_t
+    \f]
 
 SourceFiles
     scalarTransport.C
@@ -163,49 +187,52 @@ class scalarTransport
 :
     public fvMeshFunctionObject
 {
-    // Private data
+    // Private Data
+
+        //- Run-time selectable finite volume options, e.g. sources, constraints
+        fv::optionList fvOptions_;
 
         //- Name of the transport field.
         word fieldName_;
 
-        //- Name of flux field (optional)
+        //- Name of field whose schemes are used
+        word schemesField_;
+
+        //- Name of flux field
         word phiName_;
 
-        //- Name of density field (optional)
+        //- Name of density field
         word rhoName_;
 
-        //- Name of turbulent viscosity field (optional)
+        //- Name of turbulent viscosity field
         word nutName_;
 
-        //- Name of phase field (optional)
+        //- Name of phase field
         word phaseName_;
 
-        //- Name of phase field compressed flux (optional)
+        //- Name of phase field compressed flux
         word phasePhiCompressedName_;
 
-        //- Diffusion coefficient (optional)
+        //- Diffusion coefficient
         scalar D_;
 
-        //- Flag to indicate whether a constant, uniform D_ is specified
-        bool constantD_;
-
-        //- Laminar diffusion coefficient (optional)
+        //- Laminar diffusion coefficient
         scalar alphaD_;
 
-        //- Turbulent diffusion coefficient (optional)
+        //- Turbulent diffusion coefficient
         scalar alphaDt_;
 
-        //- Number of corrector iterations (optional)
-        label nCorr_;
+        //- Outer-loop initial-residual tolerance
+        scalar tol_;
+
+        //- Number of corrector iterations
+        int nCorr_;
 
         //- Flag to reset the scalar to zero on start-up
         bool resetOnStartUp_;
 
-        //- Name of field whose schemes are used (optional)
-        word schemesField_;
-
-        //- Run-time selectable finite volume options, e.g. sources, constraints
-        fv::optionList fvOptions_;
+        //- Flag to indicate whether a constant, uniform D_ is specified
+        bool constantD_;
 
         //- Bound scalar between 0-1 using MULES for multiphase case
         bool bounded01_;
@@ -223,12 +250,6 @@ class scalarTransport
             const surfaceScalarField& phi
         ) const;
 
-        //- No copy construct
-        scalarTransport(const scalarTransport&) = delete;
-
-        //- No copy assignment
-        void operator=(const scalarTransport&) = delete;
-
 
 public:
 
@@ -248,7 +269,7 @@ public:
 
 
     //- Destructor
-    virtual ~scalarTransport();
+    virtual ~scalarTransport() = default;
 
 
     // Member Functions

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C

@@ -2054,7 +2054,6 @@ Foam::labelList Foam::meshRefinement::intersectedPoints() const
 
     // Mark all points on faces that will become baffles
     bitSet isBoundaryPoint(mesh_.nPoints());
-    label nBoundaryPoints = 0;
 
     const labelList& surfIndex = surfaceIndex();
 
@@ -2062,15 +2061,7 @@ Foam::labelList Foam::meshRefinement::intersectedPoints() const
     {
         if (surfIndex[facei] != -1)
         {
-            const face& f = faces[facei];
-
-            forAll(f, fp)
-            {
-                if (isBoundaryPoint.set(f[fp]))
-                {
-                    nBoundaryPoints++;
-                }
-            }
+            isBoundaryPoint.set(faces[facei]);
         }
     }
 
@@ -2083,37 +2074,17 @@ Foam::labelList Foam::meshRefinement::intersectedPoints() const
     //    if (patchi != -1)
     //    {
     //        const polyPatch& pp = mesh_.boundaryMesh()[patchi];
-    //
-    //        label facei = pp.start();
-    //
     //        forAll(pp, i)
     //        {
-    //            const face& f = faces[facei];
-    //
-    //            forAll(f, fp)
-    //            {
-    //                if (isBoundaryPoint.set(f[fp]))
-    //                    nBoundaryPoints++;
-    //                }
-    //            }
-    //            facei++;
+    //            isBoundaryPoint.set(faces[pp.start()+i]);
     //        }
     //    }
     //}
 
+    // Make sure all processors have the same data
+    syncTools::syncPointList(mesh_, isBoundaryPoint, orEqOp<unsigned int>(), 0);
 
-    // Pack
-    labelList boundaryPoints(nBoundaryPoints);
-    nBoundaryPoints = 0;
-    forAll(isBoundaryPoint, pointi)
-    {
-        if (isBoundaryPoint.test(pointi))
-        {
-            boundaryPoints[nBoundaryPoints++] = pointi;
-        }
-    }
-
-    return boundaryPoints;
+    return isBoundaryPoint.sortedToc();
 }
 
 

src/meshTools/cellDist/cellDistFuncs.C

@@ -256,7 +256,7 @@ void Foam::cellDistFuncs::correctBoundaryFaceCells
             {
                 if (areaFraction && (areaFraction()[patchFacei] <= 0.5))
                 {
-                    // is mostly coupled
+                    // For cyclicACMI: more cyclic than wall
                     continue;
                 }
 
@@ -311,9 +311,9 @@ void Foam::cellDistFuncs::correctBoundaryPointCells
                 // Check cells with face on wall
                 forAll(patch, patchFacei)
                 {
-                    if (!areaFraction || (areaFraction()[patchFacei] <= 0.5))
+                    if (areaFraction && (areaFraction()[patchFacei] <= 0.5))
                     {
-                        // face mostly coupled
+                        // For cyclicACMI: more cyclic than wall
                         isWallPoint.unset(localFaces[patchFacei]);
                     }
                 }

src/meshTools/coordSet/writers/nastran/nastranCoordSetWriter.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -250,12 +250,12 @@ Foam::fileName Foam::coordSetWriters::nastranWriter::writeTemplate
             Info<< "Writing nastran geometry to " << outputFile << endl;
         }
 
-        if (!isDir(outputFile.path()))
+        if (!Foam::isDir(outputFile.path()))
         {
-            mkDir(outputFile.path());
+            Foam::mkDir(outputFile.path());
         }
 
-        OFstream os(outputFile);
+        OFstream os(IOstreamOption::ATOMIC, outputFile);
         fileFormats::NASCore::setPrecision(os, writeFormat_);
 
         os  << "TITLE=OpenFOAM " << outputFile.stem()
@@ -293,12 +293,12 @@ Foam::fileName Foam::coordSetWriters::nastranWriter::writeTemplate
             Info<< "Writing nastran geometry to " << outputFile << endl;
         }
 
-        if (!isDir(outputFile.path()))
+        if (!Foam::isDir(outputFile.path()))
         {
-            mkDir(outputFile.path());
+            Foam::mkDir(outputFile.path());
         }
 
-        OFstream os(outputFile);
+        OFstream os(IOstreamOption::ATOMIC, outputFile);
         fileFormats::NASCore::setPrecision(os, writeFormat_);
 
         os  << "TITLE=OpenFOAM " << outputFile.stem()

src/overset/cellCellStencil/cellCellStencil/cellCellStencilObject.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,7 +26,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "cellCellStencilObject.H"
-
+#include "dynamicOversetFvMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -36,4 +36,25 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::cellCellStencilObject::movePoints()
+{
+    // Update (if needed) underlying stencil
+    const bool changed = stencilPtr_().update();
+
+    if (changed)
+    {
+        const auto* oversetMeshPtr = isA<oversetFvMeshBase>(mesh());
+        if (oversetMeshPtr)
+        {
+            // Clear out any additional (ldu)addressing
+            const_cast<oversetFvMeshBase&>(*oversetMeshPtr).clearOut();
+        }
+    }
+
+    return changed;
+}
+
+
 // ************************************************************************* //

src/overset/cellCellStencil/cellCellStencil/cellCellStencilObject.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2019 OpenCFD Ltd.
+    Copyright (C) 2017-2019,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,19 +54,9 @@ typedef MeshObject<fvMesh, MoveableMeshObject, cellCellStencilObject> Stencil;
 
 class cellCellStencilObject
 :
-    public MeshObject<fvMesh, MoveableMeshObject, cellCellStencilObject>,
+    public Stencil,
     public cellCellStencil
 {
-    // Private Typedefs
-
-        typedef MeshObject
-        <
-            fvMesh,
-            MoveableMeshObject,
-            cellCellStencilObject
-        > MeshObject_type;
-
-
     // Private Data
 
         autoPtr<cellCellStencil> stencilPtr_;
@@ -86,8 +76,7 @@ public:
             const bool update = true
         )
         :
-            MeshObject_type(mesh),
-
+            Stencil(mesh),
             cellCellStencil(mesh),
             stencilPtr_
             (
@@ -111,10 +100,7 @@ public:
     // Member Functions
 
         //- Callback for geometry motion
-        virtual bool movePoints()
-        {
-            return stencilPtr_().update();
-        }
+        virtual bool movePoints();
 
         //- Update stencils. Return false if nothing changed.
         virtual bool update()

src/overset/oversetFvMesh/dynamicOversetFvMesh/dynamicOversetFvMesh.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2022,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -65,6 +65,10 @@ bool Foam::dynamicOversetFvMesh::update()
         return false;
     }
 
+    //Note: too late to do oversetFvMeshBase::clearOut() to get it
+    //      consistent with any new cell-cell stencil since 
+    //      dynamicMotionSolverListFvMesh already triggers
+    //      meshObject::movePoints on cellCellStencilObject
     oversetFvMeshBase::update();
 
     return true;

src/overset/oversetFvMesh/oversetFvMeshBase.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2014-2022 OpenCFD Ltd.
+    Copyright (C) 2014-2022,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -578,6 +578,14 @@ Foam::oversetFvMeshBase::primitiveLduAddr() const
 }
 
 
+void Foam::oversetFvMeshBase::clearOut()
+{
+    // Cell-cell stencil is already mesh object. Clear out local
+    // addressing to force rebuilding addressing
+    lduPtr_.clear();
+}
+
+
 bool Foam::oversetFvMeshBase::update()
 {
     {

src/overset/oversetFvMesh/oversetFvMeshBase.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2022,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,7 +37,6 @@ SourceFiles
 #ifndef oversetFvMeshBase_H
 #define oversetFvMeshBase_H
 
-#include "fvMeshPrimitiveLduAddressing.H"
 #include "fvMeshPrimitiveLduAddressing.H"
 #include "lduInterfaceFieldPtrsList.H"
 #include "volFieldsFwd.H"
@@ -238,6 +237,9 @@ public:
             ) const;
 
 
+        //- Clear out local storage
+        void clearOut();
+
         //- Update the mesh for both mesh motion and topology change
         virtual bool update();
 

src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C

@@ -75,6 +75,7 @@ void Foam::RBD::rigidBodySolvers::symplectic::solve
     qDot() = qDot0() + 0.5*deltaT0()*qDdot();
     q() = q0() + deltaT()*qDot();
 
+    // Update joints from new locations (q, q0)
     correctQuaternionJoints();
 
     // Update the body-state prior to the evaluation of the restraints

src/surfMesh/writers/nastran/nastranSurfaceWriter.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -307,16 +307,10 @@ void Foam::surfaceWriters::nastranWriter::writeGeometry
 Foam::surfaceWriters::nastranWriter::nastranWriter()
 :
     surfaceWriter(),
-    writeFormat_(fieldFormat::SHORT),
-    fieldMap_(),
+    writeFormat_(fieldFormat::FREE),
     commonGeometry_(false),
-    separator_()
-{
-    // if (writeFormat_ == fieldFormat::FREE)
-    // {
-    //     separator_ = ",";
-    // }
-}
+    separator_(",")  // FREE format
+{}
 
 
 Foam::surfaceWriters::nastranWriter::nastranWriter
@@ -331,29 +325,40 @@ Foam::surfaceWriters::nastranWriter::nastranWriter
         (
             "format",
             options,
-            fieldFormat::LONG
+            fieldFormat::FREE
         )
     ),
-    fieldMap_(),
-    commonGeometry_(options.getOrDefault("commonGeometry", false)),
-    separator_()
+    commonGeometry_(options.getOrDefault("commonGeometry", false))
 {
     if (writeFormat_ == fieldFormat::FREE)
     {
         separator_ = ",";
     }
 
+    // Explicit PLOAD2, PLOAD4 field specification
+    options.readIfPresent("PLOAD2", pload2_);
+    options.readIfPresent("PLOAD4", pload4_);
+
+    // Compatibility:
+    // Optional pairs of (field name => load format)
+    // Could make conditional on (pload2_.empty() && pload4_.empty())
+
     List<Pair<word>> fieldPairs;
-    options.readEntry("fields", fieldPairs);
+    options.readIfPresent("fields", fieldPairs);
 
     for (const Pair<word>& item : fieldPairs)
     {
-        // (field name => load format)
-        fieldMap_.insert
-        (
-            item.first(),
-            fileFormats::NASCore::loadFormatNames[item.second()]
-        );
+        const loadFormat format =
+            fileFormats::NASCore::loadFormatNames[item.second()];
+
+        if (format == loadFormat::PLOAD2)
+        {
+            pload2_.push_back(item.first());
+        }
+        else if (format == loadFormat::PLOAD4)
+        {
+            pload4_.push_back(item.first());
+        }
     }
 }
 
@@ -414,12 +419,12 @@ Foam::fileName Foam::surfaceWriters::nastranWriter::write()
 
     if (UPstream::master() || !parallel_)
     {
-        if (!isDir(outputFile.path()))
+        if (!Foam::isDir(outputFile.path()))
         {
-            mkDir(outputFile.path());
+            Foam::mkDir(outputFile.path());
         }
 
-        OFstream os(outputFile);
+        OFstream os(IOstreamOption::ATOMIC, outputFile);
         fileFormats::NASCore::setPrecision(os, writeFormat_);
 
         os  << "TITLE=OpenFOAM " << outputPath_.name() << " geometry" << nl

src/surfMesh/writers/nastran/nastranSurfaceWriter.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,13 +33,15 @@ Description
     The formatOptions for nastran:
     \table
         Property    | Description                           | Reqd | Default
-        fields      | Field pairs for PLOAD2/PLOAD4         | yes |
-        format      | Nastran format (short/long/free)      | no  | long
+        format      | Nastran format (short/long/free)      | no  | free
         scale       | Output geometry scaling               | no  | 1
         transform   | Output coordinate transform           | no  |
         fieldLevel  | Subtract field level before scaling   | no  | empty dict
         fieldScale  | Output field scaling                  | no  | empty dict
         commonGeometry | use separate geometry files        | no  | false
+        PLOAD2      | Field selection (words/regex) for PLOAD2 | no  |
+        PLOAD4      | Field selection (words/regex) for PLOAD4 | no  |
+        fields      | Compat: Field pairs for PLOAD2/PLOAD4 | no  |
     \endtable
 
     For example,
@@ -48,13 +50,6 @@ Description
     {
         nastran
         {
-            // OpenFOAM field name to NASTRAN load types
-            fields
-            (
-                (pMean  PLOAD2)
-                (p      PLOAD4)
-            );
-
             format  free;     // format type
 
             scale   1000;     // [m] -> [mm]
@@ -62,10 +57,27 @@ Description
             {
                "p.*"   0.01;  // [Pa] -> [mbar]
             }
+
+            // Specific NASTRAN load types
+            PLOAD2 ( pMean );
+            PLOAD4 ( "p.*" );
+
+            // old style specification
+            fields
+            (
+                (pMean  PLOAD2)
+                (p      PLOAD4)
+            );
         }
     }
     \endverbatim
 
+    Unless otherwise specified, all fields will be treated as PLOAD2
+    output. Can optionally specify PLOAD4 output using a combination
+    of \c PLOAD4 (accept) and \c PLOAD2 (deny) entries. The older \c fields
+    specification is also accepted and will be transcribed to
+    corresponding PLOAD4, PLOAD2 entries.
+
     \section Output file locations
 
     The \c rootdir normally corresponds to something like
@@ -93,7 +105,6 @@ Description
 
 Note
     Output variable scaling does not apply to integer types such as Ids.
-    Field pairs default to PLOAD2 for scalars and PLOAD4 for vectors etc.
 
 SourceFiles
     nastranSurfaceWriter.C
@@ -106,7 +117,7 @@ SourceFiles
 
 #include "surfaceWriter.H"
 #include "NASCore.H"
-#include "HashTable.H"
+#include "wordRes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -139,15 +150,18 @@ private:
         //- Field format (width and separator)
         fieldFormat writeFormat_;
 
-        //- Mapping from field name to data format enumeration
-        HashTable<loadFormat> fieldMap_;
-
         //- Use common geometry file
         bool commonGeometry_;
 
         //- Separator (used for free format)
         word separator_;
 
+        //- Explicit selection for PLOAD2 output (deselects for PLOAD4)
+        wordRes pload2_;
+
+        //- Explicit selection for PLOAD4 output
+        wordRes pload4_;
+
 
     // Private Member Functions
 
@@ -221,10 +235,10 @@ public:
 
     // Constructors
 
-        //- Default construct. Default SHORT format
+        //- Default construct. Default FREE format
         nastranWriter();
 
-        //- Construct with some output options. Default LONG format
+        //- Construct with some output options. Default FREE format
         explicit nastranWriter(const dictionary& options);
 
         //- Construct from components

src/surfMesh/writers/nastran/nastranSurfaceWriterImpl.C

@@ -103,13 +103,13 @@ Foam::Ostream& Foam::surfaceWriters::nastranWriter::writeFaceValue
     {
         case loadFormat::PLOAD2 :
         {
-            if (pTraits<Type>::nComponents == 1)
+            if (pTraits<Type>::nComponents > 1)
             {
-                writeValue(os, value);
+                writeValue(os, Foam::mag(value));
             }
             else
             {
-                writeValue(os, mag(value));
+                writeValue(os, value);
             }
 
             os << separator_;
@@ -121,6 +121,9 @@ Foam::Ostream& Foam::surfaceWriters::nastranWriter::writeFaceValue
         {
             writeValue(os, elemId);
 
+            // NOTE: these should actually be vertex values,
+            // but misused here to provide vector quantities!
+
             for (direction d = 0; d < pTraits<Type>::nComponents; ++d)
             {
                 os  << separator_;
@@ -156,43 +159,23 @@ Foam::fileName Foam::surfaceWriters::nastranWriter::writeTemplate
 
     checkOpen();
 
-    const loadFormat format
-    (
-        fieldMap_.lookup
-        (
-            fieldName,
-            // Default format
-            (
-                pTraits<Type>::nComponents == 1
-              ? loadFormat::PLOAD2
-              : loadFormat::PLOAD4
-            )
-        )
-    );
+    // Default is PLOAD2
+    loadFormat format = loadFormat::PLOAD2;
 
     if
     (
         !std::is_integral<Type>::value  // Handle 'Ids' etc silently
-     && !fieldMap_.empty()
-     && !fieldMap_.found(fieldName)
+     && !pload4_.empty()
     )
     {
-        WarningInFunction
-            << "No mapping found between field " << fieldName
-            << " and corresponding Nastran field.  Available types:"
-            << fieldMap_ << nl;
-    }
+        const wordRes::filter matcher(pload4_, pload2_);
 
-    // Emit any common warnings
-    if (format == loadFormat::PLOAD2 && pTraits<Type>::nComponents != 1)
-    {
-        WarningInFunction
-            << fileFormats::NASCore::loadFormatNames[format]
-            << " cannot be used for higher rank values"
-            << " - reverting to mag()" << endl;
+        if (matcher(fieldName))
+        {
+            format = loadFormat::PLOAD4;
+        }
     }
 
-
     // Common geometry
     // Field:  rootdir/<TIME>/<field>_surfaceName.bdf
 
@@ -239,9 +222,9 @@ Foam::fileName Foam::surfaceWriters::nastranWriter::writeTemplate
     {
         const auto& values = tfield();
 
-        if (!isDir(outputFile.path()))
+        if (!Foam::isDir(outputFile.path()))
         {
-            mkDir(outputFile.path());
+            Foam::mkDir(outputFile.path());
         }
 
         const scalar timeValue(0);
@@ -251,7 +234,7 @@ Foam::fileName Foam::surfaceWriters::nastranWriter::writeTemplate
         DynamicList<face> decompFaces;
 
 
-        OFstream os(outputFile);
+        OFstream os(IOstreamOption::ATOMIC, outputFile);
         fileFormats::NASCore::setPrecision(os, writeFormat_);
 
         os  << "TITLE=OpenFOAM " << outputFile.name()

src/topoChangerFvMesh/rawTopoChangerFvMesh/rawTopoChangerFvMesh.C

@@ -65,12 +65,6 @@ Foam::rawTopoChangerFvMesh::rawTopoChangerFvMesh
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::rawTopoChangerFvMesh::~rawTopoChangerFvMesh()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::rawTopoChangerFvMesh::update()
@@ -121,30 +115,24 @@ bool Foam::rawTopoChangerFvMesh::update()
             }
         }
 
-        const List<objectMap>& fromFaces = topoChangeMap().facesFromFacesMap();
-
-        forAll(fromFaces, i)
+        for (const auto& map : topoChangeMap().facesFromFacesMap())
         {
-            mappedFace.set(fromFaces[i].index());
+            mappedFace.set(map.index());
         }
 
-        const List<objectMap>& fromEdges = topoChangeMap().facesFromEdgesMap();
-
-        forAll(fromEdges, i)
+        for (const auto& map : topoChangeMap().facesFromEdgesMap())
         {
-            mappedFace.set(fromEdges[i].index());
+            mappedFace.set(map.index());
         }
 
-        const List<objectMap>& fromPts = topoChangeMap().facesFromPointsMap();
-
-        forAll(fromPts, i)
+        for (const auto& map : topoChangeMap().facesFromPointsMap())
         {
-            mappedFace.set(fromPts[i].index());
+            mappedFace.set(map.index());
         }
 
         // Set unmapped faces to zero
-        Info<< "rawTopoChangerFvMesh : zeroing unmapped boundary values."
-            << endl;
+        Info<< "rawTopoChangerFvMesh : zeroing unmapped boundary values." << nl;
+
         zeroUnmappedValues<scalar, fvPatchField, volMesh>(mappedFace);
         zeroUnmappedValues<vector, fvPatchField, volMesh>(mappedFace);
         zeroUnmappedValues<sphericalTensor, fvPatchField, volMesh>(mappedFace);
@@ -155,8 +143,8 @@ bool Foam::rawTopoChangerFvMesh::update()
         Info<< "rawTopoChangerFvMesh :"
             << " recreating phi for unmapped boundary values." << endl;
 
-        const volVectorField& U = lookupObject<volVectorField>("U");
-        surfaceScalarField& phi = lookupObjectRef<surfaceScalarField>("phi");
+        const auto& U = lookupObject<volVectorField>("U");
+        auto& phi = lookupObjectRef<surfaceScalarField>("phi");
 
         setUnmappedValues
         (

src/topoChangerFvMesh/rawTopoChangerFvMesh/rawTopoChangerFvMesh.H

@@ -38,8 +38,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef rawTopoChangerFvMesh_H
-#define rawTopoChangerFvMesh_H
+#ifndef Foam_rawTopoChangerFvMesh_H
+#define Foam_rawTopoChangerFvMesh_H
 
 #include "topoChangerFvMesh.H"
 #include "bitSet.H"
@@ -49,8 +49,6 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
-
 /*---------------------------------------------------------------------------*\
                       Class rawTopoChangerFvMesh Declaration
 \*---------------------------------------------------------------------------*/
@@ -94,8 +92,9 @@ public:
             const bool doInit=true
         );
 
+
     //- Destructor
-    virtual ~rawTopoChangerFvMesh();
+    virtual ~rawTopoChangerFvMesh() = default;
 
 
     // Member Functions

src/topoChangerFvMesh/rawTopoChangerFvMesh/rawTopoChangerFvMeshTemplates.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,9 +26,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "rawTopoChangerFvMesh.H"
-#include "Time.H"
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type, template<class> class PatchField, class GeoMesh>
@@ -42,10 +40,7 @@ void Foam::rawTopoChangerFvMesh::setUnmappedValues
 
     forAll(fld.boundaryField(), patchi)
     {
-        PatchField<Type>& fvp = const_cast<PatchField<Type>&>
-        (
-            fld.boundaryField()[patchi]
-        );
+        auto& fvp = const_cast<PatchField<Type>&>(fld.boundaryField()[patchi]);
 
         const label start = fvp.patch().start();
         forAll(fvp, i)
@@ -71,33 +66,27 @@ void Foam::rawTopoChangerFvMesh::zeroUnmappedValues
 {
     typedef GeometricField<Type, PatchField, GeoMesh> FieldType;
 
-    const wordList fldNames(names(FieldType::typeName));
+    std::unique_ptr<FieldType> zeroFieldPtr;
 
-    forAll(fldNames, i)
+    for (const word& fldName : names<FieldType>())
     {
-        //Pout<< "Checking field " << fldNames[i] << endl;
+        FieldType& fld = lookupObjectRef<FieldType>(fldName);
+        //Pout<< "Checking field " << fld.name() << endl;
 
-        FieldType& fld = lookupObjectRef<FieldType>(fldNames[i]);
-
-        setUnmappedValues
-        (
-            fld,
-            mappedFace,
-            FieldType
+        if (!zeroFieldPtr)
+        {
+            zeroFieldPtr = std::make_unique<FieldType>
             (
-                IOobject
-                (
-                    "zero",
-                    time().timeName(),
-                    *this,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    IOobject::NO_REGISTER
-                ),
+                this->newIOobject("zero"),
                 *this,
-                dimensioned<Type>(fld.dimensions(), Zero)
-            )
-        );
+                Foam::zero{},
+                dimless
+            );
+        }
+
+        zeroFieldPtr->dimensions().reset(fld.dimensions());
+
+        setUnmappedValues(fld, mappedFace, *zeroFieldPtr);
     }
 }
 

tutorials/mesh/blockMesh/sphere/system/blockMeshDict

@@ -35,7 +35,7 @@ geometry
 }
 
 // Box size
-vo   #eval{ sqrt($outerRadius/3) };
+vo   #eval{ $outerRadius/sqrt(3) };
 
 vertices
 (

tutorials/mesh/blockMesh/sphere7/system/blockMeshDict

@@ -38,7 +38,7 @@ geometry
 }
 
 // Outer box size
-vo   #eval{ sqrt($outerRadius/3) };
+vo   #eval{ $outerRadius/sqrt(3) };
 
 // Inner box size - % of overall dimension
 vi   #eval{ $vo * $innerRatio };

tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/blockMeshDict

@@ -45,7 +45,7 @@ geometry
 
 
 // Outer box size (approximate)
-vo   #eval{ sqrt($outerRadius/3) };
+vo   #eval{ $outerRadius/sqrt(3) };
 
 // Inner box size - % of overall dimension
 vi   #eval{ $vo * $innerRatio };