ENH: Updated incompressible omega wall functions to include averaging and to handle ACMI patches

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C

@@ -26,11 +26,11 @@ License
 #include "omegaWallFunctionFvPatchScalarField.H"
 #include "incompressible/turbulenceModel/turbulenceModel.H"
 #include "fvPatchFieldMapper.H"
+#include "fvMatrix.H"
 #include "volFields.H"
-#include "addToRunTimeSelectionTable.H"
 #include "wallFvPatch.H"
 #include "nutkWallFunctionFvPatchScalarField.H"
-
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,6 +64,196 @@ void omegaWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
 }
 
 
+void omegaWallFunctionFvPatchScalarField::setMaster()
+{
+    if (master_ != -1)
+    {
+        return;
+    }
+
+    const volScalarField& omega =
+        static_cast<const volScalarField&>(this->dimensionedInternalField());
+
+    const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
+
+    label master = -1;
+    forAll(bf, patchI)
+    {
+        if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchI]))
+        {
+            omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
+
+            if (master == -1)
+            {
+                master = patchI;
+            }
+
+            epf.master() = master;
+        }
+    }
+}
+
+
+void omegaWallFunctionFvPatchScalarField::createAveragingWeights()
+{
+    if (initialised_)
+    {
+        return;
+    }
+
+    const volScalarField& omega =
+        static_cast<const volScalarField&>(this->dimensionedInternalField());
+
+    const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
+
+    const fvMesh& mesh = omega.mesh();
+
+    volScalarField weights
+    (
+        IOobject
+        (
+            "weights",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false // do not register
+        ),
+        mesh,
+        dimensionedScalar("zero", dimless, 0.0)
+    );
+
+    DynamicList<label> omegaPatches(bf.size());
+    forAll(bf, patchI)
+    {
+        if (isA<omegaWallFunctionFvPatchScalarField>(bf[patchI]))
+        {
+            omegaPatches.append(patchI);
+
+            const labelUList& faceCells = bf[patchI].patch().faceCells();
+            forAll(faceCells, i)
+            {
+                label cellI = faceCells[i];
+                weights[cellI]++;
+            }
+        }
+    }
+
+    cornerWeights_.setSize(bf.size());
+    forAll(omegaPatches, i)
+    {
+        label patchI = omegaPatches[i];
+        const fvPatchField& wf = weights.boundaryField()[patchI];
+        cornerWeights_[patchI] = 1.0/wf.patchInternalField();
+    }
+
+    G_.setSize(dimensionedInternalField().size(), 0.0);
+    omega_.setSize(dimensionedInternalField().size(), 0.0);
+
+    initialised_ = true;
+}
+
+
+omegaWallFunctionFvPatchScalarField&
+omegaWallFunctionFvPatchScalarField::omegaPatch(const label patchI)
+{
+    const volScalarField& omega =
+        static_cast<const volScalarField&>(this->dimensionedInternalField());
+
+    const volScalarField::GeometricBoundaryField& bf = omega.boundaryField();
+
+    const omegaWallFunctionFvPatchScalarField& epf =
+        refCast<const omegaWallFunctionFvPatchScalarField>(bf[patchI]);
+
+    return const_cast<omegaWallFunctionFvPatchScalarField&>(epf);
+}
+
+
+void omegaWallFunctionFvPatchScalarField::calculateTurbulenceFields
+(
+    const turbulenceModel& turbulence,
+    scalarField& G0,
+    scalarField& omega0
+)
+{
+    // accumulate all of the G and omega contributions
+    forAll(cornerWeights_, patchI)
+    {
+        if (!cornerWeights_[patchI].empty())
+        {
+            omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
+
+            const List<scalar>& w = cornerWeights_[patchI];
+
+            epf.calculate(turbulence, w, epf.patch(), G0, omega0);
+        }
+    }
+
+    // apply zero-gradient condition for omega
+    forAll(cornerWeights_, patchI)
+    {
+        if (!cornerWeights_[patchI].empty())
+        {
+            omegaWallFunctionFvPatchScalarField& epf = omegaPatch(patchI);
+
+            epf == scalarField(omega0, epf.patch().faceCells());
+        }
+    }
+}
+
+
+void omegaWallFunctionFvPatchScalarField::calculate
+(
+    const turbulenceModel& turbulence,
+    const List<scalar>& cornerWeights,
+    const fvPatch& patch,
+    scalarField& G,
+    scalarField& omega
+)
+{
+    const label patchI = patch.index();
+
+    const scalarField& y = turbulence.y()[patchI];
+
+    const scalar Cmu25 = pow025(Cmu_);
+
+    const tmp<volScalarField> tk = turbulence.k();
+    const volScalarField& k = tk();
+
+    const tmp<volScalarField> tnu = turbulence.nu();
+    const scalarField& nuw = tnu().boundaryField()[patchI];
+
+    const tmp<volScalarField> tnut = turbulence.nut();
+    const volScalarField& nut = tnut();
+    const scalarField& nutw = nut.boundaryField()[patchI];
+
+    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
+
+    const scalarField magGradUw(mag(Uw.snGrad()));
+
+    // Set omega and G
+    forAll(nutw, faceI)
+    {
+        label cellI = patch.faceCells()[faceI];
+
+        scalar w = cornerWeights[faceI];
+
+        scalar omegaVis = 6.0*nuw[faceI]/(beta1_*sqr(y[faceI]));
+
+        scalar omegaLog = sqrt(k[cellI])/(Cmu25*kappa_*y[faceI]);
+
+        omega[cellI] = w*sqrt(sqr(omegaVis) + sqr(omegaLog));
+
+        G[cellI] =
+            w
+           *(nutw[faceI] + nuw[faceI])
+           *magGradUw[faceI]
+           *Cmu25*sqrt(k[cellI])
+           /(kappa_*y[faceI]);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
@@ -72,12 +262,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedInternalValueFvPatchField<scalar>(p, iF),
+    fixedValueFvPatchField<scalar>(p, iF),
     Cmu_(0.09),
     kappa_(0.41),
     E_(9.8),
     beta1_(0.075),
-    yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_))
+    yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -91,12 +286,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    fixedInternalValueFvPatchField<scalar>(ptf, p, iF, mapper),
+    fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
     Cmu_(ptf.Cmu_),
     kappa_(ptf.kappa_),
     E_(ptf.E_),
     beta1_(ptf.beta1_),
