Merge branch '2770-wall-functions' into 'develop'

ENH: wall functions: swap the order of switch statements and for loops

See merge request Development/openfoam!635

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

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

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

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

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

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

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

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016, 2019 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -188,102 +188,146 @@ void Foam::omegaWallFunctionFvPatchScalarField::calculate
 {
     const label patchi = patch.index();
 
-    const tmp<scalarField> tnutw = turbModel.nut(patchi);
+    const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
-    const scalarField& nutw = tnutw();
+    const scalar kappa = wallCoeffs_.kappa();
+    const scalar yPlusLam = wallCoeffs_.yPlusLam();
 
     const scalarField& y = turbModel.y()[patchi];
 
+    const labelUList& faceCells = patch.faceCells();
+
+    const tmp<scalarField> tnutw = turbModel.nut(patchi);
+    const scalarField& nutw = tnutw();
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
     const tmp<volScalarField> tk = turbModel.k();
     const volScalarField& k = tk();
 
-    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
+    // Calculate y-plus
-
+    const auto yPlus = [&](const label facei) -> scalar
-    const scalarField magGradUw(mag(Uw.snGrad()));
-
-    const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
-    const scalar kappa = wallCoeffs_.kappa();
-    const scalar yPlusLam = wallCoeffs_.yPlusLam();
-
-    // Set omega and G
-    forAll(nutw, facei)
     {
-        const label celli = patch.faceCells()[facei];
+        return
-        const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei];
+        (
-        const scalar w = cornerWeights[facei];
+            Cmu25*y[facei]*sqrt(k[faceCells[facei]])/nuw[facei]
+        );
+    };
 
-        // Contribution from the viscous sublayer
+    // Contribution from the viscous sublayer
-        const scalar omegaVis = 6.0*nuw[facei]/(beta1_*sqr(y[facei]));
+    const auto omegaVis = [&](const label facei) -> scalar
+    {
+        return
+        (
+            cornerWeights[facei]*6.0*nuw[facei]/(beta1_*sqr(y[facei]))
+        );
+    };
 
-        // Contribution from the inertial sublayer
+    // Contribution from the inertial sublayer
-        const scalar omegaLog = sqrt(k[celli])/(Cmu25*kappa*y[facei]);
+    const auto omegaLog = [&](const label facei) -> scalar
+    {
+        return
+        (
+            cornerWeights[facei]*sqrt(k[faceCells[facei]])
+          / (Cmu25*kappa*y[facei])
+        );
+    };
 
