STYLE: wall functions: use more modern coding practices

The most frequent changes have been as follows.

from:

    tmp<scalarField> tuTau(new scalarField(patch().size(), Zero));
    scalarField& uTau = tuTau.ref();

to:

    auto tuTau = tmp<scalarField>::New(patch().size(), Zero);
    auto& uTau = tuTau.ref();

- Other changes involved the addition of - wherever approapriate -:

    const
    noexcept
    auto

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2020 OpenCFD Ltd
+    Copyright (C) 2017-2022 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -66,12 +66,13 @@ tmp<scalarField> alphatJayatillekeWallFunctionFvPatchScalarField::yPlus
 ) const
 {
     const label patchi = patch().index();
+
     const tmp<volScalarField> tnut = turbModel.nut();
     const volScalarField& nut = tnut();
 
     if (isA<nutWallFunctionFvPatchScalarField>(nut.boundaryField()[patchi]))
     {
-        const nutWallFunctionFvPatchScalarField& nutPf =
+        const auto& nutPf =
             dynamic_cast<const nutWallFunctionFvPatchScalarField&>
             (
                 nut.boundaryField()[patchi]
@@ -106,13 +107,13 @@ scalar alphatJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
     const scalar Prat
 ) const
 {
-    scalar ypt = 11.0;
+    scalar ypt = 11;
 
-    for (int i=0; i<maxIters_; i++)
+    for (int iter = 0; iter < maxIters_; ++iter)
     {
-        scalar f = ypt - (log(E_*ypt)/kappa_ + P)/Prat;
+        const scalar f = ypt - (log(E_*ypt)/kappa_ + P)/Prat;
-        scalar df = 1.0 - 1.0/(ypt*kappa_*Prat);
+        const scalar df = 1.0 - 1.0/(ypt*kappa_*Prat);
-        scalar yptNew = ypt - f/df;
+        const scalar yptNew = ypt - f/df;
 
         if (yptNew < VSMALL)
         {
@@ -226,15 +227,14 @@ void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
     const label patchi = patch().index();
 
     // Retrieve turbulence properties from model
-    const compressible::turbulenceModel& turbModel =
+    const auto& turbModel = db().lookupObject<compressible::turbulenceModel>
-        db().lookupObject<compressible::turbulenceModel>
+    (
+        IOobject::groupName
         (
-            IOobject::groupName
+            compressible::turbulenceModel::propertiesName,
-            (
+            internalField().group()
-                compressible::turbulenceModel::propertiesName,
+        )
-                internalField().group()
+    );
-            )
-        );
 
     const scalarField yPlusp(yPlus(turbModel));
 
@@ -265,42 +265,45 @@ void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
     // Populate boundary values
     forAll(alphatw, facei)
     {
-        scalar yPlus = yPlusp[facei];
+        const scalar yPlus = yPlusp[facei];
 
-        scalar uTau = yPlus/y[facei]*(muw[facei]/rhow[facei]);
+        const scalar uTau = yPlus/y[facei]*(muw[facei]/rhow[facei]);
 
         // Molecular Prandtl number
-        scalar Pr = muw[facei]/alphaw[facei];
+        const scalar Pr = muw[facei]/alphaw[facei];
 
         // Molecular-to-turbulent Prandtl number ratio
-        scalar Prat = Pr/Prt_;
+        const scalar Prat = Pr/Prt_;
 
         // Thermal sublayer thickness
-        scalar P = Psmooth(Prat);
+        const scalar P = Psmooth(Prat);
-        scalar yPlusTherm = this->yPlusTherm(P, Prat);
+        const scalar yPlusTherm = this->yPlusTherm(P, Prat);
 
         // Evaluate new effective thermal diffusivity
-        scalar alphaEff = 0.0;
+        scalar alphaEff = 0;
         if (yPlus < yPlusTherm)
         {
-            scalar A = qDot[facei]*rhow[facei]*uTau*y[facei];
+            const scalar A = qDot[facei]*rhow[facei]*uTau*y[facei];
-            scalar B = qDot[facei]*Pr*yPlus;
+            const scalar B = qDot[facei]*Pr*yPlus;
-            scalar C = Pr*0.5*rhow[facei]*uTau*sqr(magUp[facei]);
+            const scalar C = Pr*0.5*rhow[facei]*uTau*sqr(magUp[facei]);
+
             alphaEff = A/(B + C + VSMALL);
         }
         else
         {
-            scalar A = qDot[facei]*rhow[facei]*uTau*y[facei];
+            const scalar A = qDot[facei]*rhow[facei]*uTau*y[facei];
-            scalar B = qDot[facei]*Prt_*(1.0/kappa_*log(E_*yPlus) + P);
+            const scalar B = qDot[facei]*Prt_*(1.0/kappa_*log(E_*yPlus) + P);
-            scalar magUc = uTau/kappa_*log(E_*yPlusTherm) - mag(Uw[facei]);
+            const scalar magUc =
-            scalar C =
+                uTau/kappa_*log(E_*yPlusTherm) - mag(Uw[facei]);
+            const scalar C =
                 0.5*rhow[facei]*uTau
                *(Prt_*sqr(magUp[facei]) + (Pr - Prt_)*sqr(magUc));
+
             alphaEff = A/(B + C + VSMALL);
         }
 
         // Update turbulent thermal diffusivity
-        alphatw[facei] = max(0.0, alphaEff - alphaw[facei]);
+        alphatw[facei] = max(scalar(0), alphaEff - alphaw[facei]);
 
         if (debug)
         {

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -110,15 +110,14 @@ void alphatWallFunctionFvPatchScalarField::updateCoeffs()
     const label patchi = patch().index();
 
     // Retrieve turbulence properties from model
-    const compressibleTurbulenceModel& turbModel =
+    const auto& turbModel = db().lookupObject<compressibleTurbulenceModel>
-        db().lookupObject<compressibleTurbulenceModel>
+    (
+        IOobject::groupName
         (
-            IOobject::groupName
+            compressibleTurbulenceModel::propertiesName,
-            (
+            internalField().group()
-                compressibleTurbulenceModel::propertiesName,
+        )
-                internalField().group()
+    );
-            )
-        );
 
     const scalarField& rhow = turbModel.rho().boundaryField()[patchi];
     const tmp<scalarField> tnutw = turbModel.nut(patchi);

src/TurbulenceModels/incompressible/turbulentTransportModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2020 OpenCFD Ltd
+    Copyright (C) 2017-2022 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -72,7 +72,7 @@ tmp<scalarField> alphatJayatillekeWallFunctionFvPatchScalarField::yPlus
 
     if (isA<nutWallFunctionFvPatchScalarField>(nut.boundaryField()[patchi]))
     {
-        const nutWallFunctionFvPatchScalarField& nutPf =
+        const auto& nutPf =
             dynamic_cast<const nutWallFunctionFvPatchScalarField&>
             (
                 nut.boundaryField()[patchi]
@@ -107,13 +107,13 @@ scalar alphatJayatillekeWallFunctionFvPatchScalarField::yPlusTherm
     const scalar Prat
 ) const
 {
-    scalar ypt = 11.0;
+    scalar ypt = 11;
 
-    for (int i=0; i<maxIters_; i++)
+    for (int iter = 0; iter < maxIters_; ++iter)
     {
-        scalar f = ypt - (log(E_*ypt)/kappa_ + P)/Prat;
+        const scalar f = ypt - (log(E_*ypt)/kappa_ + P)/Prat;
-        scalar df = 1.0 - 1.0/(ypt*kappa_*Prat);
+        const scalar df = 1.0 - 1.0/(ypt*kappa_*Prat);
-        scalar yptNew = ypt - f/df;
+        const scalar yptNew = ypt - f/df;
 
         if (yptNew < VSMALL)
         {
@@ -230,7 +230,7 @@ void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
 
     // Retrieve turbulence properties from model
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -245,7 +245,7 @@ void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
     const volScalarField& nu = tnu();
     const scalarField& nuw = nu.boundaryField()[patchi];
 
-    const IOdictionary& transportProperties =
+    const auto& transportProperties =
         db().lookupObject<IOdictionary>("transportProperties");
 
     // Molecular Prandtl number
@@ -258,21 +258,23 @@ void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
     scalarField& alphatw = *this;
     forAll(alphatw, facei)
     {
-        scalar yPlus = yPlusp[facei];
+        const scalar yPlus = yPlusp[facei];
 
         // Molecular-to-turbulent Prandtl number ratio
-        scalar Prat = Pr/Prt_;
+        const scalar Prat = Pr/Prt_;
 
         // Thermal sublayer thickness
-        scalar P = Psmooth(Prat);
+        const scalar P = Psmooth(Prat);
-        scalar yPlusTherm = this->yPlusTherm(P, Prat);
+        const scalar yPlusTherm = this->yPlusTherm(P, Prat);
 
         // Update turbulent thermal conductivity
         if (yPlus > yPlusTherm)
         {
-            scalar nu = nuw[facei];
+            const scalar nu = nuw[facei];
-            scalar kt = nu*(yPlus/(Prt_*(log(E_*yPlus)/kappa_ + P)) - 1/Pr);
+            const scalar kt =
-            alphatw[facei] = max(0.0, kt);
+                nu*(yPlus/(Prt_*(log(E_*yPlus)/kappa_ + P)) - 1.0/Pr);
+
+            alphatw[facei] = max(scalar(0), kt);
         }
         else
         {

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

@@ -45,7 +45,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::setMaster()
         return;
     }
 
-    const volScalarField& epsilon =
+    const auto& epsilon =
         static_cast<const volScalarField&>(this->internalField());
 
     const volScalarField::Boundary& bf = epsilon.boundaryField();
@@ -70,7 +70,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::setMaster()
 
 void Foam::epsilonWallFunctionFvPatchScalarField::createAveragingWeights()
 {
-    const volScalarField& epsilon =
+    const auto& epsilon =
         static_cast<const volScalarField&>(this->internalField());
 
     const volScalarField::Boundary& bf = epsilon.boundaryField();
@@ -133,12 +133,12 @@ Foam::epsilonWallFunctionFvPatchScalarField::epsilonPatch
     const label patchi
 )
 {
-    const volScalarField& epsilon =
+    const auto& epsilon =
         static_cast<const volScalarField&>(this->internalField());
 
     const volScalarField::Boundary& bf = epsilon.boundaryField();
 
-    const epsilonWallFunctionFvPatchScalarField& epf =
+    const auto& epf =
         refCast<const epsilonWallFunctionFvPatchScalarField>(bf[patchi]);
 
     return const_cast<epsilonWallFunctionFvPatchScalarField&>(epf);
@@ -222,8 +222,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate
         const scalar epsilonVis = w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
 
         // Contribution from the inertial sublayer
-        const scalar epsilonLog =
+        const scalar epsilonLog = w*Cmu75*pow(k[celli], 1.5)/(kappa*y[facei]);
-            w*Cmu75*pow(k[celli], 1.5)/(kappa*y[facei]);
 
         switch (blender_)
         {
@@ -455,7 +454,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -503,7 +502,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::updateWeightedCoeffs
         return;
     }
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kLowReWallFunction/kLowReWallFunctionFvPatchScalarField.C

@@ -144,7 +144,7 @@ void Foam::kLowReWallFunctionFvPatchScalarField::updateCoeffs()
 
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutLowReWallFunction/nutLowReWallFunctionFvPatchScalarField.C

@@ -107,7 +107,8 @@ Foam::tmp<Foam::scalarField> Foam::nutLowReWallFunctionFvPatchScalarField::
 yPlus() const
 {
     const label patchi = patch().index();
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -115,6 +116,7 @@ yPlus() const
             internalField().group()
         )
     );
+
     const scalarField& y = turbModel.y()[patchi];
 
     const tmp<scalarField> tnuw = turbModel.nu(patchi);

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUBlendedWallFunction/nutUBlendedWallFunctionFvPatchScalarField.C

@@ -38,7 +38,7 @@ Foam::nutUBlendedWallFunctionFvPatchScalarField::calcNut() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -46,8 +46,10 @@ Foam::nutUBlendedWallFunctionFvPatchScalarField::calcNut() const
             internalField().group()
         )
     );
+
     const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
     const scalarField magGradU(mag(Uw.snGrad()));
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
@@ -67,7 +69,7 @@ Foam::nutUBlendedWallFunctionFvPatchScalarField::calcUTau
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -87,8 +89,8 @@ Foam::nutUBlendedWallFunctionFvPatchScalarField::calcUTau
     Up -= n*(n & Up);
     const scalarField magUp(mag(Up));
 
-    tmp<scalarField> tuTaup(new scalarField(patch().size(), Zero));
+    auto tuTaup = tmp<scalarField>::New(patch().size(), Zero);
-    scalarField& uTaup = tuTaup.ref();
+    auto& uTaup = tuTaup.ref();
 
     const scalarField& nutw = *this;
 
@@ -202,7 +204,7 @@ Foam::nutUBlendedWallFunctionFvPatchScalarField::yPlus() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -210,9 +212,12 @@ Foam::nutUBlendedWallFunctionFvPatchScalarField::yPlus() const
             internalField().group()
         )
     );
+
     const scalarField& y = turbModel.y()[patchi];
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
+
     const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
     const scalarField magGradU(mag(Uw.snGrad()));
 

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutURoughWallFunction/nutURoughWallFunctionFvPatchScalarField.C

@@ -40,7 +40,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcNut() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -48,8 +48,11 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcNut() const
             internalField().group()
         )
     );
+
     const scalarField& y = turbModel.y()[patchi];
+
     const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
@@ -61,8 +64,8 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcNut() const
     tmp<scalarField> tyPlus = calcYPlus(magUp);
     scalarField& yPlus = tyPlus.ref();
 
-    tmp<scalarField> tnutw(new scalarField(patch().size(), Zero));
+    auto tnutw = tmp<scalarField>::New(patch().size(), Zero);
-    scalarField& nutw = tnutw.ref();
+    auto& nutw = tnutw.ref();
 
     forAll(yPlus, facei)
     {
@@ -85,7 +88,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -93,7 +96,9 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
             internalField().group()
         )
     );
+
     const scalarField& y = turbModel.y()[patchi];
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
@@ -101,18 +106,17 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
     const scalar E = wallCoeffs_.E();
     const scalar yPlusLam = wallCoeffs_.yPlusLam();
 
-    tmp<scalarField> tyPlus(new scalarField(patch().size(), Zero));
+    auto tyPlus = tmp<scalarField>::New(patch().size(), Zero);
-    scalarField& yPlus = tyPlus.ref();
+    auto& yPlus = tyPlus.ref();
 
-    if (0.0 < roughnessHeight_)
+    if (roughnessHeight_ > 0.0)
     {
         // Rough Walls
-        const scalar c_1 = 1/(90 - 2.25) + roughnessConstant_;
+        const scalar c_1 = 1.0/(90.0 - 2.25) + roughnessConstant_;
-        static const scalar c_2 = 2.25/(90 - 2.25);
+        static const scalar c_2 = 2.25/(90.0 - 2.25);
-        static const scalar c_3 = 2.0*atan(1.0)/log(90/2.25);
+        static const scalar c_3 = 2.0*atan(1.0)/log(90.0/2.25);
         static const scalar c_4 = c_3*log(2.25);
 
-        //if (KsPlusBasedOnYPlus_)
         {
             // If KsPlus is based on YPlus the extra term added to the law
             // of the wall will depend on yPlus
@@ -137,14 +141,14 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
                     yPlusLast = yp;
 
                     // The non-dimensional roughness height
-                    scalar KsPlus = yp*dKsPlusdYPlus;
+                    const scalar KsPlus = yp*dKsPlusdYPlus;
 
                     // The extra term in the law-of-the-wall
-                    scalar G = 0.0;
+                    scalar G = 0;
 
-                    scalar yPlusGPrime = 0.0;
+                    scalar yPlusGPrime = 0;
 
-                    if (KsPlus >= 90)
+                    if (KsPlus >= 90.0)
                     {
                         const scalar t_1 = 1 + roughnessConstant_*KsPlus;
                         G = log(t_1);
@@ -161,7 +165,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
                             (c_1*sint_2*KsPlus/t_1) + (c_3*logt_1*cos(t_2));
                     }
 
-                    scalar denom = 1.0 + log(E*yp) - G - yPlusGPrime;
+                    const scalar denom = 1.0 + log(E*yp) - G - yPlusGPrime;
                     if (mag(denom) > VSMALL)
                     {
                         yp = (kappaRe + yp*(1 - yPlusGPrime))/denom;
@@ -173,7 +177,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
                  && yp > VSMALL
                 );
 
-                yPlus[facei] = max(0.0, yp);
+                yPlus[facei] = max(scalar(0), yp);
             }
         }
     }
@@ -190,7 +194,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
             const scalar ryPlusLam = 1.0/yp;
 
             int iter = 0;
-            scalar yPlusLast = 0.0;
+            scalar yPlusLast = 0;
 
             do
             {
@@ -204,7 +208,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::calcYPlus
              && ++iter < maxIter_
             );
 
-            yPlus[facei] = max(0.0, yp);
+            yPlus[facei] = max(scalar(0), yp);
         }
     }
 