-    yPlusLam_(ptf.yPlusLam_)
+    yPlusLam_(ptf.yPlusLam_),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -109,12 +309,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     const dictionary& dict
 )
 :
-    fixedInternalValueFvPatchField<scalar>(p, iF, dict),
+    fixedValueFvPatchField<scalar>(p, iF, dict),
     Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
     E_(dict.lookupOrDefault<scalar>("E", 9.8)),
     beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075)),
-    yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_))
+    yPlusLam_(nutkWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -125,12 +330,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     const omegaWallFunctionFvPatchScalarField& owfpsf
 )
 :
-    fixedInternalValueFvPatchField<scalar>(owfpsf),
+    fixedValueFvPatchField<scalar>(owfpsf),
     Cmu_(owfpsf.Cmu_),
     kappa_(owfpsf.kappa_),
     E_(owfpsf.E_),
     beta1_(owfpsf.beta1_),
-    yPlusLam_(owfpsf.yPlusLam_)
+    yPlusLam_(owfpsf.yPlusLam_),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -142,12 +352,17 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedInternalValueFvPatchField<scalar>(owfpsf, iF),
+    fixedValueFvPatchField<scalar>(owfpsf, iF),
     Cmu_(owfpsf.Cmu_),
     kappa_(owfpsf.kappa_),
     E_(owfpsf.E_),
     beta1_(owfpsf.beta1_),
-    yPlusLam_(owfpsf.yPlusLam_)
+    yPlusLam_(owfpsf.yPlusLam_),
+    G_(),
+    omega_(),
+    initialised_(false),
+    master_(-1),
+    cornerWeights_()
 {
     checkType();
 }