-        switch (blender_)
+    switch (blender_)
+    {
+        case blenderType::STEPWISE:
         {
-            case blenderType::STEPWISE:
+            forAll(faceCells, facei)
             {
-                if (yPlus > yPlusLam)
+                if (yPlus(facei) > yPlusLam)
                 {
-                    omega0[celli] += w*omegaLog;
+                    omega0[faceCells[facei]] += omegaLog(facei);
                 }
                 else
                 {
-                    omega0[celli] += w*omegaVis;
+                    omega0[faceCells[facei]] += omegaVis(facei);
                 }
-                break;
-            }
-
-            case blenderType::BINOMIAL:
-            {
-                omega0[celli] +=
-                    w*pow
-                    (
-                        pow(omegaVis, n_) + pow(omegaLog, n_),
-                        scalar(1)/n_
-                    );
-                break;
-            }
-
-            case blenderType::MAX:
-            {
-                // (PH:Eq. 27)
-                omega0[celli] += max(omegaVis, omegaLog);
-                break;
-            }
-
-            case blenderType::EXPONENTIAL:
-            {
-                // (PH:Eq. 31)
-                const scalar Gamma = 0.01*pow4(yPlus)/(1 + 5*yPlus);
-                const scalar invGamma = scalar(1)/(Gamma + ROOTVSMALL);
-
-                omega0[celli] +=
-                    w*(omegaVis*exp(-Gamma) + omegaLog*exp(-invGamma));
-                break;
-            }
-
-            case blenderType::TANH:
-            {
-                // (KAS:Eqs. 33-34)
-                const scalar phiTanh = tanh(pow4(0.1*yPlus));
-                const scalar b1 = omegaVis + omegaLog;
-                const scalar b2 =
-                    pow(pow(omegaVis, 1.2) + pow(omegaLog, 1.2), 1.0/1.2);
-
-                omega0[celli] += phiTanh*b1 + (1 - phiTanh)*b2;
-                break;
             }
+            break;
         }
 
-        if (!(blender_ == blenderType::STEPWISE) || yPlus > yPlusLam)
+        case blenderType::BINOMIAL:
         {
-            G0[celli] +=
+            forAll(faceCells, facei)
-                w
+            {
+                omega0[faceCells[facei]] +=
+                    pow
+                    (
+                        pow(omegaVis(facei), n_) + pow(omegaLog(facei), n_),
+                        scalar(1)/n_
+                    );
+            }
+            break;
+        }
+
+        case blenderType::MAX:
+        {
+            forAll(faceCells, facei)
+            {
+                // (PH:Eq. 27)
+                omega0[faceCells[facei]] +=
+                    max(omegaVis(facei), omegaLog(facei));
+            }
+            break;
+        }
+
+        case blenderType::EXPONENTIAL:
+        {
+            forAll(faceCells, facei)
+            {
+                // (PH:Eq. 31)
+                const scalar yPlusFace = yPlus(facei);
+                const scalar Gamma = 0.01*pow4(yPlusFace)/(1 + 5*yPlusFace);
+                const scalar invGamma = scalar(1)/(Gamma + ROOTVSMALL);
+
+                omega0[faceCells[facei]] +=
+                (
+                    omegaVis(facei)*exp(-Gamma)
+                  + omegaLog(facei)*exp(-invGamma)
+                );
+            }
+            break;
+        }
+
+        case blenderType::TANH:
+        {
+            forAll(faceCells, facei)
+            {
+                // (KAS:Eqs. 33-34)
+                const scalar omegaVisFace = omegaVis(facei);
+                const scalar omegaLogFace = omegaLog(facei);
+                const scalar b1 = omegaVisFace + omegaLogFace;
+                const scalar b2 =
+                    pow
+                    (
+                        pow(omegaVisFace, 1.2) + pow(omegaLogFace, 1.2),
+                        1.0/1.2
+                    );
+                const scalar phiTanh = tanh(pow4(0.1*yPlus(facei)));
+
+                omega0[faceCells[facei]] += phiTanh*b1 + (1 - phiTanh)*b2;
+            }
+            break;
+        }
+    }
+
+    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
+    const scalarField magGradUw(mag(Uw.snGrad()));
+
+    forAll(faceCells, facei)
+    {
+        if (!(blender_ == blenderType::STEPWISE) || yPlus(facei) > yPlusLam)
+        {
+            G0[faceCells[facei]] +=
+                cornerWeights[facei]
                *(nutw[facei] + nuw[facei])
                *magGradUw[facei]
-               *Cmu25*sqrt(k[celli])
+               *Cmu25*sqrt(k[faceCells[facei]])
                /(kappa*y[facei]);
         }
     }