@@ -316,7 +320,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::yPlus() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -324,6 +328,7 @@ Foam::nutURoughWallFunctionFvPatchScalarField::yPlus() const
             internalField().group()
         )
     );
+
     const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
     tmp<scalarField> magUp = mag(Uw.patchInternalField() - Uw);
 

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutURoughWallFunction/nutURoughWallFunctionFvPatchScalarField.H

@@ -205,38 +205,38 @@ public:
         // Access
 
             //- Return the roughness height
-            scalar roughnessHeight() const
+            scalar roughnessHeight() const noexcept
             {
                 return roughnessHeight_;
             }
 
             //- Return reference to the roughness height to allow adjustment
-            scalar& roughnessHeight()
+            scalar& roughnessHeight() noexcept
             {
                 return roughnessHeight_;
             }
 
             //- Return the roughness constant scale
-            scalar roughnessConstant() const
+            scalar roughnessConstant() const noexcept
             {
                 return roughnessConstant_;
             }
 
             //- Return reference to the roughness constant to allow adjustment
-            scalar& roughnessConstant()
+            scalar& roughnessConstant() noexcept
             {
                 return roughnessConstant_;
             }
 
             //- Return the roughness scale factor
-            scalar roughnessFactor() const
+            scalar roughnessFactor() const noexcept
             {
                 return roughnessFactor_;
             }
 
             //- Return reference to the roughness scale factor to allow
             //- adjustment