@@ -155,6 +370,38 @@ omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+scalarField& omegaWallFunctionFvPatchScalarField::G(bool init)
+{
+    if (patch().index() == master_)
+    {
+        if (init)
+        {
+            G_ = 0.0;
+        }
+
+        return G_;
+    }
+
+    return omegaPatch(master_).G();
+}
+
+
+scalarField& omegaWallFunctionFvPatchScalarField::omega(bool init)
+{
+    if (patch().index() == master_)
+    {
+        if (init)
+        {
+            omega_ = 0.0;
+        }
+
+        return omega_;
+    }
+
+    return omegaPatch(master_).omega(init);
+}
+
+
 void omegaWallFunctionFvPatchScalarField::updateCoeffs()
 {
     if (updated())
@@ -162,70 +409,168 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const label patchI = patch().index();
-
     const turbulenceModel& turbulence =
-        db().lookupObject<turbulenceModel>("turbulenceModel");
-    const scalarField& y = turbulence.y()[patchI];
+        db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
 
-    const scalar Cmu25 = pow025(Cmu_);
+    setMaster();
+
+    if (patch().index() == master_)
+    {
+        createAveragingWeights();
+        calculateTurbulenceFields(turbulence, G(true), omega(true));
+    }
+
+    const scalarField& G0 = this->G();
+    const scalarField& omega0 = this->omega();
+
+    typedef DimensionedField<scalar, volMesh> FieldType;
 
-    volScalarField& G =
-        const_cast<volScalarField&>
+    FieldType& G =
+        const_cast<FieldType&>
         (
-            db().lookupObject<volScalarField>
-            (
-                turbulence.GName()
-            )
+            db().lookupObject<FieldType>(turbulence.GName())
         );
 
-    DimensionedField<scalar, volMesh>& omega =
-        const_cast<DimensionedField<scalar, volMesh>&>
+    FieldType& omega = const_cast<FieldType&>(dimensionedInternalField());
+
+    forAll(*this, faceI)
+    {
+        label cellI = patch().faceCells()[faceI];
+
+        G[cellI] = G0[cellI];
+        omega[cellI] = omega0[cellI];
+    }
+
+    fvPatchField<scalar>::updateCoeffs();
+}
+
+
+void omegaWallFunctionFvPatchScalarField::updateCoeffs
+(
+    const scalarField& weights
+)
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const turbulenceModel& turbulence =
+        db().lookupObject<turbulenceModel>(turbulenceModel::typeName);
+
+    setMaster();
+
+    if (patch().index() == master_)
+    {
+        createAveragingWeights();
+        calculateTurbulenceFields(turbulence, G(true), omega(true));
+    }
+
+    const scalarField& G0 = this->G();
+    const scalarField& omega0 = this->omega();
+
+    typedef DimensionedField<scalar, volMesh> FieldType;
+
+    FieldType& G =
+        const_cast<FieldType&>
         (
-            dimensionedInternalField()
+            db().lookupObject<FieldType>(turbulence.GName())
         );
 
-    const tmp<volScalarField> tk = turbulence.k();
-    const volScalarField& k = tk();
+    FieldType& omega = const_cast<FieldType&>(dimensionedInternalField());
 
-    const tmp<volScalarField> tnu = turbulence.nu();
-    const scalarField& nuw = tnu().boundaryField()[patchI];
+    // only set the values if the weights are < 1 - tolerance
+    forAll(weights, faceI)
+    {
+        scalar w = weights[faceI];
 
-    const tmp<volScalarField> tnut = turbulence.nut();
-    const volScalarField& nut = tnut();
-    const scalarField& nutw = nut.boundaryField()[patchI];
+        if (w < 1.0 - 1e-6)
+        {
+            label cellI = patch().faceCells()[faceI];
 
-    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
+            G[cellI] = w*G[cellI] + (1.0 - w)*G0[cellI];
+            omega[cellI] = w*omega[cellI] + (1.0 - w)*omega0[cellI];
+        }
+    }
 
-    const scalarField magGradUw(mag(Uw.snGrad()));
+    fvPatchField<scalar>::updateCoeffs();
+}
 
-    // Set omega and G
-    forAll(nutw, faceI)
+
+void omegaWallFunctionFvPatchScalarField::manipulateMatrix
+(
+    fvMatrix<scalar>& matrix
+)
+{
+    if (manipulatedMatrix())
     {
-        label faceCellI = patch().faceCells()[faceI];
+        return;
+    }
 
-        scalar omegaVis = 6.0*nuw[faceI]/(beta1_*sqr(y[faceI]));
+    matrix.setValues(patch().faceCells(), patchInternalField());
 
-        scalar omegaLog = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]);
+    fvPatchField<scalar>::manipulateMatrix(matrix);
+}
 
-        omega[faceCellI] = sqrt(sqr(omegaVis) + sqr(omegaLog));
 
