Code adjustments to avoid gcc-13 warnings about dangling references

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C

@@ -69,7 +69,8 @@ Foam::XiEqModels::Gulder::~Gulder()
 Foam::tmp<Foam::volScalarField> Foam::XiEqModels::Gulder::XiEq() const
 {
     volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
-    const volScalarField& epsilon = turbulence_.epsilon();
+    const tmp<volScalarField> tepsilon(turbulence_.epsilon());
+    const volScalarField& epsilon = tepsilon();
 
     if (subGridSchelkin_)
     {

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C

@@ -75,8 +75,10 @@ Foam::XiEqModels::SCOPEXiEq::~SCOPEXiEq()
 
 Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
 {
-    const volScalarField& k = turbulence_.k();
+    const tmp<volScalarField> tk(turbulence_.k());
-    const volScalarField& epsilon = turbulence_.epsilon();
+    const volScalarField& k = tk();
+    const tmp<volScalarField> tepsilon(turbulence_.epsilon());
+    const volScalarField& epsilon = tepsilon();
 
     volScalarField up(sqrt((2.0/3.0)*k));
     if (subGridSchelkin_)

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/KTS/KTS.C

@@ -66,7 +66,8 @@ Foam::XiGModels::KTS::~KTS()
 Foam::tmp<Foam::volScalarField> Foam::XiGModels::KTS::G() const
 {
     volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
-    const volScalarField& epsilon = turbulence_.epsilon();
+    const tmp<volScalarField> tepsilon(turbulence_.epsilon());
+    const volScalarField& epsilon = tepsilon();
 
     volScalarField tauEta(sqrt(mag(thermo_.muu()/(thermo_.rhou()*epsilon))));
 

applications/solvers/compressible/rhoCentralFoam/createFieldRefs.H

@@ -1,10 +1,9 @@
 volScalarField& p = thermo.p();
 const volScalarField& T = thermo.T();
 const volScalarField& psi = thermo.psi();
-const volScalarField& mu = thermo.mu();
 
 bool inviscid(true);
-if (max(mu.primitiveField()) > 0.0)
+if (max(thermo.mu().cref().primitiveField()) > 0.0)
 {
     inviscid = false;
 }

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/VoFPatchTransfer/VoFPatchTransfer.C

@@ -160,7 +160,8 @@ void VoFPatchTransfer::correct
     const volScalarField& heVoF = thermo.thermo1().he();
     const volScalarField& TVoF = thermo.thermo1().T();
     const volScalarField CpVoF(thermo.thermo1().Cp());
-    const volScalarField& rhoVoF = thermo.thermo1().rho()();
+    const tmp<volScalarField> trhoVoF(thermo.thermo1().rho());
+    const volScalarField& rhoVoF = trhoVoF();
     const volScalarField& alphaVoF = thermo.alpha1();
 
     forAll(patchIDs_, pidi)

applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H

@@ -28,7 +28,8 @@
     forAllConstIters(mixture.phases(), phase)
     {
         const rhoThermo& thermo = phase().thermo();
-        const volScalarField& rho = thermo.rho()();
+        const tmp<volScalarField> trho(thermo.rho());
+        const volScalarField& rho = trho();
 
         p_rghEqnComps.set
         (

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C

@@ -482,8 +482,8 @@ void Foam::searchableSurfaceControl::cellSizeFunctionVertices
     DynamicList<scalar>& sizes
 ) const
 {
-    const tmp<pointField> tmpPoints = searchableSurface_.points();
+    const tmp<pointField> tpoints(searchableSurface_.points());
-    const pointField& points = tmpPoints();
+    const pointField& points = tpoints();
 
     const scalar nearFeatDistSqrCoeff = 1e-8;
 

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/nonUniformField/nonUniformField.C

@@ -110,7 +110,8 @@ Foam::scalar Foam::nonUniformField::interpolate
 {
     const face& faceHitByPt = surfaceTriMesh_.triSurface::operator[](index);
 
-    const pointField& pts = surfaceTriMesh_.points();
+    const tmp<pointField> tpoints(surfaceTriMesh_.points());
+    const pointField& pts = tpoints();
 //    const Map<label>& pMap = surfaceTriMesh_.meshPointMap();
 
     triPointRef tri

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C

@@ -982,7 +982,7 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::findOffsetPatchFaces
         offsetBoundaryCells.write();
     }
 
-    return std::move(offsetBoundaryCells);
+    return offsetBoundaryCells;
 }
 
 

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C

@@ -1371,7 +1371,7 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::findRemainingProtrusionSet
         protrudingCells.write();
     }
 
-    return std::move(protrudingCells);
+    return protrudingCells;
 }
 
 

applications/utilities/surface/surfaceInflate/surfaceInflate.C

@@ -221,9 +221,10 @@ void detectSelfIntersections
     const edgeList& edges = s.edges();
     const indexedOctree<treeDataTriSurface>& tree = s.tree();
     const labelList& meshPoints = s.meshPoints();
-    const pointField& points = s.points();
+    const tmp<pointField> tpoints(s.points());
+    const pointField& points = tpoints();
 
-    isEdgeIntersecting.setSize(edges.size());
+    isEdgeIntersecting.resize_nocopy(edges.size());
     isEdgeIntersecting = false;
 
     forAll(edges, edgeI)
@@ -311,7 +312,8 @@ label detectIntersectionPoints
         detectSelfIntersections(s, isEdgeIntersecting);
 
         const edgeList& edges = s.edges();
-        const pointField& points = s.points();
+        const tmp<pointField> tpoints(s.points());
+        const pointField& points = tpoints();
 
         forAll(edges, edgeI)
         {
@@ -836,9 +838,10 @@ int main(int argc, char *argv[])
         // Do some smoothing (Lloyds algorithm)
         lloydsSmoothing(nSmooth, s, isFeaturePoint, edgeStat, isAffectedPoint);
 
-
         // Update pointDisplacement
-        const pointField& pts = s.points();
+        const tmp<pointField> tpoints(s.points());
+        const pointField& pts = tpoints();
+
         forAll(meshPoints, i)
         {
             label meshPointI = meshPoints[i];

src/OpenFOAM/fields/cloud/cloud.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -172,6 +172,21 @@ public:
             {
                 return obr.lookupObject<IOField<Type>>(fieldName);
             }
+
+            //- Lookup an IOField within object registry
+            //  Fatal if not found or wrong type
+            //
+            //  Note: const char signature to avoid spurious
+            //  -Wdangling-reference with gcc-13
+            template<class Type>
+            static const IOField<Type>& lookupIOField
+            (
+                const char* fieldName,
+                const objectRegistry& obr
+            )
+            {
+                return obr.lookupObject<IOField<Type>>(word(fieldName));
+            }
 };
 
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C

@@ -289,10 +289,12 @@ void Foam::processorPolyPatch::calcGeometry(PstreamBuffers& pBufs)
 
                 label vertI = 0;
 
-                forAll(faceCentres(), facej)
+                const vectorField::subField ownFaceCentres = faceCentres();
+
+                forAll(ownFaceCentres, facej)
                 {
                     const point& c0 = neighbFaceCentres_[facej];
-                    const point& c1 = faceCentres()[facej];
+                    const point& c1 = ownFaceCentres[facej];
 
                     writeOBJ(ccStr, c0, c1, vertI);
                 }

src/atmosphericModels/fvOptions/atmAmbientTurbSource/atmAmbientTurbSourceTemplates.C

@@ -45,7 +45,8 @@ void Foam::fv::atmAmbientTurbSource::atmAmbientTurbSourceEpsilon
         (
             turbulenceModel::propertiesName
         );
-    const volScalarField& epsilon = turbPtr->epsilon();
+    const tmp<volScalarField> tepsilon(turbPtr->epsilon());
+    const volScalarField& epsilon = tepsilon();
 
     // (Heuristically derived from RS:Eq. 4, rhs-term:5)
     eqn +=
@@ -67,7 +68,9 @@ void Foam::fv::atmAmbientTurbSource::atmAmbientTurbSourceOmega
         (
             turbulenceModel::propertiesName
         );
-    const volScalarField& omega = turbPtr->omega();
+    const tmp<volScalarField> tomega(turbPtr->omega());
+    const volScalarField& omega = tomega();
+
     const volScalarField::Internal& beta =
         mesh_.lookupObjectRef<volScalarField::Internal>
         (
@@ -93,7 +96,8 @@ void Foam::fv::atmAmbientTurbSource::atmAmbientTurbSourceK
         (
             turbulenceModel::propertiesName
         );
-    const volScalarField& k = turbPtr->k();
+    const tmp<volScalarField> tk(turbPtr->k());
+    const volScalarField& k = tk();
 
     if (isEpsilon_)
     {

src/atmosphericModels/fvOptions/atmBuoyancyTurbSource/atmBuoyancyTurbSourceTemplates.C

@@ -47,8 +47,10 @@ void Foam::fv::atmBuoyancyTurbSource::atmBuoyancyTurbSourceEpsilon
         );
 
     // Fetch required fields from the epsilon-based model
-    const volScalarField& k = turbPtr->k();
+    const tmp<volScalarField> tk(turbPtr->k());
-    const volScalarField& epsilon = turbPtr->epsilon();
+    const volScalarField& k = tk();
+    const tmp<volScalarField> tepsilon(turbPtr->epsilon());
+    const volScalarField& epsilon = tepsilon();
     const volScalarField::Internal& GbyNu =
         mesh_.lookupObjectRef<volScalarField::Internal>
         (
@@ -77,8 +79,10 @@ void Foam::fv::atmBuoyancyTurbSource::atmBuoyancyTurbSourceOmega
         );
 
     // Fetch required fields from the omega-based model
-    const volScalarField& k = turbPtr->k();
+    const tmp<volScalarField> tk(turbPtr->k());
-    const volScalarField& omega = turbPtr->omega();
+    const volScalarField& k = tk();
+    const tmp<volScalarField> tomega(turbPtr->omega());
+    const volScalarField& omega = tomega();
     const volScalarField::Internal& GbyNu =
         mesh_.lookupObjectRef<volScalarField::Internal>
         (
@@ -121,7 +125,8 @@ void Foam::fv::atmBuoyancyTurbSource::atmBuoyancyTurbSourceK
             turbulenceModel::propertiesName
         );
 
-    const volScalarField& k = turbPtr->k();
+    const tmp<volScalarField> tk(turbPtr->k());
+    const volScalarField& k = tk();
 
     eqn += fvm::Sp(alpha()*rho()*B_/k(), k);
 }

src/atmosphericModels/fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSourceTemplates.C

@@ -45,7 +45,8 @@ void Foam::fv::atmPlantCanopyTurbSource::atmPlantCanopyTurbSourceEpsilon
         (
             turbulenceModel::propertiesName
         );
-    const volScalarField& epsilon = turbPtr->epsilon();
+    const tmp<volScalarField> tepsilon(turbPtr->epsilon());
+    const volScalarField& epsilon = tepsilon();
     const volVectorField::Internal& U = turbPtr->U()();
 
     eqn -= fvm::Sp(alpha()*rho()*(C1_ - C2_)*calcPlantCanopyTerm(U), epsilon);
@@ -66,7 +67,8 @@ void Foam::fv::atmPlantCanopyTurbSource::atmPlantCanopyTurbSourceOmega
         (
             turbulenceModel::propertiesName
         );
-    const volScalarField& omega = turbPtr->omega();
+    const tmp<volScalarField> tomega(turbPtr->omega());
+    const volScalarField& omega = tomega();
     const volVectorField::Internal& U = turbPtr->U()();
     const volScalarField::Internal& gamma =
         mesh_.lookupObjectRef<volScalarField::Internal>

src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.H

@@ -38,8 +38,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef singleCellFvMesh_H
+#ifndef Foam_singleCellFvMesh_H
-#define singleCellFvMesh_H
+#define Foam_singleCellFvMesh_H
 
 #include "fvPatchFieldMapper.H"
 #include "fvMesh.H"
@@ -58,7 +58,7 @@ class singleCellFvMesh
 :
     public fvMesh
 {
-    // Private data
+    // Private Data
 
         const labelListIOList patchFaceAgglomeration_;
 
@@ -179,36 +179,39 @@ public:
 
     // Member Functions
 
-        bool agglomerate() const
+        bool agglomerate() const noexcept
         {
-            return patchFaceAgglomeration_.size() > 0;
+            return !patchFaceAgglomeration_.empty();
         }
 
         //- From patchFace on this back to original mesh or agglomeration
-        const labelListList& patchFaceMap() const
+        const labelListList& patchFaceMap() const noexcept
         {
             return patchFaceMap_;
         }
 
         //- From point on this back to original mesh
-        const labelList& pointMap() const
+        const labelList& pointMap() const noexcept
         {
             return pointMap_;
         }
 
         //- From face on original mesh to face on this
-        const labelList& reverseFaceMap() const
+        const labelList& reverseFaceMap() const noexcept
         {
             return reverseFaceMap_;
         }
 
         //- From point on original mesh to point on this (or -1 for removed
         //- points)
-        const labelList& reversePointMap() const
+        const labelList& reversePointMap() const noexcept
         {
             return reversePointMap_;
         }
 
+        //- Interpolate for overset (unused)
+        using fvMesh::interpolate;
+
         //- Map volField. Internal field set to average, patch fields straight
         //- copies.
         template<class Type>
@@ -217,7 +220,6 @@ public:
         (
             const GeometricField<Type, fvPatchField, volMesh>&
         ) const;
-
 };
 
 

src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.C

@@ -132,9 +132,9 @@ void Foam::patchTransformedInterpolation::interpolate
     labelList pointDisplacementNSum(nPoints, Zero);
     vectorField pointDisplacementSum(nPoints, Zero);
 
-    forAll(patches_, patchI)
+    for (const label patchi : patches_)
     {
-        const polyPatch& patch(mesh().boundaryMesh()[patches_[patchI]]);
+        const polyPatch& patch = mesh().boundaryMesh()[patchi];
 
         forAll(patch, pFaceI)
         {
@@ -145,7 +145,7 @@ void Foam::patchTransformedInterpolation::interpolate
             const labelList cPoints(c.labels(mesh().faces()));
 
             // Consider movement around the face centre
-            const point& xOrigin(patch.faceCentres()[pFaceI]);
+            const point xOrigin(patch.faceCentres()[pFaceI]);
 
             // Mean translation
             const vector uMean(f.average(points, pointDisplacement));

src/fvOptions/sources/derived/buoyancyTurbSource/buoyancyTurbSourceTemplates.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd
+    Copyright (C) 2020,2023 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,7 +46,8 @@ void Foam::fv::buoyancyTurbSource::buoyancyTurbSourceK
         (
             turbulenceModel::propertiesName
         );
-    const volScalarField& nut = turbPtr->nut();
+    const tmp<volScalarField> tnut(turbPtr->nut());
+    const volScalarField& nut = tnut();
 
     const dictionary& turbDict = turbPtr->coeffDict();
     const scalar Prt

src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C

@@ -231,7 +231,7 @@ void Foam::FreeStream<CloudType>::inflow()
 
             label celli = mesh.faceOwner()[globalFaceIndex];
 
-            const vector& fC = patch.faceCentres()[pFI];
+            const vector fC = patch.faceCentres()[pFI];
 
             scalar fA = mag(patch.faceAreas()[pFI]);
 

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C

@@ -3805,26 +3805,27 @@ const Foam::dictionary& Foam::meshRefinement::subDict
     enum keyType::option matchOpt
 )
 {
-    const auto finder(dict.csearch(keyword, matchOpt));
+    const dictionary* dictptr = dict.findDict(keyword, matchOpt);
 
-    if (!finder.good())
+    if (!dictptr)
     {
-        auto& err = FatalIOErrorInFunction(dict);
+        FatalIOErrorInFunction(dict)
-
+            << "Entry '" << keyword
-        err << "Entry '" << keyword << "' not found in dictionary "
+            << "' not found (or not a dictionary) in dictionary "
             << dict.relativeName() << nl;
 
         if (noExit)
         {
+            // Dummy return
             return dictionary::null;
         }
         else
         {
-            err << exit(FatalIOError);
+            FatalIOError << exit(FatalIOError);
         }
     }
 
-    return finder.dict();
+    return *dictptr;
 }
 
 
@@ -3840,19 +3841,18 @@ Foam::ITstream& Foam::meshRefinement::lookup
 
     if (!eptr)
     {
-        auto& err = FatalIOErrorInFunction(dict);
+        FatalIOErrorInFunction(dict)
-
+            << "Entry '" << keyword << "' not found in dictionary "
-        err << "Entry '" << keyword << "' not found in dictionary "
             << dict.relativeName() << nl;
 
         if (noExit)
         {
-            // Fake entry
+            // Dummy return
-            return dict.first()->stream();
+            return ITstream::empty_stream();
         }
         else
         {
-            err << exit(FatalIOError);
+            FatalIOError << exit(FatalIOError);
         }
     }
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/adjointSensitivity/adjointEikonalSolver/adjointEikonalSolver.C

@@ -83,7 +83,8 @@ void adjointEikonalSolver::read()
 tmp<surfaceScalarField> adjointEikonalSolver::computeYPhi()
 {
     // Primal distance field
-    const volScalarField& d = adjointSolver_.yWall();
+    const tmp<volScalarField> td(adjointSolver_.yWall());
+    const volScalarField& d = td();
 
     volVectorField ny
     (
@@ -198,7 +199,8 @@ void adjointEikonalSolver::solve()
     read();
 
     // Primal distance field
-    const volScalarField& d = adjointSolver_.yWall();
+    const tmp<volScalarField> td(adjointSolver_.yWall());
+    const volScalarField& d = td();
 
     // Convecting flux
     tmp<surfaceScalarField> tyPhi = computeYPhi();
@@ -275,7 +277,9 @@ boundaryVectorField& adjointEikonalSolver::distanceSensitivities()
 
     boundaryVectorField& distanceSens = distanceSensPtr_();
 
-    const volScalarField& d = adjointSolver_.yWall();
+    const tmp<volScalarField> td(adjointSolver_.yWall());
+    const volScalarField& d = td();
+
     for (const label patchi : sensitivityPatchIDs_)
     {
         vectorField nf(mesh_.boundary()[patchi].nf());
@@ -294,7 +298,8 @@ tmp<volTensorField> adjointEikonalSolver::getFISensitivityTerm() const
 {
     Info<< "Calculating distance sensitivities " << endl;
 
-    const volScalarField& d = adjointSolver_.yWall();
+    const tmp<volScalarField> td(adjointSolver_.yWall());
+    const volScalarField& d = td();
     const volVectorField gradD(fvc::grad(d));
 
     auto gradDDa
@@ -354,7 +359,8 @@ tmp<scalarField> adjointEikonalSolver::topologySensitivities
     const word& designVarsName
 ) const
 {
-    const volScalarField& d = adjointSolver_.yWall();
+    const tmp<volScalarField> td(adjointSolver_.yWall());
+    const volScalarField& d = td();
 
     auto tres(tmp<scalarField>::New(d.primitiveField().size(), Zero));
     scalarField dSens(d.primitiveField()*da_.primitiveField());
@@ -377,7 +383,9 @@ const volScalarField& adjointEikonalSolver::da()
 
 tmp<volVectorField> adjointEikonalSolver::gradEikonal()
 {
-    const volScalarField& d = adjointSolver_.yWall();
+    const tmp<volScalarField> td(adjointSolver_.yWall());
+    const volScalarField& d = td();
+
     volVectorField gradD(fvc::grad(d));
     return tmp<volVectorField>::New("gradEikonal", 2*gradD & fvc::grad(gradD));
 }

src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/SpalartAllmaras/SpalartAllmaras.C

@@ -95,7 +95,8 @@ tmp<volScalarField> SpalartAllmaras::nutJacobianVar1
     );
     auto& nutJacobian = tnutJacobian.ref();
 
-    const volScalarField& nu = laminarTransport.nu();
+    const tmp<volScalarField> tnu(laminarTransport.nu());
+    const volScalarField& nu = tnu();
     const volScalarField& nuTilda = TMVar1();
 
     volScalarField chi(nuTilda/nu);

src/phaseSystemModels/multiphaseEuler/multiphaseSystem/phaseModel/phaseModel.C

@@ -241,7 +241,8 @@ void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const
 {
     surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
     const volScalarField::Boundary& alphaBf = boundaryField();
-    const surfaceScalarField::Boundary& phiBf = phi().boundaryField();
+    const tmp<surfaceScalarField> tphi(phi());
+    const auto& phiBf = tphi().boundaryField();
 
     forAll(alphaPhiBf, patchi)
     {

src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C

@@ -355,7 +355,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer()
                 fvVectorMatrix& eqn = *eqns[phase.name()];
 
                 const volVectorField& U = eqn.psi();
-                const surfaceScalarField& phi = phase.phi();
+                const tmp<surfaceScalarField> tphi(phase.phi());
+                const surfaceScalarField& phi = tphi();
 
                 eqn -=
                     Vm
@@ -408,7 +409,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransferf()
 
         if (!phase.stationary())
         {
-            const volVectorField& U = phase.U();
+            const tmp<volVectorField> tU(phase.U());
+            const volVectorField& U = tU();
 
             UgradUs.set
             (

src/phaseSystemModels/reactingEuler/multiphaseSystem/interfacialModels/dragModels/segregated/segregated.C

@@ -85,8 +85,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
     const volScalarField& alpha1(pair_.phase1());
     const volScalarField& alpha2(pair_.phase2());
 
-    const volScalarField& rho1(pair_.phase1().rho());
+    const tmp<volScalarField> trho1(pair_.phase1().rho());
-    const volScalarField& rho2(pair_.phase2().rho());
+    const tmp<volScalarField> trho2(pair_.phase2().rho());
+
+    const volScalarField& rho1 = trho1();
+    const volScalarField& rho2 = trho2();
 
     tmp<volScalarField> tnu1(pair_.phase1().nu());
     tmp<volScalarField> tnu2(pair_.phase2().nu());

src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C

@@ -93,7 +93,8 @@ Foam::tmp<Foam::fvScalarMatrix>
 Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
 {
     const volScalarField& alpha = *this;
-    const volScalarField& rho = this->rho();
+    const tmp<volScalarField> trho(this->rho());
+    const volScalarField& rho(trho());
 
     const tmp<volVectorField> tU(this->U());
     const volVectorField& U(tU());

src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MovingPhaseModel/MovingPhaseModel.C

@@ -283,7 +283,8 @@ Foam::tmp<Foam::fvVectorMatrix>
 Foam::MovingPhaseModel<BasePhaseModel>::UEqn()
 {
     const volScalarField& alpha = *this;
-    const volScalarField& rho = this->thermo().rho();
+    const tmp<volScalarField> trho = this->thermo().rho();
+    const volScalarField& rho = trho();
 
     return
     (
@@ -303,14 +304,13 @@ Foam::MovingPhaseModel<BasePhaseModel>::UfEqn()
     // As the "normal" U-eqn but without the ddt terms
 
     const volScalarField& alpha = *this;
-    const volScalarField& rho = this->thermo().rho();
 
     return
     (
         fvm::div(alphaRhoPhi_, U_)
       - fvm::Sp(fvc::div(alphaRhoPhi_), U_)
       + fvm::SuSp(- this->continuityErrorSources(), U_)
-      + this->fluid().MRF().DDt(alpha*rho, U_)
+      + this->fluid().MRF().DDt(alpha*this->thermo().rho(), U_)
       + turbulence_->divDevRhoReff(U_)
     );
 }

src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C

@@ -152,7 +152,8 @@ Foam::MultiComponentPhaseModel<BasePhaseModel>::YiEqn(volScalarField& Yi)
 {
     const volScalarField& alpha = *this;
     const surfaceScalarField alphaRhoPhi(this->alphaRhoPhi());
-    const volScalarField& rho = this->thermo().rho();
+    const tmp<volScalarField> trho(this->thermo().rho());
+    const volScalarField& rho = trho();
 
     return
     (

src/phaseSystemModels/reactingEuler/multiphaseSystem/phaseModel/phaseModel/phaseModel.C

@@ -207,7 +207,8 @@ void Foam::phaseModel::correctInflowOutflow(surfaceScalarField& alphaPhi) const
 {
     surfaceScalarField::Boundary& alphaPhiBf = alphaPhi.boundaryFieldRef();
     const volScalarField::Boundary& alphaBf = boundaryField();
-    const surfaceScalarField::Boundary& phiBf = phi()().boundaryField();
+    const tmp<surfaceScalarField> tphi(phi());
+    const auto& phiBf = tphi().boundaryField();
 
     forAll(alphaPhiBf, patchi)
     {

src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/nucleationModels/wallBoiling/wallBoiling.C

@@ -154,7 +154,8 @@ Foam::diameterModels::nucleationModels::wallBoiling::addToNucleationRate
 {
     const sizeGroup& fi = popBal_.sizeGroups()[i];
     const phaseModel& phase = fi.phase();
-    const volScalarField& rho = phase.rho();
+    const tmp<volScalarField> trho(phase.rho());
+    const volScalarField& rho = trho();
     const tmp<volScalarField> talphat(turbulence_.alphat());
     const volScalarField::Boundary& alphatBf = talphat().boundaryField();
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/populationBalanceModel/populationBalanceModel/populationBalanceModel.C

@@ -1273,7 +1273,8 @@ void Foam::diameterModels::populationBalanceModel::solve()
                 const phaseModel& phase = fi.phase();
                 const volScalarField& alpha = phase;
                 const dimensionedScalar& residualAlpha = phase.residualAlpha();
-                const volScalarField& rho = phase.thermo().rho();
+                const tmp<volScalarField> trho(phase.thermo().rho());
+                const volScalarField& rho = trho();
 
                 fvScalarMatrix sizeGroupEqn
                 (

src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJacksonSchaeffer/JohnsonJacksonSchaefferFrictionalStress.C

@@ -156,7 +156,8 @@ JohnsonJacksonSchaeffer::nu
     }
 
     const fvPatchList& patches = phase.mesh().boundary();
-    const volVectorField& U = phase.U();
+    const tmp<volVectorField> tU(phase.U());
+    const volVectorField& U = tU();
 
     volScalarField::Boundary& nufBf = nuf.boundaryFieldRef();
 

src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C

@@ -154,7 +154,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::nu
     }
 
     const fvPatchList& patches = phase.mesh().boundary();
-    const volVectorField& U = phase.U();
+    const tmp<volVectorField> tU(phase.U());
+    const volVectorField& U = tU();
 
     volScalarField::Boundary& nufBf = nuf.boundaryFieldRef();
 

src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C

@@ -266,7 +266,8 @@ Foam::RASModels::kineticTheoryModel::R() const
 Foam::tmp<Foam::volScalarField>
 Foam::RASModels::kineticTheoryModel::pPrime() const
 {
-    const volScalarField& rho = phase_.rho();
+    const tmp<volScalarField> trho(phase_.rho());
+    const volScalarField& rho = trho();
 
     tmp<volScalarField> tpPrime
     (
@@ -365,11 +366,13 @@ void Foam::RASModels::kineticTheoryModel::correct()
 {
     // Local references
     volScalarField alpha(max(alpha_, scalar(0)));
-    const volScalarField& rho = phase_.rho();
+    const tmp<volScalarField> trho(phase_.rho());
+    const volScalarField& rho = trho();
     const surfaceScalarField& alphaRhoPhi = alphaRhoPhi_;
     const volVectorField& U = U_;
-    const volVectorField& Uc_ =
+    const tmp<volVectorField> tUc =
         refCast<const twoPhaseSystem>(phase_.fluid()).otherPhase(phase_).U();
+    const volVectorField& Uc = tUc();
 
     const scalar sqrtPi = sqrt(constant::mathematical::pi);
     dimensionedScalar ThetaSmall("ThetaSmall", Theta_.dimensions(), 1e-6);
@@ -421,7 +424,7 @@ void Foam::RASModels::kineticTheoryModel::correct()
         volScalarField J2
         (
             "J2",
-            0.25*sqr(beta)*da*magSqr(U - Uc_)
+            0.25*sqr(beta)*da*magSqr(U - Uc)
            /(
                max(alpha, residualAlpha_)*rho
               *sqrtPi*(ThetaSqrt + ThetaSmallSqrt)

src/phaseSystemModels/reactingEuler/twoPhaseSystem/twoPhaseSystem.C

@@ -136,8 +136,10 @@ void Foam::twoPhaseSystem::solve()
     word alphaScheme("div(phi," + alpha1.name() + ')');
     word alpharScheme("div(phir," + alpha1.name() + ')');
 
-    const surfaceScalarField& phi1 = phase1_.phi();
+    const tmp<surfaceScalarField> tphi1(phase1_.phi());
-    const surfaceScalarField& phi2 = phase2_.phi();
+    const surfaceScalarField& phi1 = tphi1();
+    const tmp<surfaceScalarField> tphi2(phase2_.phi());
+    const surfaceScalarField& phi2 = tphi2();
 
     // Construct the dilatation rate source term
     tmp<volScalarField::Internal> tdgdt;

src/phaseSystemModels/twoPhaseEuler/phaseCompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C

@@ -154,7 +154,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::nu
     }
 
     const fvPatchList& patches = phase.mesh().boundary();
-    const volVectorField& U = phase.U();
+    const tmp<volVectorField> tU(phase.U());
+    const volVectorField& U = tU();
 
     volScalarField::Boundary& nufBf = nuf.boundaryFieldRef();
 

src/phaseSystemModels/twoPhaseEuler/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C

@@ -80,11 +80,12 @@ Foam::diameterModels::IATEsources::turbulentBreakUp::R() const
 
     volScalarField R = tR();
 
-    scalar Cti = Cti_.value();
+    const scalar Cti = Cti_.value();
-    scalar WeCr = WeCr_.value();
+    const scalar WeCr = WeCr_.value();
     volScalarField Ut(this->Ut());
     volScalarField We(this->We());
-    const volScalarField& d(iate_.d()());
+    const tmp<volScalarField> td(iate_.d());
+    const volScalarField& d = td();
 
     forAll(R, celli)
     {

src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.C

@@ -192,7 +192,7 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::updateCoeffs()
     }
 
     const scalarField UAvePyr(-phiPyr/patch().magSf());
-    const vectorField& nf = patch().nf();
+    tmp<vectorField> nf(patch().nf());
 
 
     // Evaluate velocity

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/multiBandZoneAbsorptionEmission/multiBandZoneAbsorptionEmission.C

@@ -61,8 +61,8 @@ multiBandZoneAbsorptionEmission
     emiCoeffs_(maxBands_),
     nBands_(0),
     zoneAbsorptivity_(),
-    zoneEmisivity_(),
+    zoneEmissivity_(),
-    zoneCells_()
+    zoneIds_()
 {
     coeffsDict_.readEntry("absorptivity", absCoeffs_);
     coeffsDict_.readEntry("emissivity", emiCoeffs_);
@@ -71,11 +71,11 @@ multiBandZoneAbsorptionEmission
     const dictionary& zoneDict = coeffsDict_.subDict("zones");
 
     zoneDict.readEntry("absorptivity", zoneAbsorptivity_);
-    zoneDict.readEntry("emissivity", zoneEmisivity_);
+    zoneDict.readEntry("emissivity", zoneEmissivity_);
 
-    zoneCells_.setSize(zoneAbsorptivity_.size(), -1);
+    zoneIds_.resize(zoneAbsorptivity_.size(), -1);
 
-    label i = 0;
+    label numZones = 0;
     forAllConstIters(zoneAbsorptivity_, iter)
     {
         label zoneID = mesh.cellZones().findZoneID(iter.key());
@@ -86,9 +86,10 @@ multiBandZoneAbsorptionEmission
                 << "Valid cellZones are " << mesh.cellZones().names()
                 << exit(FatalError);
         }
-        zoneCells_[i++] = zoneID;
+        zoneIds_[numZones] = zoneID;
+        ++numZones;
     }
-
+    // zoneIds_.resize(numZones);
 }
 
 
@@ -124,24 +125,19 @@ Foam::radiation::multiBandZoneAbsorptionEmission::aCont
         )
     );
 
-    volScalarField& a = ta.ref();
+    scalarField& a = ta.ref().primitiveFieldRef();
 
-    forAll(zoneCells_, zonei)
+    for (const label zonei : zoneIds_)
     {
-        const cellZone& cZone = mesh().cellZones()[zoneCells_[zonei]];
+        const cellZone& zn = mesh().cellZones()[zonei];
+        const auto iter = zoneAbsorptivity_.cfind(zn.name());
 
-        tmp<volScalarField> tzoneAbs(a*0.0);
+        if (iter.good())  // Check is redundant (cannot fail)
-        volScalarField& zoneAbs = tzoneAbs.ref();
-
-        const scalarList& abs = zoneAbsorptivity_.find(cZone.name())();
-
-        forAll(cZone, i)
         {
-            label cellId = cZone[i];
+            const scalarList& absorb = iter.val();
-            zoneAbs[cellId] = abs[bandI] - absCoeffs_[bandI];
+
+            UIndirectList<scalar>(a, zn) = absorb[bandI];
         }
-
-        a += zoneAbs;
     }
 
     return ta;
@@ -171,25 +167,20 @@ Foam::radiation::multiBandZoneAbsorptionEmission::eCont
         )
     );
 
-    volScalarField& e = te.ref();
+    scalarField& e = te.ref().primitiveFieldRef();
 
-    forAll(zoneCells_, zonei)
+    for (const label zonei : zoneIds_)
     {
-        const cellZone& cZone = mesh().cellZones()[zoneCells_[zonei]];
+        const cellZone& zn = mesh().cellZones()[zonei];
+        const auto iter = zoneEmissivity_.cfind(zn.name());
 
-        tmp<volScalarField> tzoneEm(e*0.0);
+        if (iter.good())  // Check is redundant (cannot fail)
-        volScalarField& zoneEm = tzoneEm.ref();
-
-        const scalarList& emi = zoneEmisivity_.find(cZone.name())();
-
-        forAll(cZone, i)
         {
-            label cellId = cZone[i];
+            const scalarList& emit = iter.val();
-            zoneEm[cellId] = emi[bandI] - emiCoeffs_[bandI];
-        }
-        e += zoneEm;
-    }
 
+            UIndirectList<scalar>(e, zn) = emit[bandI];
+        }
+    }
 
     return te;
 }

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/multiBandZoneAbsorptionEmission/multiBandZoneAbsorptionEmission.H

@@ -86,11 +86,11 @@ private:
         //- Cell zones absorptivity
         HashTable<scalarList> zoneAbsorptivity_;
 
-        //- Cell zones emisivity
+        //- Cell zones emissivity
-        HashTable<scalarList> zoneEmisivity_;
+        HashTable<scalarList> zoneEmissivity_;
 
-        //- Cells for each zone
+        //- The cellZones selected
-        labelList zoneCells_;
+        labelList zoneIds_;
 
 
 

src/waveModels/waveModel/waveModel.C

@@ -90,8 +90,7 @@ void Foam::waveModel::initialiseGeometry()
     }
 
     // Local face centres
-    const vectorField& Cf = patch_.faceCentres();
+    const vectorField CfLocal(Rgl_ & patch_.faceCentres());
-    const vectorField CfLocal(Rgl_ & Cf);
     z_ = CfLocal.component(2);
 
     // Local face extents in z-direction