-            scalar& roughnessFactor()
+            scalar& roughnessFactor() noexcept
             {
                 return roughnessFactor_;
             }

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUSpaldingWallFunction/nutUSpaldingWallFunctionFvPatchScalarField.C

@@ -40,7 +40,7 @@ Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcNut() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -48,8 +48,10 @@ Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcNut() const
             internalField().group()
         )
     );
+
     const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
     const scalarField magGradU(mag(Uw.snGrad()));
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
@@ -110,7 +112,7 @@ Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -118,6 +120,7 @@ Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
             internalField().group()
         )
     );
+
     const scalarField& y = turbModel.y()[patchi];
 
     const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
@@ -131,8 +134,8 @@ Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
 
     const scalarField& nutw = *this;
 
-    tmp<scalarField> tuTau(new scalarField(patch().size(), Zero));
+    auto tuTau = tmp<scalarField>::New(patch().size(), Zero);
-    scalarField& uTau = tuTau.ref();
+    auto& uTau = tuTau.ref();
 
     err.setSize(uTau.size());
     err = 0.0;
@@ -149,20 +152,20 @@ Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
 
             do
             {
-                scalar kUu = min(kappa*magUp[facei]/ut, 50);
+                const scalar kUu = min(kappa*magUp[facei]/ut, scalar(50));
-                scalar fkUu = exp(kUu) - 1 - kUu*(1 + 0.5*kUu);
+                const scalar fkUu = exp(kUu) - 1 - kUu*(1 + 0.5*kUu);
 
-                scalar f =
+                const scalar f =
                     - ut*y[facei]/nuw[facei]
                     + magUp[facei]/ut
-                    + 1/E*(fkUu - 1.0/6.0*kUu*sqr(kUu));
+                    + 1.0/E*(fkUu - 1.0/6.0*kUu*sqr(kUu));
 
-                scalar df =
+                const scalar df =
                     y[facei]/nuw[facei]
                   + magUp[facei]/sqr(ut)
-                  + 1/E*kUu*fkUu/ut;
+                  + 1.0/E*kUu*fkUu/ut;
 
-                scalar uTauNew = ut + f/df;
+                const scalar uTauNew = ut + f/df;
                 err[facei] = mag((ut - uTauNew)/ut);
                 ut = uTauNew;
 
@@ -175,7 +178,7 @@ Foam::nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
              && ++iter < maxIter
             );
 