-        G[faceCellI] =
-            (nutw[faceI] + nuw[faceI])
-           *magGradUw[faceI]
-           *Cmu25*sqrt(k[faceCellI])
-           /(kappa_*y[faceI]);
+void omegaWallFunctionFvPatchScalarField::manipulateMatrix
+(
+    fvMatrix<scalar>& matrix,
+    const Field<scalar>& weights
+)
+{
+    if (manipulatedMatrix())
+    {
+        return;
+    }
+
+    // filter weights so that we only apply the constraint where the
+    // weight > SMALL
+    DynamicList<label> constraintCells(weights.size());
+    DynamicList<scalar> constraintomega(weights.size());
+    const labelUList& faceCells = patch().faceCells();
+
+    const DimensionedField<scalar, volMesh>& omega
+        = dimensionedInternalField();
+
+    label nConstrainedCells = 0;
+
+
+    forAll(weights, faceI)
+    {
+        // only set the values if the weights are < 1 - tolerance
+        if (weights[faceI] < (1.0 - 1e-6))
+        {
+            nConstrainedCells++;
+
+            label cellI = faceCells[faceI];
+
+            constraintCells.append(cellI);
+            constraintomega.append(omega[cellI]);
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "Patch: " << patch().name()
+            << ": number of constrained cells = " << nConstrainedCells
+            << " out of " << patch().size()
+            << endl;
     }
 
-    fixedInternalValueFvPatchField<scalar>::updateCoeffs();
+    matrix.setValues
+    (
+        constraintCells,
+        scalarField(constraintomega.xfer())
+    );
 
-    // TODO: perform averaging for cells sharing more than one boundary face
+    fvPatchField<scalar>::manipulateMatrix(matrix);
 }
 
 
 void omegaWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
-    fixedInternalValueFvPatchField<scalar>::write(os);
+    fixedValueFvPatchField<scalar>::write(os);
     writeLocalEntries(os);
     writeEntry("value", os);
 }

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,7 +76,7 @@ SourceFiles
 #ifndef omegaWallFunctionFvPatchScalarField_H
 #define omegaWallFunctionFvPatchScalarField_H
 
-#include "fixedInternalValueFvPatchField.H"
+#include "fixedValueFvPatchField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -85,13 +85,15 @@ namespace Foam
 namespace incompressible
 {
 
+class turbulenceModel;
+
 /*---------------------------------------------------------------------------*\
-           Class omegaWallFunctionFvPatchScalarField Declaration
+             Class omegaWallFunctionFvPatchScalarField Declaration
 \*---------------------------------------------------------------------------*/
 
 class omegaWallFunctionFvPatchScalarField
 :
-    public fixedInternalValueFvPatchField<scalar>
+    public fixedValueFvPatchField<scalar>
 {
 protected:
 
@@ -112,6 +114,21 @@ protected:
         //- Y+ at the edge of the laminar sublayer
         scalar yPlusLam_;
 
+        //- Local copy of turbulence G field
+        scalarField G_;
+
+        //- Local copy of turbulence omega field
+        scalarField omega_;
+
+        //- Initialised flag
+        bool initialised_;
+
+        //- Master patch ID
+        label master_;
+
+        //- List of averaging corner weights
+        List<List<scalar> > cornerWeights_;
+
 
     // Protected Member Functions
 
@@ -121,6 +138,44 @@ protected:
         //- Write local wall function variables
         virtual void writeLocalEntries(Ostream&) const;
 
+        //- Set the master patch - master is responsible for updating all
+        //  wall function patches
+        virtual void setMaster();
+
+        //- Create the averaging weights for cells which are bounded by
+        //  multiple wall function faces
+        virtual void createAveragingWeights();
+
+        //- Helper function to return non-const access to an omega patch
+        virtual omegaWallFunctionFvPatchScalarField& omegaPatch
+        (
+            const label patchI
+        );
+
+        //- Main driver to calculate the turbulence fields
+        virtual void calculateTurbulenceFields
+        (
+            const turbulenceModel& turbulence,
+            scalarField& G0,
+            scalarField& omega0
+        );
+
+        //- Calculate the omega and G
+        virtual void calculate
+        (
+            const turbulenceModel& turbulence,
+            const List<scalar>& cornerWeights,
+            const fvPatch& patch,
+            scalarField& G,
+            scalarField& omega
+        );
+
+        //- Return non-const access to the master patch ID
+        virtual label& master()
+        {
+            return master_;
+        }
+
 
 public:
 
@@ -193,11 +248,33 @@ public:
 
     // Member functions
 
+        // Access
+
+            //- Return non-const access to the master's G field
+            scalarField& G(bool init = false);
+
+            //- Return non-const access to the master's omega field
+            scalarField& omega(bool init = false);
+
+
         // Evaluation functions
 
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs(const scalarField& weights);
+
+            //- Manipulate matrix
+            virtual void manipulateMatrix(fvMatrix<scalar>& matrix);
+
+            //- Manipulate matrix with given weights
+            virtual void manipulateMatrix
+            (
+                fvMatrix<scalar>& matrix,
+                const scalarField& weights
+            );
+
 
         // I-O