src/thermoTools/derivedFvPatchFields/wallFunctions/sorptionWallFunction/sorptionWallFunctionFvPatchScalarField.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2022-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -64,9 +64,16 @@ static scalar calcYStarLam
 
 tmp<scalarField> sorptionWallFunctionFvPatchScalarField::yPlus() const
 {
-    // Calculate fields of interest
     const label patchi = patch().index();
 
+    // Calculate the empirical constant given by (Jayatilleke, 1966) (FDC:Eq. 6)
+    const scalar Pc =
+        9.24*(pow(Sc_/Sct_, 0.75) - 1)*(1 + 0.28*exp(-0.007*Sc_/Sct_));
+    const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
+    const scalar kappa = wallCoeffs_.kappa();
+    const scalar E = wallCoeffs_.E();
+
+    // Calculate fields of interest
     const auto& k = db().lookupObject<volScalarField>(kName_);
     tmp<scalarField> tkwc = k.boundaryField()[patchi].patchInternalField();
     const scalarField& kwc = tkwc.cref();
@@ -79,86 +86,116 @@ tmp<scalarField> sorptionWallFunctionFvPatchScalarField::yPlus() const
     tmp<scalarField> tywc = y.boundaryField()[patchi].patchInternalField();
     const scalarField& ywc = tywc.cref();
 
+    // (FDC:Eq. 3)
+    const auto yStar = [&](const label facei) -> scalar
+    {
+        return
+        (
+            Cmu25*sqrt(kwc[facei])*ywc[facei]/nuwc[facei]
+        );
+    };
 
-    // Calculate the empirical constant given by (Jayatilleke, 1966) (FDC:Eq. 6)
+    // (FDC:Eq. 4)
-    const scalar Pc =
+    const auto yPlusVis = [&](const label facei) -> scalar
-        9.24*(pow(Sc_/Sct_, 0.75) - 1)*(1 + 0.28*exp(-0.007*Sc_/Sct_));
+    {
-    const scalar Cmu25 = pow025(wallCoeffs_.Cmu());
+        return (Sc_*yStar(facei));
-    const scalar kappa = wallCoeffs_.kappa();
+    };
-    const scalar E = wallCoeffs_.E();
 
-    auto tyPlus = tmp<scalarField>::New(patch().size(), Zero);
+    // (FDC:Eq. 5)
+    const auto yPlusLog = [&](const label facei) -> scalar
+    {
+        return
+        (
+            Sct_*(log(max(E*yStar(facei), 1 + 1e-4))/kappa + Pc)
+        );
+    };
+
+    auto tyPlus = tmp<scalarField>::New(patch().size());
     auto& yPlus = tyPlus.ref();
 
-    forAll(yPlus, facei)
+    switch (blender_)
     {
-        // (FDC:Eq. 3)
+        case blenderType::EXPONENTIAL:
-        const scalar yStar = Cmu25*sqrt(kwc[facei])*ywc[facei]/nuwc[facei];
-
-        // (FDC:Eq. 4)
-        const scalar yPlusVis = Sc_*yStar;
-
-        // (FDC:Eq. 5)
-        const scalar yPlusLog = Sct_*(log(max(E*yStar, 1 + 1e-4))/kappa + Pc);
-
-        switch (blender_)
         {
-            case blenderType::EXPONENTIAL:
+            forAll(yPlus, facei)
             {
                 // (FDC:Eq. 2)
+                const scalar yStarFace = yStar(facei);
                 const scalar Gamma =
-                    0.01*pow4(Sc_*yStar)/(1 + 5*pow3(Sc_)*yStar);
+                    0.01*pow4(Sc_*yStarFace)/(1 + 5*pow3(Sc_)*yStarFace);
                 const scalar invGamma = scalar(1)/max(Gamma, ROOTVSMALL);
 
                 // (FDC:Eq. 1)
-                yPlus[facei] = yPlusVis*exp(-Gamma) + yPlusLog*exp(-invGamma);
+                yPlus[facei] =
-                break;
+                (
+                    yPlusVis(facei)*exp(-Gamma)
+                  + yPlusLog(facei)*exp(-invGamma)
+                );
             }
+            break;
+        }
 
-            case blenderType::STEPWISE:
+        case blenderType::STEPWISE:
+        {
+            static const scalar yStarLam =
+                calcYStarLam(kappa, E, Sc_, Sct_, Pc);
+
+            forAll(yPlus, facei)
             {
-                static const scalar yStarLam =
-                    calcYStarLam(kappa, E, Sc_, Sct_, Pc);
-
                 // (F:Eq. 5.3)
-                if (yStar < yStarLam)
+                if (yStar(facei) < yStarLam)
                 {
-                    yPlus[facei] = yPlusVis;
+                    yPlus[facei] = yPlusVis(facei);
                 }
                 else
                 {
-                    yPlus[facei] = yPlusLog;
+                    yPlus[facei] = yPlusLog(facei);
                 }
-                break;
             }
+            break;
+        }
 
-            case blenderType::BINOMIAL:
+        case blenderType::BINOMIAL:
+        {
+            forAll(yPlus, facei)
             {
                 yPlus[facei] =
                     pow
                     (
-                        pow(yPlusVis, n_) + pow(yPlusLog, n_),
+                        pow(yPlusVis(facei), n_) + pow(yPlusLog(facei), n_),
                         scalar(1)/n_
                     );
-                break;
             }
+            break;
+        }
 
-            case blenderType::MAX:
+        case blenderType::MAX:
+        {
+            forAll(yPlus, facei)
             {
-                yPlus[facei] = max(yPlusVis, yPlusLog);
+                yPlus[facei] = max(yPlusVis(facei), yPlusLog(facei));
-                break;
             }
+            break;
+        }
 
-            case blenderType::TANH:
+        case blenderType::TANH:
+        {
+            forAll(yPlus, facei)
             {
-                const scalar phiTanh = tanh(pow4(0.1*yStar));
+                const scalar yPlusVisFace = yPlusVis(facei);
-                const scalar b1 = yPlusVis + yPlusLog;
+                const scalar yPlusLogFace = yPlusLog(facei);
+                const scalar b1 = yPlusVisFace + yPlusLogFace;
                 const scalar b2 =
-                    pow(pow(yPlusVis, 1.2) + pow(yPlusLog, 1.2), 1.0/1.2);
+                    pow
+                    (
+                        pow(yPlusVisFace, 1.2) + pow(yPlusLogFace, 1.2),
+                        1.0/1.2
+                    );
+                const scalar phiTanh = tanh(pow4(0.1*yStar(facei)));
 
                 yPlus[facei] = phiTanh*b1 + (1 - phiTanh)*b2;
-                break;
             }
+            break;
         }
     }