-            uTau[facei] = max(0.0, ut);
+            uTau[facei] = max(scalar(0), ut);
 
             //invocations_++;
             //if (iter > 1)
@@ -320,7 +323,7 @@ Foam::nutUSpaldingWallFunctionFvPatchScalarField::yPlus() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -328,8 +331,11 @@ Foam::nutUSpaldingWallFunctionFvPatchScalarField::yPlus() const
             internalField().group()
         )
     );
+
     const scalarField& y = turbModel.y()[patchi];
+
     const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUTabulatedWallFunction/nutUTabulatedWallFunctionFvPatchScalarField.C

@@ -39,7 +39,7 @@ Foam::nutUTabulatedWallFunctionFvPatchScalarField::calcNut() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -47,10 +47,13 @@ Foam::nutUTabulatedWallFunctionFvPatchScalarField::calcNut() const
             internalField().group()
         )
     );
+
     const scalarField& y = turbModel.y()[patchi];
+
     const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
     const scalarField magUp(mag(Uw.patchInternalField() - Uw));
     const scalarField magGradU(mag(Uw.snGrad()));
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
@@ -71,8 +74,8 @@ Foam::nutUTabulatedWallFunctionFvPatchScalarField::calcUPlus
     const scalarField& Rey
 ) const
 {
-    tmp<scalarField> tuPlus(new scalarField(patch().size(), Zero));
+    auto tuPlus = tmp<scalarField>::New(patch().size(), Zero);
-    scalarField& uPlus = tuPlus.ref();
+    auto& uPlus = tuPlus.ref();
 
     forAll(uPlus, facei)
     {
@@ -192,7 +195,7 @@ Foam::nutUTabulatedWallFunctionFvPatchScalarField::yPlus() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -200,11 +203,15 @@ Foam::nutUTabulatedWallFunctionFvPatchScalarField::yPlus() const
             internalField().group()
         )
     );
+
     const scalarField& y = turbModel.y()[patchi];
+
     const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
     const scalarField magUp(mag(Uw.patchInternalField() - Uw));
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
+
     const scalarField Rey(magUp*y/nuw);
 
     return Rey/(calcUPlus(Rey) + ROOTVSMALL);

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

@@ -40,7 +40,7 @@ Foam::nutUWallFunctionFvPatchScalarField::calcNut() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -48,8 +48,10 @@ Foam::nutUWallFunctionFvPatchScalarField::calcNut() const
             internalField().group()
         )
     );
+
     const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
     const scalarField magUp(mag(Uw.patchInternalField() - Uw));
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
@@ -60,8 +62,8 @@ Foam::nutUWallFunctionFvPatchScalarField::calcNut() const
     tmp<scalarField> tyPlus = calcYPlus(magUp);
     const scalarField& yPlus = tyPlus();
 
-    tmp<scalarField> tnutw(new scalarField(patch().size(), Zero));
+    auto tnutw = tmp<scalarField>::New(patch().size(), Zero);
-    scalarField& nutw = tnutw.ref();
+    auto& nutw = tnutw.ref();
 
     forAll(yPlus, facei)
     {
@@ -144,7 +146,7 @@ Foam::nutUWallFunctionFvPatchScalarField::calcYPlus
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -152,7 +154,9 @@ Foam::nutUWallFunctionFvPatchScalarField::calcYPlus
             internalField().group()
         )
     );
+
     const scalarField& y = turbModel.y()[patchi];
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
@@ -160,8 +164,8 @@ Foam::nutUWallFunctionFvPatchScalarField::calcYPlus
     const scalar E = wallCoeffs_.E();
     const scalar yPlusLam = wallCoeffs_.yPlusLam();
 
-    tmp<scalarField> tyPlus(new scalarField(patch().size(), Zero));
+    auto tyPlus = tmp<scalarField>::New(patch().size(), Zero);
-    scalarField& yPlus = tyPlus.ref();
+    auto& yPlus = tyPlus.ref();
 
     forAll(yPlus, facei)
     {
@@ -180,7 +184,7 @@ Foam::nutUWallFunctionFvPatchScalarField::calcYPlus
 
         } while (mag(ryPlusLam*(yp - yPlusLast)) > 0.01 && ++iter < 10 );
 
-        yPlus[facei] = max(0.0, yp);
+        yPlus[facei] = max(scalar(0), yp);
     }
 
     return tyPlus;

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutWallFunction/nutWallFunctionFvPatchScalarField.C

@@ -65,10 +65,8 @@ const Foam::volVectorField& Foam::nutWallFunctionFvPatchScalarField::U
     {
         return turb.U();
     }
-    else
+
-    {
+    return db().lookupObject<volVectorField>(UName_);
-        return db().lookupObject<volVectorField>(UName_);
-    }
 }
 
 

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C

@@ -50,10 +50,8 @@ Foam::scalar Foam::nutkRoughWallFunctionFvPatchScalarField::fnRough
             sin(0.4258*(log(KsPlus) - 0.811))
         );
     }
-    else
+
-    {
+    return (1.0 + Cs*KsPlus);
-        return (1.0 + Cs*KsPlus);
-    }
 }
 
 
@@ -62,7 +60,7 @@ calcNut() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -70,9 +68,12 @@ calcNut() const
             internalField().group()
         )
     );
+
     const scalarField& y = turbModel.y()[patchi];
+
     const tmp<volScalarField> tk = turbModel.k();
     const volScalarField& k = tk();
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
@@ -80,8 +81,8 @@ calcNut() const
     const scalar kappa = wallCoeffs_.kappa();
     const scalar E = wallCoeffs_.E();
 
-    tmp<scalarField> tnutw(new scalarField(*this));
+    auto tnutw = tmp<scalarField>::New(*this);
-    scalarField& nutw = tnutw.ref();
+    auto& nutw = tnutw.ref();
 
     forAll(nutw, facei)
     {
@@ -226,7 +227,7 @@ void Foam::nutkRoughWallFunctionFvPatchScalarField::rmap
 {
     nutkWallFunctionFvPatchScalarField::rmap(ptf, addr);
 
-    const nutkRoughWallFunctionFvPatchScalarField& nrwfpsf =
+    const auto& nrwfpsf =
         refCast<const nutkRoughWallFunctionFvPatchScalarField>(ptf);
 
     Ks_.rmap(nrwfpsf.Ks_, addr);

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

@@ -40,7 +40,7 @@ calcNut() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -50,8 +50,10 @@ calcNut() const
     );
 
     const scalarField& y = turbModel.y()[patchi];
+
     const tmp<volScalarField> tk = turbModel.k();
     const volScalarField& k = tk();
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
@@ -60,8 +62,8 @@ calcNut() const
     const scalar E = wallCoeffs_.E();
     const scalar yPlusLam = wallCoeffs_.yPlusLam();
 
-    tmp<scalarField> tnutw(new scalarField(patch().size(), Zero));
+    auto tnutw = tmp<scalarField>::New(patch().size(), Zero);
-    scalarField& nutw = tnutw.ref();
+    auto& nutw = tnutw.ref();
 
     forAll(nutw, facei)
     {

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C

@@ -45,7 +45,7 @@ void Foam::omegaWallFunctionFvPatchScalarField::setMaster()
         return;
     }
 
-    const volScalarField& omega =
+    const auto& omega =
         static_cast<const volScalarField&>(this->internalField());
 
     const volScalarField::Boundary& bf = omega.boundaryField();
@@ -70,7 +70,7 @@ void Foam::omegaWallFunctionFvPatchScalarField::setMaster()
 
 void Foam::omegaWallFunctionFvPatchScalarField::createAveragingWeights()
 {
-    const volScalarField& omega =
+    const auto& omega =
         static_cast<const volScalarField&>(this->internalField());
 
     const volScalarField::Boundary& bf = omega.boundaryField();
@@ -132,12 +132,12 @@ Foam::omegaWallFunctionFvPatchScalarField::omegaPatch
     const label patchi
 )
 {
-    const volScalarField& omega =
+    const auto& omega =
         static_cast<const volScalarField&>(this->internalField());
 
     const volScalarField::Boundary& bf = omega.boundaryField();
 
-    const omegaWallFunctionFvPatchScalarField& opf =
+    const auto& opf =
         refCast<const omegaWallFunctionFvPatchScalarField>(bf[patchi]);
 
     return const_cast<omegaWallFunctionFvPatchScalarField&>(opf);
@@ -445,7 +445,7 @@ void Foam::omegaWallFunctionFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -493,7 +493,7 @@ void Foam::omegaWallFunctionFvPatchScalarField::updateWeightedCoeffs
         return;
     }
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (

src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmAlphatkWallFunction/atmAlphatkWallFunctionFvPatchScalarField.C

@@ -41,10 +41,8 @@ namespace Foam
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 scalar atmAlphatkWallFunctionFvPatchScalarField::tolerance_ = 0.01;
-
 label atmAlphatkWallFunctionFvPatchScalarField::maxIters_ = 10;
 
-
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 void atmAlphatkWallFunctionFvPatchScalarField::checkType()
@@ -206,15 +204,14 @@ void atmAlphatkWallFunctionFvPatchScalarField::updateCoeffs()
     const label patchi = patch().index();
 
     // Retrieve turbulence properties from model
-    const auto& turbModel =
+    const auto& turbModel = db().lookupObject<turbulenceModel>
-        db().lookupObject<turbulenceModel>
+    (
+        IOobject::groupName
         (
-            IOobject::groupName
+            turbulenceModel::propertiesName,
-            (
+            internalField().group()
-                turbulenceModel::propertiesName,
+        )
-                internalField().group()
+    );
-            )
-        );
 
     const scalarField& y = turbModel.y()[patchi];
 
@@ -286,8 +283,15 @@ void atmAlphatkWallFunctionFvPatchScalarField::autoMap
 )
 {
     fixedValueFvPatchScalarField::autoMap(m);
-    Prt_->autoMap(m);
+
-    z0_->autoMap(m);
+    if (Prt_)
+    {
+        Prt_->autoMap(m);
+    }
+    if (z0_)
+    {
+        z0_->autoMap(m);
+    }
 }
 
 
@@ -299,11 +303,17 @@ void atmAlphatkWallFunctionFvPatchScalarField::rmap
 {
     fixedValueFvPatchScalarField::rmap(ptf, addr);
 
-    const atmAlphatkWallFunctionFvPatchScalarField& nrwfpsf =
+    const auto& nrwfpsf =
         refCast<const atmAlphatkWallFunctionFvPatchScalarField>(ptf);
 
-    z0_->rmap(nrwfpsf.z0_(), addr);
+    if (Prt_)
-    Prt_->rmap(nrwfpsf.Prt_(), addr);
+    {
+        Prt_->rmap(nrwfpsf.Prt_(), addr);
+    }
+    if (z0_)
+    {
+        z0_->rmap(nrwfpsf.z0_(), addr);
+    }
 }
 
 

src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmEpsilonWallFunction/atmEpsilonWallFunctionFvPatchScalarField.C

@@ -204,7 +204,11 @@ void Foam::atmEpsilonWallFunctionFvPatchScalarField::autoMap
 )
 {
     epsilonWallFunctionFvPatchScalarField::autoMap(m);
-    z0_->autoMap(m);
+
+    if (z0_)
+    {
+        z0_->autoMap(m);
+    }
 }
 
 
@@ -216,10 +220,12 @@ void Foam::atmEpsilonWallFunctionFvPatchScalarField::rmap
 {
     epsilonWallFunctionFvPatchScalarField::rmap(ptf, addr);
 
-    const atmEpsilonWallFunctionFvPatchScalarField& atmpsf =
+    const auto& atmpsf =
         refCast<const atmEpsilonWallFunctionFvPatchScalarField>(ptf);
-
+    if (z0_)
-    z0_->rmap(atmpsf.z0_(), addr);
+    {
+        z0_->rmap(atmpsf.z0_(), addr);
+    }
 }
 
 

src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutUWallFunction/atmNutUWallFunctionFvPatchScalarField.C

@@ -43,15 +43,14 @@ tmp<scalarField> atmNutUWallFunctionFvPatchScalarField::calcNut() const
 {
     const label patchi = patch().index();
 
-    const auto& turbModel =
+    const auto& turbModel = db().lookupObject<turbulenceModel>
-        db().lookupObject<turbulenceModel>
+    (
+        IOobject::groupName
         (
-            IOobject::groupName
+            turbulenceModel::propertiesName,
-            (
+            internalField().group()
-                turbulenceModel::propertiesName,
+        )
-                internalField().group()
+    );
-            )
-        );
 
     const scalarField& y = turbModel.y()[patchi];
 
@@ -186,7 +185,11 @@ void atmNutUWallFunctionFvPatchScalarField::autoMap
 )
 {
     nutUWallFunctionFvPatchScalarField::autoMap(m);
-    z0_->autoMap(m);
+
+    if (z0_)
+    {
+        z0_->autoMap(m);
+    }
 }
 
 
@@ -201,7 +204,10 @@ void atmNutUWallFunctionFvPatchScalarField::rmap
     const atmNutUWallFunctionFvPatchScalarField& nrwfpsf =
         refCast<const atmNutUWallFunctionFvPatchScalarField>(ptf);
 
-    z0_->rmap(nrwfpsf.z0_(), addr);
+    if (z0_)
+    {
+        z0_->rmap(nrwfpsf.z0_(), addr);
+    }
 }
 
 

src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutWallFunction/atmNutWallFunctionFvPatchScalarField.C

@@ -44,15 +44,15 @@ tmp<scalarField> atmNutWallFunctionFvPatchScalarField::calcNut() const
 {
     const label patchi = patch().index();
 
-    const auto& turbModel =
+    const auto& turbModel = db().lookupObject<turbulenceModel>
-        db().lookupObject<turbulenceModel>
+    (
+        IOobject::groupName
         (
-            IOobject::groupName
+            turbulenceModel::propertiesName,
-            (
+            internalField().group()
-                turbulenceModel::propertiesName,
+        )
-                internalField().group()
+    );
-            )
+
-        );
     const scalarField& y = turbModel.y()[patchi];
 
     const tmp<volScalarField> tk = turbModel.k();
@@ -206,7 +206,11 @@ void atmNutWallFunctionFvPatchScalarField::autoMap
 )
 {
     nutkWallFunctionFvPatchScalarField::autoMap(m);
-    z0_->autoMap(m);
+
+    if (z0_)
+    {
+        z0_->autoMap(m);
+    }
 }
 
 
@@ -218,10 +222,13 @@ void atmNutWallFunctionFvPatchScalarField::rmap
 {
     nutkWallFunctionFvPatchScalarField::rmap(ptf, addr);
 
-    const atmNutWallFunctionFvPatchScalarField& nrwfpsf =
+    const auto& nrwfpsf =
         refCast<const atmNutWallFunctionFvPatchScalarField>(ptf);
 
-    z0_->rmap(nrwfpsf.z0_(), addr);
+    if (z0_)
+    {
+        z0_->rmap(nrwfpsf.z0_(), addr);
+    }
 }
 
 

src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmNutkWallFunction/atmNutkWallFunctionFvPatchScalarField.C

@@ -45,15 +45,15 @@ tmp<scalarField> atmNutkWallFunctionFvPatchScalarField::calcNut() const
 {
     const label patchi = patch().index();
 
-    const auto& turbModel =
+    const auto& turbModel = db().lookupObject<turbulenceModel>
-        db().lookupObject<turbulenceModel>
+    (
+        IOobject::groupName
         (
-            IOobject::groupName
+            turbulenceModel::propertiesName,
-            (
+            internalField().group()
-                turbulenceModel::propertiesName,
+        )
-                internalField().group()
+    );
-            )
+
-        );
     const scalarField& y = turbModel.y()[patchi];
 
     const tmp<volScalarField> tk = turbModel.k();
@@ -193,7 +193,11 @@ void atmNutkWallFunctionFvPatchScalarField::autoMap
 )
 {
     nutkWallFunctionFvPatchScalarField::autoMap(m);
-    z0_->autoMap(m);
+
+    if (z0_)
+    {
+        z0_->autoMap(m);
+    }
 }
 
 
@@ -205,10 +209,13 @@ void atmNutkWallFunctionFvPatchScalarField::rmap
 {
     nutkWallFunctionFvPatchScalarField::rmap(ptf, addr);
 
-    const atmNutkWallFunctionFvPatchScalarField& nrwfpsf =
+    const auto& nrwfpsf =
         refCast<const atmNutkWallFunctionFvPatchScalarField>(ptf);
 
-    z0_->rmap(nrwfpsf.z0_(), addr);
+    if (z0_)
+    {
+        z0_->rmap(nrwfpsf.z0_(), addr);
+    }
 }
 
 

src/atmosphericModels/derivedFvPatchFields/wallFunctions/atmOmegaWallFunction/atmOmegaWallFunctionFvPatchScalarField.C

@@ -187,7 +187,11 @@ void Foam::atmOmegaWallFunctionFvPatchScalarField::autoMap
 )
 {
     omegaWallFunctionFvPatchScalarField::autoMap(m);
-    z0_->autoMap(m);
+
+    if (z0_)
+    {
+        z0_->autoMap(m);
+    }
 }
 
 
@@ -199,10 +203,13 @@ void Foam::atmOmegaWallFunctionFvPatchScalarField::rmap
 {
     omegaWallFunctionFvPatchScalarField::rmap(ptf, addr);
 
-    const atmOmegaWallFunctionFvPatchScalarField& atmpsf =
+    const auto& atmpsf =
         refCast<const atmOmegaWallFunctionFvPatchScalarField>(ptf);
 
-    z0_->rmap(atmpsf.z0_(), addr);
+    if (z0_)
+    {
+        z0_->rmap(atmpsf.z0_(), addr);
+    }
 }
 
 

src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/alphatFilmWallFunction/alphatFilmWallFunctionFvPatchScalarField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -168,12 +168,12 @@ void alphatFilmWallFunctionFvPatchScalarField::updateCoeffs()
     // Retrieve phase change mass from surface film model
     const label filmPatchi = filmModel.regionPatchID(patchi);
 
-    tmp<volScalarField> mDotFilm(filmModel.primaryMassTrans());
+    tmp<volScalarField> mDotFilm = filmModel.primaryMassTrans();
     scalarField mDotFilmp = mDotFilm().boundaryField()[filmPatchi];
     filmModel.toPrimary(filmPatchi, mDotFilmp);
 
     // Retrieve RAS turbulence model
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -183,49 +183,55 @@ void alphatFilmWallFunctionFvPatchScalarField::updateCoeffs()
     );
 
     const scalarField& y = turbModel.y()[patchi];
+
     const scalarField& rhow = turbModel.rho().boundaryField()[patchi];
+
     const tmp<volScalarField> tk = turbModel.k();
     const volScalarField& k = tk();
+
     const tmp<scalarField> tmuw = turbModel.mu(patchi);
     const scalarField& muw = tmuw();
+
     const tmp<scalarField> talpha = turbModel.alpha(patchi);
     const scalarField& alphaw = talpha();
 
-    const scalar Cmu25 = pow(Cmu_, 0.25);
+    const scalar Cmu25 = pow025(Cmu_);
 
     // Populate alphat field values
     scalarField& alphat = *this;
     forAll(alphat, facei)
     {
-        label faceCelli = patch().faceCells()[facei];
+        const label faceCelli = patch().faceCells()[facei];
 
-        scalar uTau = Cmu25*sqrt(k[faceCelli]);
+        const scalar uTau = Cmu25*sqrt(k[faceCelli]);
 
-        scalar yPlus = y[facei]*uTau/(muw[facei]/rhow[facei]);
+        const scalar yPlus = y[facei]*uTau/(muw[facei]/rhow[facei]);
 
-        scalar Pr = muw[facei]/alphaw[facei];
+        const scalar Pr = muw[facei]/alphaw[facei];
 
-        scalar factor = 0.0;
+        scalar factor = 0;
-        scalar mStar = mDotFilmp[facei]/(y[facei]*uTau);
+        const scalar mStar = mDotFilmp[facei]/(y[facei]*uTau);
         if (yPlus > yPlusCrit_)
         {
-            scalar expTerm = exp(min(50.0, yPlusCrit_*mStar*Pr));
+            const scalar expTerm = exp(min(scalar(50), yPlusCrit_*mStar*Pr));
-            scalar yPlusRatio = yPlus/yPlusCrit_;
+            const scalar yPlusRatio = yPlus/yPlusCrit_;
-            scalar powTerm = mStar*Prt_/kappa_;
+            const scalar powTerm = mStar*Prt_/kappa_;
+
             factor =
                 mStar/(expTerm*(pow(yPlusRatio, powTerm)) - 1.0 + ROOTVSMALL);
         }
         else
         {
-            scalar expTerm = exp(min(50.0, yPlus*mStar*Pr));
+            const scalar expTerm = exp(min(scalar(50), yPlus*mStar*Pr));
+
             factor = mStar/(expTerm - 1.0 + ROOTVSMALL);
         }
 
-        scalar dx = patch().deltaCoeffs()[facei];
+        const scalar dx = patch().deltaCoeffs()[facei];
 
-        scalar alphaEff = dx*rhow[facei]*uTau*factor;
+        const scalar alphaEff = dx*rhow[facei]*uTau*factor;
 
-        alphat[facei] = max(alphaEff - alphaw[facei], 0.0);
+        alphat[facei] = max(alphaEff - alphaw[facei], scalar(0));
     }
 
     // Restore tag

src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/nutkFilmWallFunction/nutkFilmWallFunctionFvPatchScalarField.C

@@ -51,8 +51,8 @@ tmp<scalarField> nutkFilmWallFunctionFvPatchScalarField::calcUTau
     const scalarField& magGradU
 ) const
 {
-    tmp<scalarField> tuTau(new scalarField(patch().size(), Zero));
+    auto tuTau = tmp<scalarField>::New(patch().size(), Zero);
-    scalarField& uTau = tuTau.ref();
+    auto& uTau = tuTau.ref();
 
     const auto* filmModelPtr = db().time().findObject
         <regionModels::surfaceFilmModels::surfaceFilmRegionModel>
@@ -71,13 +71,13 @@ tmp<scalarField> nutkFilmWallFunctionFvPatchScalarField::calcUTau
 
     const label filmPatchi = filmModel.regionPatchID(patchi);
 
-    tmp<volScalarField> mDotFilm(filmModel.primaryMassTrans());
+    tmp<volScalarField> mDotFilm = filmModel.primaryMassTrans();
     scalarField mDotFilmp = mDotFilm().boundaryField()[filmPatchi];
     filmModel.toPrimary(filmPatchi, mDotFilmp);
 
 
     // Retrieve RAS turbulence model
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -87,38 +87,41 @@ tmp<scalarField> nutkFilmWallFunctionFvPatchScalarField::calcUTau
     );
 
     const scalarField& y = turbModel.y()[patchi];
+
     const tmp<volScalarField> tk = turbModel.k();
     const volScalarField& k = tk();
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
-    const scalar Cmu25 = pow(wallCoeffs_.Cmu(), 0.25);
+    const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
     const scalar kappa = wallCoeffs_.kappa();
 
     forAll(uTau, facei)
     {
-        label faceCelli = patch().faceCells()[facei];
+        const label faceCelli = patch().faceCells()[facei];
 
-        scalar ut = Cmu25*sqrt(k[faceCelli]);
+        const scalar ut = Cmu25*sqrt(k[faceCelli]);
 
-        scalar yPlus = y[facei]*ut/nuw[facei];
+        const scalar yPlus = y[facei]*ut/nuw[facei];
 
-        scalar mStar = mDotFilmp[facei]/(y[facei]*ut);
+        const scalar mStar = mDotFilmp[facei]/(y[facei]*ut);
 
-        scalar factor = 0.0;
+        scalar factor = 0;
         if (yPlus > yPlusCrit_)
         {
-            scalar expTerm = exp(min(50.0, B_*mStar));
+            const scalar expTerm = exp(min(scalar(50), B_*mStar));
-            scalar powTerm = pow(yPlus, mStar/kappa);
+            const scalar powTerm = pow(yPlus, mStar/kappa);
             factor = mStar/(expTerm*powTerm - 1.0 + ROOTVSMALL);
         }
         else
         {
-            scalar expTerm = exp(min(50.0, mStar));
+            const scalar expTerm = exp(min(scalar(50), mStar));
+
             factor = mStar/(expTerm*yPlus - 1.0 + ROOTVSMALL);
         }
 
-        uTau[facei] = sqrt(max(0, magGradU[facei]*ut*factor));
+        uTau[facei] = sqrt(max(scalar(0), magGradU[facei]*ut*factor));
     }
 
     return tuTau;
@@ -129,7 +132,7 @@ tmp<scalarField> nutkFilmWallFunctionFvPatchScalarField::calcNut() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -140,6 +143,7 @@ tmp<scalarField> nutkFilmWallFunctionFvPatchScalarField::calcNut() const
 
     const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
     const scalarField magGradU(mag(Uw.snGrad()));
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
@@ -245,7 +249,7 @@ tmp<scalarField> nutkFilmWallFunctionFvPatchScalarField::yPlus() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -255,7 +259,9 @@ tmp<scalarField> nutkFilmWallFunctionFvPatchScalarField::yPlus() const
     );
 
     const scalarField& y = turbModel.y()[patchi];
+
     const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();