Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options

@@ -2,7 +2,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModel/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \
     -IphaseIncompressibleTurbulenceModels/lnInclude \

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/Make/options

@@ -3,7 +3,7 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/transportModel \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModel/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \
     -I../twoPhaseSystem/lnInclude \

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C

@@ -24,14 +24,12 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "kineticTheoryModel.H"
-#include "surfaceInterpolate.H"
 #include "mathematicalConstants.H"
 #include "twoPhaseSystem.H"
-#include "fvcDiv.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::kineticTheoryModel::kineticTheoryModel
+Foam::RASModels::kineticTheoryModel::kineticTheoryModel
 (
     const volScalarField& alpha,
     const geometricOneField& rho,
@@ -43,7 +41,7 @@ Foam::kineticTheoryModel::kineticTheoryModel
     const word& type
 )
 :
-    RASModels::eddyViscosity<PhaseIncompressibleTurbulenceModel<phaseModel> >
+    eddyViscosity<RASModel<PhaseIncompressibleTurbulenceModel<phaseModel> > >
     (
         type,
         alpha,
@@ -131,18 +129,224 @@ Foam::kineticTheoryModel::kineticTheoryModel
         U.mesh(),
         dimensionedScalar("zero", dimensionSet(0, 2, -1, 0, 0), 0.0)
     )
-{}
+{
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::kineticTheoryModel::~kineticTheoryModel()
+Foam::RASModels::kineticTheoryModel::~kineticTheoryModel()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::kineticTheoryModel::correct()
+bool Foam::RASModels::kineticTheoryModel::read()
+{
+    if
+    (
+        eddyViscosity
+        <
+            RASModel<PhaseIncompressibleTurbulenceModel<phaseModel> >
+        >::read()
+    )
+    {
+        this->coeffDict().lookup("equilibrium") >> equilibrium_;
+        e_.readIfPresent(this->coeffDict());
+        alphaMax_.readIfPresent(this->coeffDict());
+        alphaMinFriction_.readIfPresent(this->coeffDict());
+
+        viscosityModel_->read();
+        conductivityModel_->read();
+        radialModel_->read();
+        granularPressureModel_->read();
+        frictionalStressModel_->read();
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::RASModels::kineticTheoryModel::k() const
+{
+    notImplemented("kineticTheoryModel::k()");
+    return nut_;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::RASModels::kineticTheoryModel::epsilon() const
+{
+    notImplemented("kineticTheoryModel::epsilon()");
+    return nut_;
+}
+
+
+Foam::tmp<Foam::volSymmTensorField>
+Foam::RASModels::kineticTheoryModel::R() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                IOobject::groupName("R", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+          - (this->nut_)*dev(twoSymm(fvc::grad(this->U_)))
+          - (lambda_*fvc::div(this->phi_))*symmTensor::I
+        )
+    );
+}
+
+
+/*
+Foam::tmp<Foam::volScalarField>
+Foam::RASModels::kineticTheoryModel::pp() const
+{
+
+    // Particle pressure coefficient
+    // Coefficient in front of Theta (Eq. 3.22, p. 45)
+    volScalarField PsCoeff
+    (
+        granularPressureModel_->granularPressureCoeff
+        (
+            alpha,
+            gs0,
+            rho,
+            e_
+        )
+    );
+
+    // Frictional pressure
+    volScalarField pf
+    (
+        frictionalStressModel_->frictionalPressure
+        (
+            alpha,
+            alphaMinFriction_,
+            alphaMax_
+        )
+    );
+
+    // Return total particle pressure
+    return PsCoeff*Theta_ + pf;
+}
+*/
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::RASModels::kineticTheoryModel::pPrime() const
+{
+    // Local references
+    const volScalarField& alpha = this->alpha_;
+    const volScalarField& rho = phase_.rho();
+
+    return
+    (
+        Theta_
+       *granularPressureModel_->granularPressureCoeffPrime
+        (
+            alpha,
+            radialModel_->g0(alpha, alphaMinFriction_, alphaMax_),
+            radialModel_->g0prime(alpha, alphaMinFriction_, alphaMax_),
+            rho,
+            e_
+        )
+     +  frictionalStressModel_->frictionalPressurePrime
+        (
+            alpha,
+            alphaMinFriction_,
+            alphaMax_
+        )
+    );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::RASModels::kineticTheoryModel::pPrimef() const
+{
+    // Local references
+    const volScalarField& alpha = this->alpha_;
+    const volScalarField& rho = phase_.rho();
+
+    return fvc::interpolate
+    (
+        Theta_
+       *granularPressureModel_->granularPressureCoeffPrime
+        (
+            alpha,
+            radialModel_->g0(alpha, alphaMinFriction_, alphaMax_),
+            radialModel_->g0prime(alpha, alphaMinFriction_, alphaMax_),
+            rho,
+            e_
+        )
+     +  frictionalStressModel_->frictionalPressurePrime
+        (
+            alpha,
+            alphaMinFriction_,
+            alphaMax_
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volSymmTensorField>
+Foam::RASModels::kineticTheoryModel::devRhoReff() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                IOobject::groupName("devRhoReff", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+          - (this->rho_*this->nut_)
+           *dev(twoSymm(fvc::grad(this->U_)))
+          - ((this->rho_*lambda_)*fvc::div(this->phi_))*symmTensor::I
+        )
+    );
+}
+
+
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::RASModels::kineticTheoryModel::divDevRhoReff
+(
+    volVectorField& U
+) const
+{
+    return
+    (
+      - fvm::laplacian(this->rho_*this->nut_, U)
+      - fvc::div
+        (
+            (this->rho_*this->nut_)*dev2(T(fvc::grad(U)))
+          + ((this->rho_*lambda_)*fvc::div(this->phi_))
+           *dimensioned<symmTensor>("I", dimless, symmTensor::I)
+        )
+    );
+}
+
+
+void Foam::RASModels::kineticTheoryModel::correct()
 {
     // Local references
     volScalarField alpha(max(this->alpha_, 0.0));
@@ -345,202 +549,4 @@ void Foam::kineticTheoryModel::correct()
     }
 }
 
-
-/*
-Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::pp() const
-{
-
-    // Particle pressure coefficient
-    // Coefficient in front of Theta (Eq. 3.22, p. 45)
-    volScalarField PsCoeff
-    (
-        granularPressureModel_->granularPressureCoeff
-        (
-            alpha,
-            gs0,
-            rho,
-            e_
-        )
-    );
-
-    // Frictional pressure
-    volScalarField pf
-    (
-        frictionalStressModel_->frictionalPressure
-        (
-            alpha,
-            alphaMinFriction_,
-            alphaMax_
-        )
-    );
-
-    // Return total particle pressure
-    return PsCoeff*Theta_ + pf;
-}
-*/
-
-
-Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::pPrime() const
-{
-    // Local references
-    const volScalarField& alpha = this->alpha_;
-    const volScalarField& rho = phase_.rho();
-
-    return
-    (
-        Theta_
-       *granularPressureModel_->granularPressureCoeffPrime
-        (
-            alpha,
-            radialModel_->g0(alpha, alphaMinFriction_, alphaMax_),
-            radialModel_->g0prime(alpha, alphaMinFriction_, alphaMax_),
-            rho,
-            e_
-        )
-     +  frictionalStressModel_->frictionalPressurePrime
-        (
-            alpha,
-            alphaMinFriction_,
-            alphaMax_
-        )
-    );
-}
-
-
-Foam::tmp<Foam::surfaceScalarField> Foam::kineticTheoryModel::pPrimef() const
-{
-    // Local references
-    const volScalarField& alpha = this->alpha_;
-    const volScalarField& rho = phase_.rho();
-
-    return fvc::interpolate
-    (
-        Theta_
-       *granularPressureModel_->granularPressureCoeffPrime
-        (
-            alpha,
-            radialModel_->g0(alpha, alphaMinFriction_, alphaMax_),
-            radialModel_->g0prime(alpha, alphaMinFriction_, alphaMax_),
-            rho,
-            e_
-        )
-     +  frictionalStressModel_->frictionalPressurePrime
-        (
-            alpha,
-            alphaMinFriction_,
-            alphaMax_
-        )
-    );
-}
-
-
-void Foam::kineticTheoryModel::correctNut()
-{}
-
-
-Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::k() const
-{
-    notImplemented("kineticTheoryModel::k()");
-    return nut_;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::epsilon() const
-{
-    notImplemented("kineticTheoryModel::epsilon()");
-    return nut_;
-}
-
-
-Foam::tmp<Foam::volSymmTensorField> Foam::kineticTheoryModel::R() const
-{
-    return tmp<volSymmTensorField>
-    (
-        new volSymmTensorField
-        (
-            IOobject
-            (
-                IOobject::groupName("R", this->U_.group()),
-                this->runTime_.timeName(),
-                this->mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-          - (this->nut_)*dev(twoSymm(fvc::grad(this->U_)))
-          - (lambda_*fvc::div(this->phi_))*symmTensor::I
-        )
-    );
-}
-
-
-Foam::tmp<Foam::volSymmTensorField> Foam::kineticTheoryModel::devRhoReff() const
-{
-    return tmp<volSymmTensorField>
-    (
-        new volSymmTensorField
-        (
-            IOobject
-            (
-                IOobject::groupName("devRhoReff", this->U_.group()),
-                this->runTime_.timeName(),
-                this->mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-          - (this->rho_*this->nut_)
-           *dev(twoSymm(fvc::grad(this->U_)))
-          - ((this->rho_*lambda_)*fvc::div(this->phi_))*symmTensor::I
-        )
-    );
-}
-
-
-Foam::tmp<Foam::fvVectorMatrix> Foam::kineticTheoryModel::divDevRhoReff
-(
-    volVectorField& U
-) const
-{
-    return
-    (
-      - fvm::laplacian(this->rho_*this->nut_, U)
-      - fvc::div
-        (
-            (this->rho_*this->nut_)*dev2(T(fvc::grad(U)))
-          + ((this->rho_*lambda_)*fvc::div(this->phi_))
-           *dimensioned<symmTensor>("I", dimless, symmTensor::I)
-        )
-    );
-}
-
-
-bool Foam::kineticTheoryModel::read()
-{
-    if
-    (
-        RASModels::eddyViscosity
-        <
-            PhaseIncompressibleTurbulenceModel<phaseModel>
-        >::read()
-    )
-    {
-        this->coeffDict().lookup("equilibrium") >> equilibrium_;
-        e_.readIfPresent(this->coeffDict());
-        alphaMax_.readIfPresent(this->coeffDict());
-        alphaMinFriction_.readIfPresent(this->coeffDict());
-
-        viscosityModel_->read();
-        conductivityModel_->read();
-        radialModel_->read();
-        granularPressureModel_->read();
-        frictionalStressModel_->read();
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
 // ************************************************************************* //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H

@@ -25,6 +25,17 @@ Class
     Foam::kineticTheoryModel
 
 Description
+    Kinetic theory particle phase RAS model
+
+    Reference:
+    \verbatim
+        "Derivation, implementation, and validation of computer simulation
+        models for gas-solid fluidized beds",
+        B.G.M. van Wachem,
+        Ph.D. Thesis, Delft University of Technology, Amsterdam, 2000.
+    \endverbatim
+
+    There are no default model coefficients.
 
 SourceFiles
     kineticTheoryModel.C
@@ -34,11 +45,11 @@ SourceFiles
 #ifndef kineticTheoryModel_H
 #define kineticTheoryModel_H
 
+#include "RASModel.H"
 #include "eddyViscosity.H"
 #include "PhaseIncompressibleTurbulenceModel.H"
-#include "dragModel.H"
 #include "phaseModel.H"
-#include "autoPtr.H"
+#include "dragModel.H"
 #include "viscosityModel.H"
 #include "conductivityModel.H"
 #include "radialModel.H"
@@ -50,6 +61,8 @@ SourceFiles
 
 namespace Foam
 {
+namespace RASModels
+{
 
 /*---------------------------------------------------------------------------*\
                            Class kineticTheoryModel Declaration
@@ -57,9 +70,9 @@ namespace Foam
 
 class kineticTheoryModel
 :
-    public RASModels::eddyViscosity
+    public eddyViscosity
     <
-        PhaseIncompressibleTurbulenceModel<phaseModel>
+        RASModel<PhaseIncompressibleTurbulenceModel<phaseModel> >
     >
 {
     // Private data
@@ -118,6 +131,9 @@ class kineticTheoryModel
 
     // Private Member Functions
 
+        void correctNut()
+        {}
+
         //- Disallow default bitwise copy construct
         kineticTheoryModel(const kineticTheoryModel&);
 
@@ -125,13 +141,6 @@ class kineticTheoryModel
         void operator=(const kineticTheoryModel&);
 
 
-protected:
-
-    // Protected member functions
-
-        virtual void correctNut();
-
-
 public:
 
     //- Runtime type information
@@ -160,6 +169,9 @@ public:
 
     // Member Functions
 
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
         //- Return the effective viscosity
         virtual tmp<volScalarField> nuEff() const
         {
@@ -197,14 +209,12 @@ public:
 
         //- Solve the kinetic theory equations and correct the viscosity
         virtual void correct();
-
-        //- Re-read model coefficients if they have changed
-        virtual bool read();
 };
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+} // End namespace RASModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phaseIncompressibleTurbulenceModels.C

@@ -24,149 +24,64 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "PhaseIncompressibleTurbulenceModel.H"
-#include "laminar.H"
-#include "RASModel.H"
-#include "kEpsilon.H"
-#include "LaheyKEpsilon.H"
-#include "continuousGasKEpsilon.H"
-#include "kineticTheoryModel.H"
-#include "phasePressureModel.H"
 #include "phaseModel.H"
 #include "addToRunTimeSelectionTable.H"
+#include "makeTurbulenceModel.H"
 
-namespace Foam
+#include "laminar.H"
-{
+#include "RASModel.H"
-    typedef TurbulenceModel
+#include "LESModel.H"
-    <
-        volScalarField,
-        geometricOneField,
-        incompressibleTurbulenceModel,
-        phaseModel
-    > basePhaseIncompressibleTransportTurbulenceModel;
 
-    defineTemplateRunTimeSelectionTable
+makeBaseTurbulenceModel
-    (
+(
-        basePhaseIncompressibleTransportTurbulenceModel,
+    volScalarField,
-        dictionary
+    geometricOneField,
-    );
+    incompressibleTurbulenceModel,
+    PhaseIncompressibleTurbulenceModel,
+    phaseModel
+);
 
-    typedef PhaseIncompressibleTurbulenceModel<phaseModel>
+#define makeRASModel(Type)                                                     \
-        incompressibleTransportTurbulenceModel;
+    makeTemplatedTurbulenceModel                                               \
+    (phaseModelPhaseIncompressibleTurbulenceModel, RAS, Type)
 
-    typedef laminar<incompressibleTransportTurbulenceModel>
+#define makeLESModel(Type)                                                     \
-        incompressibleLaminar;
+    makeTemplatedTurbulenceModel                                               \
+    (phaseModelPhaseIncompressibleTurbulenceModel, LES, Type)
 
-    defineNamedTemplateTypeNameAndDebug(incompressibleLaminar, 0);
+#include "kEpsilon.H"
+makeRASModel(kEpsilon);
 
-    addToRunTimeSelectionTable
+#include "LaheyKEpsilon.H"
-    (
+makeRASModel(LaheyKEpsilon);
-        basePhaseIncompressibleTransportTurbulenceModel,
+
-        incompressibleLaminar,
+#include "continuousGasKEpsilon.H"
-        dictionary
+makeRASModel(continuousGasKEpsilon);
-    );
+
+#include "Smagorinsky.H"
+makeLESModel(Smagorinsky);
+
+#include "kEqn.H"
+makeLESModel(kEqn);
+
+#include "SmagorinskyZhang.H"
+makeLESModel(SmagorinskyZhang);
+
+#include "NicenoKEqn.H"
+makeLESModel(NicenoKEqn);
+
+#include "continuousGasKEqn.H"
+makeLESModel(continuousGasKEqn);
 
 
-    typedef RASModel<incompressibleTransportTurbulenceModel>
+#include "kineticTheoryModel.H"
-        incompressibleRASModel;
+makeTurbulenceModel
-
+(phaseModelPhaseIncompressibleTurbulenceModel, RAS, kineticTheoryModel);
-    defineNamedTemplateTypeNameAndDebug(incompressibleRASModel, 0);
-
-    defineTemplateRunTimeSelectionTable(incompressibleRASModel, dictionary);
-
-    addToRunTimeSelectionTable
-    (
-        basePhaseIncompressibleTransportTurbulenceModel,
-        incompressibleRASModel,
-        dictionary
-    );
-
-    namespace RASModels
-    {
-        typedef kEpsilon<incompressibleTransportTurbulenceModel>
-            incompressiblekEpsilon;
-
-        defineNamedTemplateTypeNameAndDebug(incompressiblekEpsilon, 0);
-
-        addToRunTimeSelectionTable
-        (
-            incompressibleRASModel,
-            incompressiblekEpsilon,
-            dictionary
-        );
-    }
-
-    namespace RASModels
-    {
-        typedef LaheyKEpsilon<incompressibleTransportTurbulenceModel>
-            incompressibleLaheyKEpsilon;
-
-        defineNamedTemplateTypeNameAndDebug(incompressibleLaheyKEpsilon, 0);
-
-        addToRunTimeSelectionTable
-        (
-            incompressibleRASModel,
-            incompressibleLaheyKEpsilon,
-            dictionary
-        );
-    }
-
-    namespace RASModels
-    {
-        typedef continuousGasKEpsilon<incompressibleTransportTurbulenceModel>
-            incompressiblecontinuousGasKEpsilon;
-
-        defineNamedTemplateTypeNameAndDebug
-        (
-            incompressiblecontinuousGasKEpsilon,
-            0
-        );
-
-        addToRunTimeSelectionTable
-        (
-            incompressibleRASModel,
-            incompressiblecontinuousGasKEpsilon,
-            dictionary
-        );
-    }
-}
 
 
-namespace Foam
+#include "phasePressureModel.H"
-{
+makeTurbulenceModel
-    typedef PhaseIncompressibleTurbulenceModel<phaseModel>
+(phaseModelPhaseIncompressibleTurbulenceModel, RAS, phasePressureModel);
-        incompressibleTransportTurbulenceModel;
-
-    typedef RASModel<incompressibleTransportTurbulenceModel>
-        incompressibleRASModel;
-
-    defineTypeNameAndDebug(kineticTheoryModel, 0);
-
-    addToRunTimeSelectionTable
-    (
-        incompressibleRASModel,
-        kineticTheoryModel,
-        dictionary
-    );
-}
-
-
-namespace Foam
-{
-    typedef PhaseIncompressibleTurbulenceModel<phaseModel>
-        incompressibleTransportTurbulenceModel;
-
-    typedef RASModel<incompressibleTransportTurbulenceModel>
-        incompressibleRASModel;
-
-    defineTypeNameAndDebug(phasePressureModel, 0);
-
-    addToRunTimeSelectionTable
-    (
-        incompressibleRASModel,
-        phasePressureModel,
-        dictionary
-    );
-}
 
 
 // ************************************************************************* //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C

@@ -24,13 +24,11 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "phasePressureModel.H"
-#include "surfaceInterpolate.H"
 #include "twoPhaseSystem.H"
-#include "dimensionedType.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::phasePressureModel::phasePressureModel
+Foam::RASModels::phasePressureModel::phasePressureModel
 (
     const volScalarField& alpha,
     const geometricOneField& rho,
@@ -42,7 +40,7 @@ Foam::phasePressureModel::phasePressureModel
     const word& type
 )
 :
-    RASModels::eddyViscosity<PhaseIncompressibleTurbulenceModel<phaseModel> >
+    eddyViscosity<RASModel<PhaseIncompressibleTurbulenceModel<phaseModel> > >
     (
         type,
         alpha,
@@ -67,64 +65,64 @@ Foam::phasePressureModel::phasePressureModel
     )
 {
     this->nut_ == dimensionedScalar("zero", this->nut_.dimensions(), 0.0);
+
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
 }
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::phasePressureModel::~phasePressureModel()
+Foam::RASModels::phasePressureModel::~phasePressureModel()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::phasePressureModel::correct()
+bool Foam::RASModels::phasePressureModel::read()
-{}
-
-
-void Foam::phasePressureModel::correctNut()
-{}
-
-
-Foam::tmp<Foam::volScalarField> Foam::phasePressureModel::pPrime() const
 {
-    return
+    if
-        g0_
+    (
-       *min
+        eddyViscosity
-        (
+        <
-            exp(preAlphaExp_*(this->alpha_ - alphaMax_)),
+            RASModel<PhaseIncompressibleTurbulenceModel<phaseModel> >
-            expMax_
+        >::read()
-        );
+    )
+    {
+        this->coeffDict().lookup("alphaMax") >> alphaMax_;
+        this->coeffDict().lookup("preAlphaExp") >> preAlphaExp_;
+        this->coeffDict().lookup("expMax") >> expMax_;
+        g0_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
 }
 
 
-Foam::tmp<Foam::surfaceScalarField> Foam::phasePressureModel::pPrimef() const
+Foam::tmp<Foam::volScalarField>
-{
+Foam::RASModels::phasePressureModel::k() const
-    return
-        g0_
-       *min
-        (
-            exp(preAlphaExp_*(fvc::interpolate(this->alpha_) - alphaMax_)),
-            expMax_
-        );
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::phasePressureModel::k() const
 {
     notImplemented("phasePressureModel::k()");
     return nut_;
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::phasePressureModel::epsilon() const
+Foam::tmp<Foam::volScalarField>
+Foam::RASModels::phasePressureModel::epsilon() const
 {
     notImplemented("phasePressureModel::epsilon()");
     return nut_;
 }
 
 
-Foam::tmp<Foam::volSymmTensorField> Foam::phasePressureModel::R() const
+Foam::tmp<Foam::volSymmTensorField>
+Foam::RASModels::phasePressureModel::R() const
 {
     return tmp<volSymmTensorField>
     (
@@ -150,7 +148,34 @@ Foam::tmp<Foam::volSymmTensorField> Foam::phasePressureModel::R() const
 }
 
 
-Foam::tmp<Foam::volSymmTensorField> Foam::phasePressureModel::devRhoReff() const
+Foam::tmp<Foam::volScalarField>
+Foam::RASModels::phasePressureModel::pPrime() const
+{
+    return
+        g0_
+       *min
+        (
+            exp(preAlphaExp_*(this->alpha_ - alphaMax_)),
+            expMax_
+        );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField>
+Foam::RASModels::phasePressureModel::pPrimef() const
+{
+    return
+        g0_
+       *min
+        (
+            exp(preAlphaExp_*(fvc::interpolate(this->alpha_) - alphaMax_)),
+            expMax_
+        );
+}
+
+
+Foam::tmp<Foam::volSymmTensorField>
+Foam::RASModels::phasePressureModel::devRhoReff() const
 {
     return tmp<volSymmTensorField>
     (
@@ -176,7 +201,8 @@ Foam::tmp<Foam::volSymmTensorField> Foam::phasePressureModel::devRhoReff() const
 }
 
 
-Foam::tmp<Foam::fvVectorMatrix> Foam::phasePressureModel::divDevRhoReff
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::RASModels::phasePressureModel::divDevRhoReff
 (
     volVectorField& U
 ) const
@@ -192,28 +218,8 @@ Foam::tmp<Foam::fvVectorMatrix> Foam::phasePressureModel::divDevRhoReff
 }
 
 
-bool Foam::phasePressureModel::read()
+void Foam::RASModels::phasePressureModel::correct()
-{
+{}
-    if
-    (
-        RASModels::eddyViscosity
-        <
-            PhaseIncompressibleTurbulenceModel<phaseModel>
-        >::read()
-    )
-    {
-        this->coeffDict().lookup("alphaMax") >> alphaMax_;
-        this->coeffDict().lookup("preAlphaExp") >> preAlphaExp_;
-        this->coeffDict().lookup("expMax") >> expMax_;
-        g0_.readIfPresent(this->coeffDict());
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
 
 
 // ************************************************************************* //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phasePressureModel/phasePressureModel.H

@@ -25,6 +25,23 @@ Class
     Foam::phasePressureModel
 
 Description
+    Particle-particle phase-pressure RAS model
+
+    The derivative of the phase-pressure with respect to the phase-fraction
+    is evaluated as
+
+        g0*min(exp(preAlphaExp*(alpha - alphaMax)), expMax)
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        phasePressureCoeffs
+        {
+            preAlphaExp     500;
+            expMax          1000;
+            alphaMax        0.62;
+            g0              1000;
+        }
+    \endverbatim
 
 SourceFiles
     phasePressureModel.C
@@ -34,22 +51,17 @@ SourceFiles
 #ifndef phasePressureModel_H
 #define phasePressureModel_H
 
+#include "RASModel.H"
 #include "eddyViscosity.H"
 #include "PhaseIncompressibleTurbulenceModel.H"
-#include "dragModel.H"
 #include "phaseModel.H"
-#include "autoPtr.H"
-#include "viscosityModel.H"
-#include "conductivityModel.H"
-#include "radialModel.H"
-#include "granularPressureModel.H"
-#include "frictionalStressModel.H"
-
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+namespace RASModels
+{
 
 /*---------------------------------------------------------------------------*\
                            Class phasePressureModel Declaration
@@ -57,9 +69,9 @@ namespace Foam
 
 class phasePressureModel
 :
-    public RASModels::eddyViscosity
+    public eddyViscosity
     <
-        PhaseIncompressibleTurbulenceModel<phaseModel>
+        RASModel<PhaseIncompressibleTurbulenceModel<phaseModel> >
     >
 {
     // Private data
@@ -86,6 +98,9 @@ class phasePressureModel
 
     // Private Member Functions
 
+        void correctNut()
+        {}
+
         //- Disallow default bitwise copy construct
         phasePressureModel(const phasePressureModel&);
 
@@ -93,13 +108,6 @@ class phasePressureModel
         void operator=(const phasePressureModel&);
 
 
-protected:
-
-    // Protected member functions
-
-        virtual void correctNut();
-
-
 public:
 
     //- Runtime type information
@@ -128,6 +136,9 @@ public:
 
     // Member Functions
 
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
         //- Return the effective viscosity
         virtual tmp<volScalarField> nuEff() const
         {
@@ -165,14 +176,12 @@ public:
 
         //- Solve the kinetic theory equations and correct the viscosity
         virtual void correct();
-
-        //- Re-read model coefficients if they have changed
-        virtual bool read();
 };
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+} // End namespace RASModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H

@@ -2,6 +2,7 @@
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
+      - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
       + turbulence->divDevRhoReff(rho, U)
     );
 

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H

@@ -58,7 +58,6 @@
         {
             tphiAlphaCorr() -= tphiAlpha();
 
-
             volScalarField alpha100("alpha100", alpha10);
             alpha10 = alpha1;
 

applications/test/PisoFoam/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModel/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \

applications/test/RhoPimpleFoam/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModel/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/cfdTools \

applications/test/patchRegion/Make/files

@@ -0,0 +1,3 @@
+Test-patchRegion.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-patchRegion

applications/test/patchRegion/Make/options

@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools

applications/test/patchRegion/README

@@ -0,0 +1,2 @@
+2013-05-27 Detect 'pinches' in patch (i.e., non-manifold on point). Run on
+cavity_pinched subdirectory.

applications/test/patchRegion/Test-patchRegion.C

@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Detect point pinches
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "PatchTools.H"
+#include "Time.H"
+#include "polyMesh.H"
+#include "patchEdgeFaceRegions.H"
+#include "PatchEdgeFaceWave.H"
+#include "globalIndex.H"
+#include "syncTools.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+   argList::validArgs.append("patch");
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+    const word patchName = args[1];
+
+#   include "createPolyMesh.H"
+
+    Info<< "Mesh read in = "
+        << runTime.cpuTimeIncrement()
+        << " s\n" << endl << endl;
+
+
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+    label patchI = pbm.findPatchID(patchName);
+    const polyPatch& patch = pbm[patchI];
+
+    Info<< "Patch:" << patch.name() << endl;
+
+
+
+    // Data on all edges and faces
+    List<patchEdgeFaceRegions> allEdgeInfo(patch.nEdges());
+    List<patchEdgeFaceRegions> allFaceInfo(patch.size());
+
+    // Determine parallel global indexing
+    const globalIndex globalNumbering(patch.size());
+
+    DynamicList<label> changedEdges(4*patch.size());
+    DynamicList<patchEdgeFaceRegions> changedInfo(changedEdges.size());
+
+    forAll(patch, faceI)
+    {
+        const labelList& fEdges = patch.faceEdges()[faceI];
+
+        label globalFaceI = globalNumbering.toGlobal(faceI);
+
+        forAll(fEdges, i)
+        {
+            changedEdges.append(fEdges[i]);
+            changedInfo.append
+            (
+                patchEdgeFaceRegions(labelPair(globalFaceI, globalFaceI))
+            );
+        }
+    }
+
+    // Walk
+    PatchEdgeFaceWave
+    <
+        primitivePatch,
+        patchEdgeFaceRegions
+    > calc
+    (
+        mesh,
+        patch,
+        changedEdges,
+        changedInfo,
+        allEdgeInfo,
+        allFaceInfo,
+        returnReduce(patch.nEdges(), sumOp<label>())
+    );
+
+
+    Info<< "Time now = " << runTime.timeName() << endl;
+
+
+    // Detect points with multiple regions
+    labelList duplicateRegion(patch.nPoints(), -1);
+    {
+        labelList currentRegion(patch.nPoints(), -1);
+
+        forAll(patch.localFaces(), faceI)
+        {
+            const face& f = patch.localFaces()[faceI];
+
+            forAll(f, fp)
+            {
+                label faceRegion = allFaceInfo[faceI].regions()[fp];
+
+                label pointI = f[fp];
+
+                if (currentRegion[pointI] == -1)
+                {
+                    currentRegion[pointI] = faceRegion;
+                }
+                else if (currentRegion[pointI] != faceRegion)
+                {
+                    if (duplicateRegion[pointI] == -1)
+                    {
+                        Pout<< "Multi region point:"
+                            << patch.localPoints()[pointI]
+                            << " with region:" << currentRegion[pointI]
+                            << " and region:" << faceRegion
+                            << endl;
+                        duplicateRegion[pointI] = currentRegion[pointI];
+                    }
+                }
+            }
+        }
+    }
+
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/patchRegion/cavity_pinched/README.txt

@@ -0,0 +1,17 @@
+constant/
+    cavity blockMesh
+    point (0.05 0.05 0.01) moved to (0.05 0.05 0.001)
+
+0.005/
+    collapseEdges with
+
+        collapseEdgesCoeffs.minimumEdgeLength   0.0011;
+
+    so it collapses the one edge
+
+
+processor*/0.005/
+    decomposePar
+
+
+mpirun -np 2 Test-patchRegion -parallel

applications/test/patchRegion/cavity_pinched/constant/polyMesh/blockMeshDict

@@ -0,0 +1,75 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.1;
+
+vertices
+(
+    (0 0 0)
+    (1 0 0)
+    (1 1 0)
+    (0 1 0)
+    (0 0 0.1)
+    (1 0 0.1)
+    (1 1 0.1)
+    (0 1 0.1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (20 20 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    movingWall
+    {
+        type wall;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    fixedWalls
+    {
+        type wall;
+        faces
+        (
+            (0 4 7 3)
+            (2 6 5 1)
+            (1 5 4 0)
+        );
+    }
+    frontAndBack
+    {
+        type empty;
+        faces
+        (
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

applications/test/patchRegion/cavity_pinched/constant/transportProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+nu              nu [ 0 2 -1 0 0 0 0 ] 0.01;
+
+
+// ************************************************************************* //

applications/test/patchRegion/cavity_pinched/system/collapseDict

@@ -0,0 +1,89 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      collapseDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// If on, after collapsing check the quality of the mesh. If bad faces are
+// generated then redo the collapsing with stricter filtering.
+controlMeshQuality      off;
+
+
+collapseEdgesCoeffs
+{
+    // Edges shorter than this absolute value will be merged
+    minimumEdgeLength   0.0011;
+
+    // The maximum angle between two edges that share a point attached to
+    // no other edges
+    maximumMergeAngle   30;
+}
+
+
+//collapseFacesCoeffs
+//{
+//    // The initial face length factor
+//    initialFaceLengthFactor     0.5;
+//
+//    // If the face can't be collapsed to an edge, and it has a span less than
+//    // the target face length multiplied by this coefficient, collapse it
+//    // to a point.
+//    maxCollapseFaceToPointSideLengthCoeff   0.3;
+//
+//    // Allow early collapse of edges to a point
+//    allowEarlyCollapseToPoint               on;
+//
+//    // Fraction to premultiply maxCollapseFaceToPointSideLengthCoeff by if
+//    // allowEarlyCollapseToPoint is enabled
+//    allowEarlyCollapseCoeff                 0.2;
+//
+//    // Defining how close to the midpoint (M) of the projected
+//    // vertices line a projected vertex (X) can be before making this
+//    // an invalid edge collapse
+//    //
+//    // X---X-g----------------M----X-----------g----X--X
+//    //
+//    // Only allow a collapse if all projected vertices are outwith
+//    // guardFraction (g) of the distance form the face centre to the
+//    // furthest vertex in the considered direction
+//    guardFraction                           0.1;
+//}
+//
+//
+//controlMeshQualityCoeffs
+//{
+//    // Name of the dictionary that has the mesh quality coefficients used
+//    // by motionSmoother::checkMesh
+//    #include                    "meshQualityDict";
+//
+//    // The amount that minimumEdgeLength will be reduced by for each
+//    // edge if that edge's collapse generates a poor quality face
+//    edgeReductionFactor         0.5;
+//
+//    // The amount that initialFaceLengthFactor will be reduced by for each
+//    // face if its collapse generates a poor quality face
+//    faceReductionFactor         $initialFaceLengthFactor;
+//
+//    // Maximum number of smoothing iterations for the reductionFactors
+//    maximumSmoothingIterations  2;
+//
+//    // Maximum number of outer iterations is mesh quality checking is enabled
+//    maximumIterations           10;
+//
+//    // Maximum number of iterations deletion of a point can cause a bad face
+//    // to be constructed before it is forced to not be deleted
+//    maxPointErrorCount          5;
+//}
+
+
+// ************************************************************************* //

applications/test/patchRegion/cavity_pinched/system/controlDict

@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+DebugSwitches
+{
+    PatchEdgeFaceWave 1;
+}
+
+application     icoFoam;
+
+startFrom       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0.5;
+
+deltaT          0.005;
+
+writeControl    timeStep;
+
+writeInterval   20;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //

applications/test/patchRegion/cavity_pinched/system/decomposeParDict

@@ -0,0 +1,138 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  2;
+
+//- Keep owner and neighbour on same processor for faces in zones:
+// preserveFaceZones (heater solid1 solid3);
+
+//- Keep owner and neighbour on same processor for faces in patches:
+//  (makes sense only for cyclic patches)
+//preservePatches (cyclic_half0 cyclic_half1);
+
+//- Keep all of faceSet on a single processor. This puts all cells
+//  connected with a point, edge or face on the same processor.
+//  (just having face connected cells might not guarantee a balanced
+//  decomposition)
+// The processor can be -1 (the decompositionMethod chooses the processor
+// for a good load balance) or explicitly provided (upsets balance).
+//singleProcessorFaceSets ((f0 -1));
+
+
+//- Use the volScalarField named here as a weight for each cell in the
+//  decomposition.  For example, use a particle population field to decompose
+//  for a balanced number of particles in a lagrangian simulation.
+// weightField dsmcRhoNMean;
+
+//method          scotch;
+method          hierarchical;
+// method          simple;
+// method          metis;
+// method          manual;
+// method          multiLevel;
+// method          structured;  // does 2D decomposition of structured mesh
+
+multiLevelCoeffs
+{
+    // Decomposition methods to apply in turn. This is like hierarchical but
+    // fully general - every method can be used at every level.
+
+    level0
+    {
+        numberOfSubdomains  64;
+        //method simple;
+        //simpleCoeffs
+        //{
+        //    n           (2 1 1);
+        //    delta       0.001;
+        //}
+        method scotch;
+    }
+    level1
+    {
+        numberOfSubdomains  4;
+        method scotch;
+    }
+}
+
+// Desired output
+
+simpleCoeffs
+{
+    n           (2 1 1);
+    delta       0.001;
+}
+
+hierarchicalCoeffs
+{
+    n           (2 1 1);
+    delta       0.001;
+    order       xyz;
+}
+
+metisCoeffs
+{
+ /*
+    processorWeights
+    (
+        1
+        1
+        1
+        1
+    );
+  */
+}
+
+scotchCoeffs
+{
+    //processorWeights
+    //(
+    //    1
+    //    1
+    //    1
+    //    1
+    //);
+    //writeGraph  true;
+    //strategy "b";
+}
+
+manualCoeffs
+{
+    dataFile    "decompositionData";
+}
+
+structuredCoeffs
+{
+    // Patches to do 2D decomposition on. Structured mesh only; cells have
+    // to be in 'columns' on top of patches.
+    patches     (movingWall);
+
+    // Method to use on the 2D subset
+    method      scotch;
+}
+
+//// Is the case distributed? Note: command-line argument -roots takes
+//// precedence
+//distributed     yes;
+//// Per slave (so nProcs-1 entries) the directory above the case.
+//roots
+//(
+//    "/tmp"
+//    "/tmp"
+//);
+
+// ************************************************************************* //

applications/test/patchRegion/cavity_pinched/system/fvSchemes

@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         none;
+    laplacian(nu,U) Gauss linear orthogonal;
+    laplacian((1|A(U)),p) Gauss linear orthogonal;
+}
+
+interpolationSchemes
+{
+    default         linear;
+    interpolate(HbyA) linear;
+}
+
+snGradSchemes
+{
+    default         orthogonal;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //

applications/test/patchRegion/cavity_pinched/system/fvSolution

@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        //solver          PCG;
+        //preconditioner  DIC;
+
+        solver          GAMG;
+        smoother        GaussSeidel;
+        nPreSweeps      0;
+        nPostSweeps     2;
+        cacheAgglomeration on;
+        agglomerator    faceAreaPair;
+        nCellsInCoarsestLevel 3;
+        mergeLevels     1;
+
+        processorAgglomerator   none;   //procFaces;
+        //nAgglomeratingCells     20;
+
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0;
+    }
+}
+
+PISO
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       0;
+}
+
+
+// ************************************************************************* //

applications/utilities/parallelProcessing/decomposePar/domainDecompositionDistribute.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,328 +39,33 @@ void Foam::domainDecomposition::distributeCells()
 
     cpuTime decompositionTime;
 
-
-    // See if any faces need to have owner and neighbour on same processor
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelHashSet sameProcFaces;
-
-    if (decompositionDict_.found("preservePatches"))
-    {
-        wordList pNames(decompositionDict_.lookup("preservePatches"));
-
-        Info<< nl
-            << "Keeping owner of faces in patches " << pNames
-            << " on same processor. This only makes sense for cyclics." << endl;
-
-        const polyBoundaryMesh& patches = boundaryMesh();
-
-        forAll(pNames, i)
-        {
-            const label patchI = patches.findPatchID(pNames[i]);
-
-            if (patchI == -1)
-            {
-                FatalErrorIn("domainDecomposition::distributeCells()")
-                    << "Unknown preservePatch " << pNames[i]
-                    << endl << "Valid patches are " << patches.names()
-                    << exit(FatalError);
-            }
-
-            const polyPatch& pp = patches[patchI];
-
-            forAll(pp, i)
-            {
-                sameProcFaces.insert(pp.start() + i);
-            }
-        }
-    }
-    if (decompositionDict_.found("preserveFaceZones"))
-    {
-        wordList zNames(decompositionDict_.lookup("preserveFaceZones"));
-
-        Info<< nl
-            << "Keeping owner and neighbour of faces in zones " << zNames
-            << " on same processor" << endl;
-
-        const faceZoneMesh& fZones = faceZones();
-
-        forAll(zNames, i)
-        {
-            label zoneI = fZones.findZoneID(zNames[i]);
-
-            if (zoneI == -1)
-            {
-                FatalErrorIn("domainDecomposition::distributeCells()")
-                    << "Unknown preserveFaceZone " << zNames[i]
-                    << endl << "Valid faceZones are " << fZones.names()
-                    << exit(FatalError);
-            }
-
-            const faceZone& fz = fZones[zoneI];
-
-            forAll(fz, i)
-            {
-                sameProcFaces.insert(fz[i]);
-            }
-        }
-    }
-
-
-    // Specified processor for owner and neighbour of faces
-    Map<label> specifiedProcessorFaces;
-    List<Tuple2<word, label> > zNameAndProcs;
-
-    if (decompositionDict_.found("singleProcessorFaceSets"))
-    {
-        decompositionDict_.lookup("singleProcessorFaceSets") >> zNameAndProcs;
-
-        label nCells = 0;
-
-        Info<< endl;
-
-        forAll(zNameAndProcs, i)
-        {
-            Info<< "Keeping all cells connected to faceSet "
-                << zNameAndProcs[i].first()
-                << " on processor " << zNameAndProcs[i].second() << endl;
-
-            // Read faceSet
-            faceSet fz(*this, zNameAndProcs[i].first());
-            nCells += fz.size();
-        }
-
-
-        // Size
-        specifiedProcessorFaces.resize(2*nCells);
-
-
-        // Fill
-        forAll(zNameAndProcs, i)
-        {
-            faceSet fz(*this, zNameAndProcs[i].first());
-
-            label procI = zNameAndProcs[i].second();
-
-            forAllConstIter(faceSet, fz, iter)
-            {
-                label faceI = iter.key();
-
-                specifiedProcessorFaces.insert(faceI, procI);
-            }
-        }
-    }
-
-
-    // Construct decomposition method and either do decomposition on
-    // cell centres or on agglomeration
-
-
     autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
     (
         decompositionDict_
     );
 
-
+    scalarField cellWeights;
-    if (sameProcFaces.empty() && specifiedProcessorFaces.empty())
+    if (decompositionDict_.found("weightField"))
     {
-        if (decompositionDict_.found("weightField"))
+        word weightName = decompositionDict_.lookup("weightField");
-        {
-            word weightName = decompositionDict_.lookup("weightField");
 
-            volScalarField weights
+        volScalarField weights
-            (
-                IOobject
-                (
-                    weightName,
-                    time().timeName(),
-                    *this,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                ),
-                *this
-            );
-
-            cellToProc_ = decomposePtr().decompose
-            (
-                *this,
-                cellCentres(),
-                weights.internalField()
-            );
-        }
-        else
-        {
-            cellToProc_ = decomposePtr().decompose(*this, cellCentres());
-        }
-
-    }
-    else
-    {
-        Info<< "Constrained decomposition:" << endl
-            << "    faces with same processor owner and neighbour : "
-            << sameProcFaces.size() << endl
-            << "    faces all on same processor                   : "
-            << specifiedProcessorFaces.size() << endl << endl;
-
-        // Faces where owner and neighbour are not 'connected' (= all except
-        // sameProcFaces)
-        boolList blockedFace(nFaces(), true);
-
-        forAllConstIter(labelHashSet, sameProcFaces, iter)
-        {
-            blockedFace[iter.key()] = false;
-        }
-
-
-        // For specifiedProcessorFaces add all point connected faces
-        {
-            forAllConstIter(Map<label>, specifiedProcessorFaces, iter)
-            {
-                const face& f = faces()[iter.key()];
-                forAll(f, fp)
-                {
-                    const labelList& pFaces = pointFaces()[f[fp]];
-                    forAll(pFaces, i)
-                    {
-                        blockedFace[pFaces[i]] = false;
-                    }
-                }
-            }
-        }
-
-
-        // Connect coupled boundary faces
-        const polyBoundaryMesh& patches =  boundaryMesh();
-
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            if (pp.coupled())
-            {
-                forAll(pp, i)
-                {
-                    blockedFace[pp.start()+i] = false;
-                }
-            }
-        }
-
-        // Determine global regions, separated by blockedFaces
-        regionSplit globalRegion(*this, blockedFace);
-
-
-        // Determine region cell centres
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        // This just takes the first cell in the region. Otherwise the problem
-        // is with cyclics - if we'd average the region centre might be
-        // somewhere in the middle of the domain which might not be anywhere
-        // near any of the cells.
-
-        pointField regionCentres(globalRegion.nRegions(), point::max);
-
-        forAll(globalRegion, cellI)
-        {
-            label regionI = globalRegion[cellI];
-
-            if (regionCentres[regionI] == point::max)
-            {
-                regionCentres[regionI] = cellCentres()[cellI];
-            }
-        }
-
-        // Do decomposition on agglomeration
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        scalarField regionWeights(globalRegion.nRegions(), 0);
-
-        if (decompositionDict_.found("weightField"))
-        {
-            word weightName = decompositionDict_.lookup("weightField");
-
-            volScalarField weights
-            (
-                IOobject
-                (
-                    weightName,
-                    time().timeName(),
-                    *this,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE
-                ),
-                *this
-            );
-
-            forAll(globalRegion, cellI)
-            {
-                label regionI = globalRegion[cellI];
-
-                regionWeights[regionI] += weights.internalField()[cellI];
-            }
-        }
-        else
-        {
-            forAll(globalRegion, cellI)
-            {
-                label regionI = globalRegion[cellI];
-
-                regionWeights[regionI] += 1.0;
-            }
-        }
-
-        cellToProc_ = decomposePtr().decompose
         (
-            *this,
+            IOobject
-            globalRegion,
+            (
-            regionCentres,
+                weightName,
-            regionWeights
+                time().timeName(),
+                *this,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            ),
+            *this
         );
-
+        cellWeights = weights.internalField();
-
-        // For specifiedProcessorFaces rework the cellToProc to enforce
-        // all on one processor since we can't guarantee that the input
-        // to regionSplit was a single region.
-        // E.g. faceSet 'a' with the cells split into two regions
-        // by a notch formed by two walls
-        //
-        //          \   /
-        //           \ /
-        //    ---a----+-----a-----
-        //
-        //
-        // Note that reworking the cellToProc might make the decomposition
-        // unbalanced.
-        if (specifiedProcessorFaces.size())
-        {
-            forAll(zNameAndProcs, i)
-            {
-                faceSet fz(*this, zNameAndProcs[i].first());
-
-                if (fz.size())
-                {
-                    label procI = zNameAndProcs[i].second();
-                    if (procI == -1)
-                    {
-                        // If no processor specified use the one from the
-                        // 0th element
-                        procI = cellToProc_[faceOwner()[fz[0]]];
-                    }
-
-                    forAllConstIter(faceSet, fz, iter)
-                    {
-                        label faceI = iter.key();
-
-                        cellToProc_[faceOwner()[faceI]] = procI;
-                        if (isInternalFace(faceI))
-                        {
-                            cellToProc_[faceNeighbour()[faceI]] = procI;
-                        }
-                    }
-                }
-            }
-        }
     }
 
+    cellToProc_ = decomposePtr().decompose(*this, cellWeights);
+
     Info<< "\nFinished decomposition in "
         << decompositionTime.elapsedCpuTime()
         << " s" << endl;

applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C

@@ -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
@@ -48,6 +48,7 @@ Description
 #include "faceCoupleInfo.H"
 #include "fvMeshAdder.H"
 #include "polyTopoChange.H"
+#include "zeroGradientFvPatchFields.H"
 
 using namespace Foam;
 
@@ -297,6 +298,12 @@ int main(int argc, char *argv[])
         "fullMatch",
         "do (slower) geometric matching on all boundary faces"
     );
+    argList::addBoolOption
+    (
+        "cellDist",
+        "write cell distribution as a labelList - for use with 'manual' "
+        "decomposition method or as a volScalarField for post-processing."
+    );
 
     #include "addTimeOptions.H"
     #include "addRegionOption.H"
@@ -362,6 +369,7 @@ int main(int argc, char *argv[])
             << nl << "This assumes a correct decomposition." << endl;
     }
 
+    bool writeCellDist = args.optionFound("cellDist");
 
 
     int nProcs = 0;
@@ -507,7 +515,7 @@ int main(int argc, char *argv[])
     {
         // Construct empty mesh.
         Info<< "Constructing empty mesh to add to." << nl << endl;
-        polyMesh masterMesh
+        fvMesh masterMesh
         (
             IOobject
             (
@@ -528,7 +536,7 @@ int main(int argc, char *argv[])
                 << " for time = " << databases[procI].timeName()
                 << nl << endl;
 
-            polyMesh meshToAdd
+            fvMesh meshToAdd
             (
                 IOobject
                 (
@@ -560,7 +568,7 @@ int main(int argc, char *argv[])
             // Add elements to mesh
             Info<< "Adding to master mesh" << nl << endl;
 
-            autoPtr<mapAddedPolyMesh> map = polyMeshAdder::add
+            autoPtr<mapAddedPolyMesh> map = fvMeshAdder::add
             (
                 masterMesh,
                 meshToAdd,
@@ -608,6 +616,67 @@ int main(int argc, char *argv[])
                 << "Failed writing polyMesh."
                 << exit(FatalError);
         }
+
+
+
+        if (writeCellDist)
+        {
+            // Write the decomposition as labelList for use with 'manual'
+            // decomposition method.
+            labelIOList cellDecomposition
+            (
+                IOobject
+                (
+                    "cellDecomposition",
+                    masterMesh.facesInstance(),
+                    masterMesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                masterMesh.nCells()
+            );
+
+            forAll(cellProcAddressing, procI)
+            {
+                const labelList& pCells = cellProcAddressing[procI];
+                UIndirectList<label>(cellDecomposition, pCells) = procI;
+            }
+
+            cellDecomposition.write();
+
+            Info<< nl << "Wrote decomposition to "
+                << cellDecomposition.objectPath()
+                << " for use in manual decomposition." << endl;
+
+
+            // Write as volScalarField for postprocessing.
+            volScalarField cellDist
+            (
+                IOobject
+                (
+                    "cellDist",
+                    runTime.timeName(),
+                    masterMesh,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                masterMesh,
+                dimensionedScalar("cellDist", dimless, 0),
+                zeroGradientFvPatchScalarField::typeName
+            );
+
+            forAll(cellDecomposition, cellI)
+            {
+                cellDist[cellI] = cellDecomposition[cellI];
+            }
+
+            cellDist.write();
+
+            Info<< nl << "Wrote decomposition as volScalarField to "
+                << cellDist.name() << " for use in postprocessing."
+                << endl;
+        }
     }
 
 
@@ -770,6 +839,7 @@ int main(int argc, char *argv[])
         Info<< endl;
     }
 
+
     Info<< "End.\n" << endl;
 
     return 0;

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C

@@ -974,28 +974,6 @@ int main(int argc, char *argv[])
         const word extractionMethod = surfaceDict.lookup("extractionMethod");
 
 
-#ifndef ENABLE_CURVATURE
-        if (curvature)
-        {
-            WarningIn(args.executable())
-                << "Curvature calculation has been requested but "
-                << args.executable() << " has not " << nl
-                << "    been compiled with CGAL. "
-                << "Skipping the curvature calculation." << endl;
-        }
-#else
-        if (curvature && env("FOAM_SIGFPE"))
-        {
-            WarningIn(args.executable())
-                << "Detected floating point exception trapping (FOAM_SIGFPE)."
-                << " This might give" << nl
-                << "    problems when calculating curvature on straight angles"
-                << " (infinite curvature)" << nl
-                << "    Switch it off in case of problems." << endl;
-        }
-#endif
-
-
         Info<< nl << "Feature line extraction is only valid on closed manifold "
             << "surfaces." << endl;
 

bin/foamCleanPath

@@ -3,7 +3,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
@@ -126,7 +126,7 @@ do
     if [ -e "$dir" ]
     then
         #- no duplicate dirs
-        duplicate=$(echo " $dirList " | sed -ne "s@ $dir @DUP@p")
+        duplicate=$(echo " $dirList " | sed -ne "s: $dir :DUP:p")
 
         if [ ! "$duplicate" ]
         then

src/OpenFOAM/containers/Lists/CompactListList/CompactListListI.H

@@ -223,11 +223,7 @@ inline Foam::UList<T> Foam::CompactListList<T, Container>::operator[]
 )
 {
     label start = offsets_[i];
-    return UList<T>
+    return UList<T>(m_.begin() + start, offsets_[i+1] - start);
-    (
-        (m_.size() ? m_.begin() + start : NULL),
-        offsets_[i+1] - start
-    );
 }
 
 
@@ -241,7 +237,7 @@ Foam::CompactListList<T, Container>::operator[]
     label start = offsets_[i];
     return UList<T>
     (
-        (m_.size() ? const_cast<T*>(m_.begin() + start) : NULL),
+        const_cast<T*>(m_.begin() + start),
         offsets_[i+1] - start
     );
 }

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H

@@ -34,10 +34,11 @@ SourceFiles
     PatchToolsCheck.C
     PatchToolsEdgeOwner.C
     PatchToolsGatherAndMerge.C
+    PatchToolsMatch.C
+    PatchToolsNormals.C
     PatchToolsSearch.C
     PatchToolsSortEdges.C
-    PatchToolsNormals.C
+    PatchToolsSortPoints.C
-    PatchToolsMatch.C
 
 \*---------------------------------------------------------------------------*/
 

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSortPoints.C

@@ -49,7 +49,7 @@ Foam::PatchTools::sortedPointEdges
     const labelListList& faceEdges = p.faceEdges();
 
     // create the lists for the various results. (resized on completion)
-    labelListList sortedPointEdges(pointEdges.size());
+    labelListList sortedPointEdges(pointEdges);
 
     DynamicList<label> newEdgeList;
 
@@ -57,12 +57,16 @@ Foam::PatchTools::sortedPointEdges
     {
         const labelList& pEdges = pointEdges[pointI];
 
+        label nPointEdges = pEdges.size();
+
         label edgeI = pEdges[0];
 
         label prevFaceI = edgeFaces[edgeI][0];
 
         newEdgeList.clear();
-        newEdgeList.setCapacity(pEdges.size());
+        newEdgeList.setCapacity(nPointEdges);
+
+        label nVisitedEdges = 0;
 
         do
         {
@@ -103,10 +107,42 @@ Foam::PatchTools::sortedPointEdges
 
             prevFaceI = faceI;
 
+            nVisitedEdges++;
+            if (nVisitedEdges > nPointEdges)
+            {
+                WarningIn("Foam::PatchTools::sortedPointEdges()")
+                    << "Unable to order pointEdges as the face connections "
+                    << "are not circular" << nl
+                    << "    Original pointEdges = " << pEdges << nl
+                    << "    New pointEdges = " << newEdgeList
+                    << endl;
+
+                newEdgeList = pEdges;
+
+                break;
+            }
+
         } while (edgeI != pEdges[0]);
 
-        if (newEdgeList.size() == pEdges.size())
+        if (newEdgeList.size() == nPointEdges)
         {
+            forAll(pEdges, eI)
+            {
+                if (findIndex(newEdgeList, pEdges[eI]) == -1)
+                {
+                    WarningIn("Foam::PatchTools::sortedPointEdges()")
+                        << "Cannot find all original edges in the new list"
+                        << nl
+                        << "    Original pointEdges = " << pEdges << nl
+                        << "    New pointEdges = " << newEdgeList
+                        << endl;
+
+                    newEdgeList = pEdges;
+
+                    break;
+                }
+            }
+
             sortedPointEdges[pointI] = newEdgeList;
         }
     }

src/TurbulenceModels/Allwmake

@@ -3,7 +3,7 @@ cd ${0%/*} || exit 1    # run from this directory
 makeType=${1:-libso}
 set -x
 
-wmake libso turbulenceModel
+wmake libso turbulenceModels
 wmake libso incompressible
 wmake libso compressible
 wmakeLnInclude phaseIncompressible

src/TurbulenceModels/compressible/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I../turbulenceModel/lnInclude \
+    -I../turbulenceModels/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
@@ -9,4 +9,5 @@ LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -lfluidThermophysicalModels \
+    -lturbulenceModels \
     -lspecie

src/TurbulenceModels/compressible/compressibleTurbulenceModel.H

@@ -96,7 +96,7 @@ public:
     {}
 
 
-    // Member Functions
+    // Member functions
 
         //- Return the effective stress tensor including the laminar stress
         virtual tmp<volSymmTensorField> devRhoReff() const = 0;

src/TurbulenceModels/compressible/compressibleTurbulenceModels.C

@@ -24,74 +24,39 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "CompressibleTurbulenceModel.H"
-#include "laminar.H"
-#include "RASModel.H"
-#include "kEpsilon.H"
 #include "fluidThermo.H"
 #include "addToRunTimeSelectionTable.H"
+#include "makeTurbulenceModel.H"
 
-namespace Foam
+#include "laminar.H"
-{
+#include "RASModel.H"
-    typedef TurbulenceModel
+#include "LESModel.H"
-    <
-        geometricOneField,
-        volScalarField,
-        compressibleTurbulenceModel,
-        fluidThermo
-    > baseCompressibleFluidThermoTurbulenceModel;
 
-    defineTemplateRunTimeSelectionTable
+makeBaseTurbulenceModel
-    (
+(
-        baseCompressibleFluidThermoTurbulenceModel,
+    geometricOneField,
-        dictionary
+    volScalarField,
-    );
+    compressibleTurbulenceModel,
+    CompressibleTurbulenceModel,
+    fluidThermo
+);
 
+#define makeRASModel(Type)                                                     \
+    makeTemplatedTurbulenceModel                                               \
+    (fluidThermoCompressibleTurbulenceModel, RAS, Type)
 
-    typedef CompressibleTurbulenceModel<fluidThermo>
+#define makeLESModel(Type)                                                     \
-        compressibleFluidThermoTurbulenceModel;
+    makeTemplatedTurbulenceModel                                               \
+    (fluidThermoCompressibleTurbulenceModel, LES, Type)
 
-    typedef laminar<compressibleFluidThermoTurbulenceModel> compressibleLaminar;
+#include "kEpsilon.H"
+makeRASModel(kEpsilon);
 
-    defineNamedTemplateTypeNameAndDebug(compressibleLaminar, 0);
+#include "Smagorinsky.H"
+makeLESModel(Smagorinsky);
 
-    addToRunTimeSelectionTable
+#include "kEqn.H"
-    (
+makeLESModel(kEqn);
-        baseCompressibleFluidThermoTurbulenceModel,
-        compressibleLaminar,
-        dictionary
-    );
-
-
-    typedef RASModel<compressibleFluidThermoTurbulenceModel>
-        compressibleRASModel;
-
-    defineNamedTemplateTypeNameAndDebug(compressibleRASModel, 0);
-
-    defineTemplateRunTimeSelectionTable(compressibleRASModel, dictionary);
-
-    addToRunTimeSelectionTable
-    (
-        baseCompressibleFluidThermoTurbulenceModel,
-        compressibleRASModel,
-        dictionary
-    );
-
-
-    namespace RASModels
-    {
-        typedef kEpsilon<compressibleFluidThermoTurbulenceModel>
-            compressibleKEpsilon;
-
-        defineNamedTemplateTypeNameAndDebug(compressibleKEpsilon, 0);
-
-        addToRunTimeSelectionTable
-        (
-            compressibleRASModel,
-            compressibleKEpsilon,
-            dictionary
-        );
-    }
-}
 
 
 // ************************************************************************* //

src/TurbulenceModels/incompressible/Make/options

@@ -1,10 +1,11 @@
 EXE_INC = \
-    -I../turbulenceModel/lnInclude \
+    -I../turbulenceModels/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude
 
 LIB_LIBS = \
     -lincompressibleTransportModels \
+    -lturbulenceModels \
     -lfiniteVolume \
     -lmeshTools

src/TurbulenceModels/incompressible/incompressibleTurbulenceModels.C

@@ -24,75 +24,39 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "IncompressibleTurbulenceModel.H"
-#include "laminar.H"
-#include "RASModel.H"
-#include "kEpsilon.H"
 #include "transportModel.H"
 #include "addToRunTimeSelectionTable.H"
+#include "makeTurbulenceModel.H"
 
-namespace Foam
+#include "laminar.H"
-{
+#include "RASModel.H"
-    typedef TurbulenceModel
+#include "LESModel.H"
-    <
-        geometricOneField,
-        geometricOneField,
-        incompressibleTurbulenceModel,
-        transportModel
-    > baseIncompressibleTransportTurbulenceModel;
 
-    defineTemplateRunTimeSelectionTable
+makeBaseTurbulenceModel
-    (
+(
-        baseIncompressibleTransportTurbulenceModel,
+    geometricOneField,
-        dictionary
+    geometricOneField,
-    );
+    incompressibleTurbulenceModel,
+    IncompressibleTurbulenceModel,
+    transportModel
+);
 
-    typedef IncompressibleTurbulenceModel
+#define makeRASModel(Type)                                                     \
-    <
+    makeTemplatedTurbulenceModel                                               \
-        transportModel
+    (transportModelIncompressibleTurbulenceModel, RAS, Type)
-    > incompressibleTransportTurbulenceModel;
 
-    typedef laminar<incompressibleTransportTurbulenceModel>
+#define makeLESModel(Type)                                                     \
-        incompressibleLaminar;
+    makeTemplatedTurbulenceModel                                               \
+    (transportModelIncompressibleTurbulenceModel, LES, Type)
 
-    defineNamedTemplateTypeNameAndDebug(incompressibleLaminar, 0);
+#include "kEpsilon.H"
+makeRASModel(kEpsilon);
 
-    addToRunTimeSelectionTable
+#include "Smagorinsky.H"
-    (
+makeLESModel(Smagorinsky);
-        baseIncompressibleTransportTurbulenceModel,
-        incompressibleLaminar,
-        dictionary
-    );
 
-
+#include "kEqn.H"
-    typedef RASModel<incompressibleTransportTurbulenceModel>
+makeLESModel(kEqn);
-        incompressibleRASModel;
-
-    defineNamedTemplateTypeNameAndDebug(incompressibleRASModel, 0);
-
-    defineTemplateRunTimeSelectionTable(incompressibleRASModel, dictionary);
-
-    addToRunTimeSelectionTable
-    (
-        baseIncompressibleTransportTurbulenceModel,
-        incompressibleRASModel,
-        dictionary
-    );
-
-    namespace RASModels
-    {
-        typedef kEpsilon<incompressibleTransportTurbulenceModel>
-            incompressibleKEpsilon;
-
-        defineNamedTemplateTypeNameAndDebug(incompressibleKEpsilon, 0);
-
-        addToRunTimeSelectionTable
-        (
-            incompressibleRASModel,
-            incompressibleKEpsilon,
-            dictionary
-        );
-    }
-}
 
 
 // ************************************************************************* //

src/TurbulenceModels/phaseIncompressible/LES/Niceno/NicenoKEqn.C

@@ -0,0 +1,237 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "NicenoKEqn.H"
+#include "addToRunTimeSelectionTable.H"
+#include "twoPhaseSystem.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+NicenoKEqn<BasicTurbulenceModel>::NicenoKEqn
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName,
+    const word& type
+)
+:
+    kEqn<BasicTurbulenceModel>
+    (
+        alpha,
+        rho,
+        U,
+        alphaPhi,
+        phi,
+        transport,
+        propertiesName,
+        type
+    ),
+
+    gasTurbulencePtr_(NULL),
+
+    alphaInversion_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "alphaInversion",
+            this->coeffDict_,
+            0.3
+        )
+    ),
+
+    Cp_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cp",
+            this->coeffDict_,
+            this->Ck_.value()
+        )
+    ),
+
+    Cmub_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cmub",
+            this->coeffDict_,
+            0.6
+        )
+    )
+{
+    if (type == typeName)
+    {
+        correctNut();
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool NicenoKEqn<BasicTurbulenceModel>::read()
+{
+    if (kEqn<BasicTurbulenceModel>::read())
+    {
+        alphaInversion_.readIfPresent(this->coeffDict());
+        Cp_.readIfPresent(this->coeffDict());
+        Cmub_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+const PhaseIncompressibleTurbulenceModel
+<
+    typename BasicTurbulenceModel::transportModel
+>&
+NicenoKEqn<BasicTurbulenceModel>::gasTurbulence() const
+{
+    if (!gasTurbulencePtr_)
+    {
+        const volVectorField& U = this->U_;
+
+        const transportModel& liquid = this->transport();
+        const twoPhaseSystem& fluid = liquid.fluid();
+        const transportModel& gas = fluid.otherPhase(liquid);
+
+        gasTurbulencePtr_ =
+           &U.db()
+           .lookupObject<PhaseIncompressibleTurbulenceModel<transportModel> >
+            (
+                IOobject::groupName
+                (
+                    turbulenceModel::propertiesName,
+                    gas.name()
+                )
+            );
+    }
+
+    return *gasTurbulencePtr_;
+}
+
+
+template<class BasicTurbulenceModel>
+void NicenoKEqn<BasicTurbulenceModel>::correctNut()
+{
+    const PhaseIncompressibleTurbulenceModel<transportModel>& gasTurbulence =
+        this->gasTurbulence();
+
+    this->nut_ =
+        this->Ck_*sqrt(this->k_)*this->delta()
+      + Cmub_*gasTurbulence.transport().d()*gasTurbulence.alpha()
+       *(mag(this->U_ - gasTurbulence.U()));
+
+    this->nut_.correctBoundaryConditions();
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> NicenoKEqn<BasicTurbulenceModel>::bubbleG() const
+{
+    const PhaseIncompressibleTurbulenceModel<transportModel>& gasTurbulence =
+        this->gasTurbulence();
+
+    const transportModel& liquid = this->transport();
+    const twoPhaseSystem& fluid = liquid.fluid();
+    const transportModel& gas = fluid.otherPhase(liquid);
+
+    volScalarField magUr(mag(this->U_ - gasTurbulence.U()));
+
+    tmp<volScalarField> bubbleG
+    (
+        Cp_*gas*sqr(magUr)*fluid.drag(gas).K(magUr)/liquid.rho()
+    );
+
+    return bubbleG;
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField>
+NicenoKEqn<BasicTurbulenceModel>::phaseTransferCoeff() const
+{
+    const volVectorField& U = this->U_;
+    const alphaField& alpha = this->alpha_;
+    const rhoField& rho = this->rho_;
+
+    const turbulenceModel& gasTurbulence = this->gasTurbulence();
+
+    return
+    (
+        max(alphaInversion_ - alpha, 0.0)
+       *rho
+       *min
+        (
+            this->Ce_*sqrt(gasTurbulence.k())/this->delta(),
+            1.0/U.time().deltaT()
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<fvScalarMatrix> NicenoKEqn<BasicTurbulenceModel>::kSource() const
+{
+    const alphaField& alpha = this->alpha_;
+    const rhoField& rho = this->rho_;
+
+    const PhaseIncompressibleTurbulenceModel<transportModel>& gasTurbulence =
+        this->gasTurbulence();
+
+    const volScalarField phaseTransferCoeff(this->phaseTransferCoeff());
+
+    return
+        alpha*rho*bubbleG()
+      + phaseTransferCoeff*gasTurbulence.k()
+      - fvm::Sp(phaseTransferCoeff, this->k_);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/phaseIncompressible/LES/Niceno/NicenoKEqn.H

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESModels::NicenoKEqn
+
+Group
+    grpLESTurbulence
+
+Description
+    One-equation SGS model for the continuous phase in a two-phase system
+    including bubble-generated turbulence.
+
+    Reference:
+    \verbatim
+        "One-equation sub-grid scale (SGS) modelling for Euler-Euler
+        large eddy simulation (EELES) of dispersed bubbly flow"
+        B. Niceno,
+        M.T. Dhotre,
+        N.G. Dee
+        Chemical Engineering Science 63 (2008) pp. 3923-3931.
+    \endverbatim
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        NicenoKEqnCoeffs
+        {
+            Ck              0.094;
+            Ce              1.048;
+            alphaInversion  0.3;
+            Cp              Ck;
+            Cmub            0.6;
+        }
+    \endverbatim
+
+SourceFiles
+    NicenoKEqn.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef NicenoKEqn_H
+#define NicenoKEqn_H
+
+#include "kEqn.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class NicenoKEqn Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class NicenoKEqn
+:
+    public kEqn<BasicTurbulenceModel>
+{
+    // Private data
+
+        mutable const PhaseIncompressibleTurbulenceModel
+        <
+            typename BasicTurbulenceModel::transportModel
+        > *gasTurbulencePtr_;
+
+
+    // Private Member Functions
+
+        //- Return the turbulence model for the gas phase
+        const PhaseIncompressibleTurbulenceModel
+        <
+            typename BasicTurbulenceModel::transportModel
+        >&
+        gasTurbulence() const;
+
+        // Disallow default bitwise copy construct and assignment
+        NicenoKEqn(const NicenoKEqn&);
+        NicenoKEqn& operator=(const NicenoKEqn&);
+
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar alphaInversion_;
+            dimensionedScalar Cp_;
+            dimensionedScalar Cmub_;
+
+
+    // Protected Member Functions
+
+        virtual void correctNut();
+        tmp<volScalarField> bubbleG() const;
+        tmp<volScalarField> phaseTransferCoeff() const;
+        virtual tmp<fvScalarMatrix> kSource() const;
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("NicenoKEqn");
+
+
+    // Constructors
+
+        //- Construct from components
+        NicenoKEqn
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~NicenoKEqn()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "NicenoKEqn.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/phaseIncompressible/LES/SmagorinskyZhang/SmagorinskyZhang.C

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SmagorinskyZhang.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+SmagorinskyZhang<BasicTurbulenceModel>::SmagorinskyZhang
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName,
+    const word& type
+)
+:
+    Smagorinsky<BasicTurbulenceModel>
+    (
+        alpha,
+        rho,
+        U,
+        alphaPhi,
+        phi,
+        transport,
+        propertiesName,
+        type
+    ),
+
+    gasTurbulencePtr_(NULL),
+
+    Cmub_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cmub",
+            this->coeffDict_,
+            0.6
+        )
+    )
+{
+    if (type == typeName)
+    {
+        correctNut();
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool SmagorinskyZhang<BasicTurbulenceModel>::read()
+{
+    if (Smagorinsky<BasicTurbulenceModel>::read())
+    {
+        Cmub_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+const PhaseIncompressibleTurbulenceModel
+<
+    typename BasicTurbulenceModel::transportModel
+>&
+SmagorinskyZhang<BasicTurbulenceModel>::gasTurbulence() const
+{
+    if (!gasTurbulencePtr_)
+    {
+        const volVectorField& U = this->U_;
+
+        const transportModel& liquid = this->transport();
+        const twoPhaseSystem& fluid = liquid.fluid();
+        const transportModel& gas = fluid.otherPhase(liquid);
+
+        gasTurbulencePtr_ =
+           &U.db()
+           .lookupObject<PhaseIncompressibleTurbulenceModel<transportModel> >
+            (
+                IOobject::groupName
+                (
+                    turbulenceModel::propertiesName,
+                    gas.name()
+                )
+            );
+    }
+
+    return *gasTurbulencePtr_;
+}
+
+
+template<class BasicTurbulenceModel>
+void SmagorinskyZhang<BasicTurbulenceModel>::correctNut()
+{
+    const PhaseIncompressibleTurbulenceModel<transportModel>& gasTurbulence =
+        this->gasTurbulence();
+
+    volScalarField k(this->k(fvc::grad(this->U_)));
+
+    this->nut_ =
+        this->Ck_*sqrt(k)*this->delta()
+      + Cmub_*gasTurbulence.transport().d()*gasTurbulence.alpha()
+       *(mag(this->U_ - gasTurbulence.U()));
+
+    this->nut_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/phaseIncompressible/LES/SmagorinskyZhang/SmagorinskyZhang.H

@@ -0,0 +1,170 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESModels::SmagorinskyZhang
+
+Group
+    grpLESTurbulence
+
+Description
+    The Smagorinsky SGS model including bubble-generated turbulence
+
+    Reference:
+    \verbatim
+        "Numerical simulation of the dynamic flow behavior in a bubble column:
+        A study of closures for turbulence and interface forces"
+        D. Zhang,
+        N.G. Deen,
+        J.A.M. Kuipers,
+        Chemical Engineering Science 61 (2006) pp 7593-7608.
+    \endverbatim
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        SmagorinskyZhangCoeffs
+        {
+            Ck              0.094;
+            Ce              1.048;
+            Cmub            0.6;
+        }
+    \endverbatim
+
+SourceFiles
+    SmagorinskyZhang.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SmagorinskyZhang_H
+#define SmagorinskyZhang_H
+
+#include "LESModel.H"
+#include "eddyViscosity.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class SmagorinskyZhang Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class SmagorinskyZhang
+:
+    public Smagorinsky<BasicTurbulenceModel>
+{
+    // Private data
+
+        mutable const PhaseIncompressibleTurbulenceModel
+        <
+            typename BasicTurbulenceModel::transportModel
+        > *gasTurbulencePtr_;
+
+
+    // Private Member Functions
+
+        //- Return the turbulence model for the gas phase
+        const PhaseIncompressibleTurbulenceModel
+        <
+            typename BasicTurbulenceModel::transportModel
+        >&
+        gasTurbulence() const;
+
+        // Disallow default bitwise copy construct and assignment
+        SmagorinskyZhang(const SmagorinskyZhang&);
+        SmagorinskyZhang& operator=(const SmagorinskyZhang&);
+
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar Cmub_;
+
+
+    // Protected Member Functions
+
+        virtual void correctNut();
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("SmagorinskyZhang");
+
+
+    // Constructors
+
+        //- Construct from components
+        SmagorinskyZhang
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~SmagorinskyZhang()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "SmagorinskyZhang.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/phaseIncompressible/LES/continuousGasKEqn/continuousGasKEqn.C

@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "continuousGasKEqn.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+continuousGasKEqn<BasicTurbulenceModel>::continuousGasKEqn
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName,
+    const word& type
+)
+:
+    kEqn<BasicTurbulenceModel>
+    (
+        alpha,
+        rho,
+        U,
+        alphaPhi,
+        phi,
+        transport,
+        propertiesName,
+        type
+    ),
+
+    liquidTurbulencePtr_(NULL),
+
+    alphaInversion_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "alphaInversion",
+            this->coeffDict_,
+            0.7
+        )
+    )
+{
+    if (type == typeName)
+    {
+        kEqn<BasicTurbulenceModel>::correctNut();
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool continuousGasKEqn<BasicTurbulenceModel>::read()
+{
+    if (kEqn<BasicTurbulenceModel>::read())
+    {
+        alphaInversion_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+const turbulenceModel&
+continuousGasKEqn<BasicTurbulenceModel>::liquidTurbulence() const
+{
+    if (!liquidTurbulencePtr_)
+    {
+        const volVectorField& U = this->U_;
+
+        const transportModel& gas = this->transport();
+        const twoPhaseSystem& fluid = gas.fluid();
+        const transportModel& liquid = fluid.otherPhase(gas);
+
+        liquidTurbulencePtr_ =
+           &U.db().lookupObject<turbulenceModel>
+            (
+                IOobject::groupName
+                (
+                    turbulenceModel::propertiesName,
+                    liquid.name()
+                )
+            );
+    }
+
+    return *liquidTurbulencePtr_;
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField>
+continuousGasKEqn<BasicTurbulenceModel>::phaseTransferCoeff() const
+{
+    const volVectorField& U = this->U_;
+    const alphaField& alpha = this->alpha_;
+    const rhoField& rho = this->rho_;
+
+    const turbulenceModel& liquidTurbulence = this->liquidTurbulence();
+
+    return
+    (
+        max(alphaInversion_ - alpha, 0.0)
+       *rho
+       *min
+        (
+            this->Ce_*sqrt(liquidTurbulence.k())/this->delta(),
+            1.0/U.time().deltaT()
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<fvScalarMatrix>
+continuousGasKEqn<BasicTurbulenceModel>::kSource() const
+{
+    const turbulenceModel& liquidTurbulence = this->liquidTurbulence();
+    const volScalarField phaseTransferCoeff(this->phaseTransferCoeff());
+
+    return
+        phaseTransferCoeff*liquidTurbulence.k()
+      - fvm::Sp(phaseTransferCoeff, this->k_);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/phaseIncompressible/LES/continuousGasKEqn/continuousGasKEqn.H

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESModels::continuousGasKEqn
+
+Group
+    grpLESTurbulence
+
+Description
+    One-equation SGS model for the gas-phase in a two-phase system
+    supporting phase-inversion.
+
+    In the limit that the gas-phase fraction approaches zero a contribution from
+    the other phase is blended into the k-equation up to the phase-fraction of
+    alphaInversion at which point phase-inversion is considered to have occurred
+    and the model reverts to the pure single-phase form.
+
+    This model is unpublished and is provided as a stable numerical framework
+    on which a more physical model may be built.
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        continuousKEqnCoeffs
+        {
+            Ck              0.094;
+            Ce              1.048;
+            alphaInversion  0.7;
+        }
+    \endverbatim
+
+SourceFiles
+    continuousGasKEqn.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef continuousGasKEqn_H
+#define continuousGasKEqn_H
+
+#include "kEqn.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class continuousGasKEqn Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class continuousGasKEqn
+:
+    public kEqn<BasicTurbulenceModel>
+{
+    // Private data
+
+        mutable const turbulenceModel *liquidTurbulencePtr_;
+
+
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        continuousGasKEqn(const continuousGasKEqn&);
+        continuousGasKEqn& operator=(const continuousGasKEqn&);
+
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar alphaInversion_;
+
+
+    // Protected Member Functions
+
+        //- Return the turbulence model for the liquid phase
+        const turbulenceModel& liquidTurbulence() const;
+
+        tmp<volScalarField> phaseTransferCoeff() const;
+        virtual tmp<fvScalarMatrix> kSource() const;
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("continuousGasKEqn");
+
+
+    // Constructors
+
+        //- Construct from components
+        continuousGasKEqn
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~continuousGasKEqn()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "continuousGasKEqn.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/phaseIncompressible/RAS/LaheyKEpsilon/LaheyKEpsilon.H

@@ -28,6 +28,31 @@ Group
     grpRASTurbulence
 
 Description
+    Continuous-phase k-epsilon model including bubble-generated turbulence.
+
+    Reference:
+    \verbatim
+        "The simulation of multidimensional multiphase flows",
+        Lahey R.T.,
+        Nucl. Eng. & Design
+        2005 (235) pp.1043-1060.
+    \endverbatim
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        LaheyKEpsilonCoeffs
+        {
+            Cmu             0.09;
+            C1              1.44;
+            C2              1.92;
+            C3              -0.33;
+            sigmak          1.0;
+            sigmaEps        1.3;
+            Cp              0.25;
+            Cmub            0.6;
+            alphaInversion  0.3;
+        }
+    \endverbatim
 
 SourceFiles
     LaheyKEpsilon.C
@@ -63,6 +88,20 @@ class LaheyKEpsilon
         > *gasTurbulencePtr_;
 
 
+    // Private Member Functions
+
+        //- Return the turbulence model for the gas phase
+        const PhaseIncompressibleTurbulenceModel
+        <
+            typename BasicTurbulenceModel::transportModel
+        >&
+        gasTurbulence() const;
+
+        // Disallow default bitwise copy construct and assignment
+        LaheyKEpsilon(const LaheyKEpsilon&);
+        LaheyKEpsilon& operator=(const LaheyKEpsilon&);
+
+
 protected:
 
     // Protected data
@@ -74,7 +113,7 @@ protected:
             dimensionedScalar Cmub_;
 
 
-    // Protected member functions
+    // Protected Member Functions
 
         virtual void correctNut();
         tmp<volScalarField> bubbleG() const;
@@ -82,6 +121,7 @@ protected:
         virtual tmp<fvScalarMatrix> kSource() const;
         virtual tmp<fvScalarMatrix> epsilonSource() const;
 
+
 public:
 
     typedef typename BasicTurbulenceModel::alphaField alphaField;
@@ -116,15 +156,11 @@ public:
 
     // Member Functions
 
-        //- Return the turbulence model for the gas phase
+        //- Read model coefficients if they have changed
-        const PhaseIncompressibleTurbulenceModel<transportModel>&
+        virtual bool read();
-        gasTurbulence() const;
 
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
-
-        //- Read RASProperties dictionary
-        virtual bool read();
 };
 
 

src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.H

@@ -28,6 +28,31 @@ Group
     grpRASTurbulence
 
 Description
+    k-epsilon model for the gas-phase in a two-phase system
+    supporting phase-inversion.
+
+    In the limit that the gas-phase fraction approaches zero a contribution from
+    the other phase is blended into the k and epsilon equations up to the
+    phase-fraction of alphaInversion at which point phase-inversion is
+    considered to have occurred and the model reverts to the pure single-phase
+    form.
+
+    This model is unpublished and is provided as a stable numerical framework
+    on which a more physical model may be built.
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        continuousGasKEpsilonCoeffs
+        {
+            Cmu             0.09;
+            C1              1.44;
+            C2              1.92;
+            C3              -0.33;
+            sigmak          1.0;
+            sigmaEps        1.3;
+            alphaInversion  0.7;
+        }
+    \endverbatim
 
 SourceFiles
     continuousGasKEpsilon.C
@@ -62,6 +87,13 @@ class continuousGasKEpsilon
         volScalarField nutEff_;
 
 
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        continuousGasKEpsilon(const continuousGasKEpsilon&);
+        continuousGasKEpsilon& operator=(const continuousGasKEpsilon&);
+
+
 protected:
 
     // Protected data
@@ -71,7 +103,7 @@ protected:
             dimensionedScalar alphaInversion_;
 
 
-    // Protected member functions
+    // Protected Member Functions
 
         virtual void correctNut();
         tmp<volScalarField> phaseTransferCoeff() const;
@@ -113,6 +145,9 @@ public:
 
     // Member Functions
 
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
         //- Return the turbulence model for the liquid phase
         const turbulenceModel& liquidTurbulence() const;
 
@@ -124,9 +159,6 @@ public:
 
         //- Return the Reynolds stress tensor
         virtual tmp<volSymmTensorField> R() const;
-
-        //- Read RASProperties dictionary
-        virtual bool read();
 };
 
 

src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.C

@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "LESModel.H"
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+void Foam::LESModel<BasicTurbulenceModel>::printCoeffs(const word& type)
+{
+    if (printCoeffs_)
+    {
+        Info<< type << "Coeffs" << coeffDict_ << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+Foam::LESModel<BasicTurbulenceModel>::LESModel
+(
+    const word& type,
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName
+)
+:
+    BasicTurbulenceModel
+    (
+        alpha,
+        rho,
+        U,
+        alphaPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    LESDict_(this->subOrEmptyDict("LES")),
+    turbulence_(LESDict_.lookup("turbulence")),
+    printCoeffs_(LESDict_.lookupOrDefault<Switch>("printCoeffs", false)),
+    coeffDict_(LESDict_.subOrEmptyDict(type + "Coeffs")),
+
+    kMin_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "kMin",
+            LESDict_,
+            SMALL,
+            sqr(dimVelocity)
+        )
+    ),
+
+    delta_(LESdelta::New("delta", U.mesh(), LESDict_))
+{
+    // Force the construction of the mesh deltaCoeffs which may be needed
+    // for the construction of the derived models and BCs
+    this->mesh_.deltaCoeffs();
+}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+Foam::autoPtr<Foam::LESModel<BasicTurbulenceModel> >
+Foam::LESModel<BasicTurbulenceModel>::New
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName
+)
+{
+    // get model name, but do not register the dictionary
+    // otherwise it is registered in the database twice
+    const word modelType
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                IOobject::groupName(propertiesName, U.group()),
+                U.time().constant(),
+                U.db(),
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).subDict("LES").lookup("LESModel")
+    );
+
+    Info<< "Selecting LES turbulence model " << modelType << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "LESModel::New"
+            "("
+                "const volScalarField&, "
+                "const volVectorField&, "
+                "const surfaceScalarField&, "
+                "transportModel&, "
+                "const word&"
+            ")"
+        )   << "Unknown LESModel type "
+            << modelType << nl << nl
+            << "Valid LESModel types:" << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<LESModel>
+    (
+        cstrIter()(alpha, rho, U, alphaPhi, phi, transport, propertiesName)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+void Foam::LESModel<BasicTurbulenceModel>::correct()
+{
+    BasicTurbulenceModel::correct();
+}
+
+
+template<class BasicTurbulenceModel>
+bool Foam::LESModel<BasicTurbulenceModel>::read()
+{
+    if (turbulenceModel::read())
+    {
+        LESDict_ <<= this->subDict("LES");
+        LESDict_.lookup("turbulence") >> turbulence_;
+
+        if (const dictionary* dictPtr = LESDict_.subDictPtr(type() + "Coeffs"))
+        {
+            coeffDict_ <<= *dictPtr;
+        }
+
+        kMin_.readIfPresent(LESDict_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESModel/LESModel.H

@@ -0,0 +1,240 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Namespace
+    Foam::LESModels
+
+Description
+    Namespace for LES SGS models.
+
+Class
+    Foam::LESModel
+
+Description
+    Templated abstract base class for LES SGS models
+
+SourceFiles
+    LESModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LESModel_H
+#define LESModel_H
+
+#include "TurbulenceModel.H"
+#include "LESdelta.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class LESModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class LESModel
+:
+    public BasicTurbulenceModel
+{
+
+protected:
+
+    // Protected data
+
+        //- LES coefficients dictionary
+        dictionary LESDict_;
+
+        //- Turbulence on/off flag
+        Switch turbulence_;
+
+        //- Flag to print the model coeffs at run-time
+        Switch printCoeffs_;
+
+        //- Model coefficients dictionary
+        dictionary coeffDict_;
+
+        //- Lower limit of k
+        dimensionedScalar kMin_;
+
+        //- Run-time selectable delta model
+        autoPtr<Foam::LESdelta> delta_;
+
+
+    // Protected Member Functions
+
+        //- Print model coefficients
+        virtual void printCoeffs(const word& type);
+
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        LESModel(const LESModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const LESModel&);
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("LES");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            LESModel,
+            dictionary,
+            (
+                const alphaField& alpha,
+                const rhoField& rho,
+                const volVectorField& U,
+                const surfaceScalarField& alphaPhi,
+                const surfaceScalarField& phi,
+                const transportModel& transport,
+                const word& propertiesName
+            ),
+            (alpha, rho, U, alphaPhi, phi, transport, propertiesName)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        LESModel
+        (
+            const word& type,
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected LES model
+        static autoPtr<LESModel> New
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName
+        );
+
+
+    //- Destructor
+    virtual ~LESModel()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+
+        // Access
+
+            //- Const access to the coefficients dictionary
+            virtual const dictionary& coeffDict() const
+            {
+                return coeffDict_;
+            }
+
+            //- Return the lower allowable limit for k (default: SMALL)
+            const dimensionedScalar& kMin() const
+            {
+                return kMin_;
+            }
+
+            //- Allow kMin to be changed
+            dimensionedScalar& kMin()
+            {
+                return kMin_;
+            }
+
+            //- Access function to filter width
+            inline const volScalarField& delta() const
+            {
+                return delta_();
+            }
+
+
+        //- Return the effective viscosity
+        virtual tmp<volScalarField> nuEff() const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    IOobject::groupName("nuEff", this->U_.group()),
+                    this->nut() + this->nu()
+                )
+            );
+        }
+
+        //- Return the effective viscosity on patch
+        virtual tmp<scalarField> nuEff(const label patchi) const
+        {
+            return this->nut(patchi) + this->nu(patchi);
+        }
+
+        //- Solve the turbulence equations and correct the turbulence viscosity
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "LESModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/LESdelta/LESdelta.C

@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "LESdelta.H"
+#include "calculatedFvPatchFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(LESdelta, 0);
+    defineRunTimeSelectionTable(LESdelta, dictionary);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::LESdelta::LESdelta(const word& name, const fvMesh& mesh)
+:
+    mesh_(mesh),
+    delta_
+    (
+        IOobject
+        (
+            name,
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar(name, dimLength, SMALL),
+        calculatedFvPatchScalarField::typeName
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::LESdelta> Foam::LESdelta::New
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+{
+    const word deltaType(dict.lookup("delta"));
+
+    Info<< "Selecting LES delta type " << deltaType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(deltaType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "LESdelta::New(const fvMesh&, const dictionary&)"
+        )   << "Unknown LESdelta type "
+            << deltaType << nl << nl
+            << "Valid LESdelta types are :" << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<LESdelta>(cstrIter()(name, mesh, dict));
+}
+
+
+Foam::autoPtr<Foam::LESdelta> Foam::LESdelta::New
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict,
+    const dictionaryConstructorTable& additionalConstructors
+)
+{
+    const word deltaType(dict.lookup("delta"));
+
+    Info<< "Selecting LES delta type " << deltaType << endl;
+
+    // First on additional ones
+    dictionaryConstructorTable::const_iterator cstrIter =
+        additionalConstructors.find(deltaType);
+
+    if (cstrIter != additionalConstructors.end())
+    {
+        return autoPtr<LESdelta>(cstrIter()(name, mesh, dict));
+    }
+    else
+    {
+        dictionaryConstructorTable::const_iterator cstrIter =
+            dictionaryConstructorTablePtr_->find(deltaType);
+
+        if (cstrIter == dictionaryConstructorTablePtr_->end())
+        {
+            FatalErrorIn
+            (
+                "LESdelta::New(const fvMesh&, const dictionary&)"
+            )   << "Unknown LESdelta type "
+                << deltaType << nl << nl
+                << "Valid LESdelta types are :" << endl
+                << additionalConstructors.sortedToc()
+                << " and "
+                << dictionaryConstructorTablePtr_->sortedToc()
+                << exit(FatalError);
+            return autoPtr<LESdelta>();
+        }
+        else
+        {
+            return autoPtr<LESdelta>(cstrIter()(name, mesh, dict));
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/LESdelta/LESdelta.H

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESdelta
+
+Description
+    Abstract base class for LES deltas
+
+SourceFiles
+    LESdelta.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LESdelta_H
+#define LESdelta_H
+
+#include "volFields.H"
+#include "typeInfo.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+                          Class LESdelta Declaration
+\*---------------------------------------------------------------------------*/
+
+class LESdelta
+{
+
+protected:
+
+    // Protected data
+
+        const fvMesh& mesh_;
+
+        volScalarField delta_;
+
+
+private:
+
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        LESdelta(const LESdelta&);
+        void operator=(const LESdelta&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("LESdelta");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            LESdelta,
+            dictionary,
+            (
+                const word& name,
+                const fvMesh& mesh,
+                const dictionary& LESdeltaDict
+            ),
+            (name, mesh, LESdeltaDict)
+        );
+
+
+    // Constructors
+
+        //- Construct from name and mesh
+        LESdelta(const word& name, const fvMesh&);
+
+
+    // Selectors
+
+        //- Return a reference to the selected LES delta
+        static autoPtr<LESdelta> New
+        (
+            const word& name,
+            const fvMesh&,
+            const dictionary&
+        );
+
+        //- Return a reference to the selected LES delta
+        static autoPtr<LESdelta> New
+        (
+            const word& name,
+            const fvMesh&,
+            const dictionary&,
+            const dictionaryConstructorTable&
+        );
+
+
+    //- Destructor
+    virtual ~LESdelta()
+    {}
+
+
+    // Member Functions
+
+        //- Return mesh reference
+        const fvMesh& mesh() const
+        {
+            return mesh_;
+        }
+
+        //- Read the LESdelta dictionary
+        virtual void read(const dictionary&) = 0;
+
+        // Correct values
+        virtual void correct() = 0;
+
+
+    // Member Operators
+
+        virtual operator const volScalarField&() const
+        {
+            return delta_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/PrandtlDelta/PrandtlDelta.C

@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PrandtlDelta.H"
+#include "wallDist.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(PrandtlDelta, 0);
+    addToRunTimeSelectionTable(LESdelta, PrandtlDelta, dictionary);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::PrandtlDelta::calcDelta()
+{
+    delta_ = min
+    (
+        static_cast<const volScalarField&>(geometricDelta_()),
+        (kappa_/Cdelta_)*wallDist(mesh_).y()
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::PrandtlDelta::PrandtlDelta
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    LESdelta(name, mesh),
+    geometricDelta_(LESdelta::New(name, mesh, dict.subDict(type() + "Coeffs"))),
+    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
+    Cdelta_
+    (
+        dict.subDict(type() + "Coeffs").lookupOrDefault<scalar>("Cdelta", 0.158)
+    )
+{
+    calcDelta();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::PrandtlDelta::read(const dictionary& dict)
+{
+    const dictionary& coeffDict(dict.subDict(type() + "Coeffs"));
+
+    geometricDelta_().read(coeffDict);
+    dict.readIfPresent<scalar>("kappa", kappa_);
+    coeffDict.readIfPresent<scalar>("Cdelta", Cdelta_);
+    calcDelta();
+}
+
+
+void Foam::PrandtlDelta::correct()
+{
+    geometricDelta_().correct();
+
+    if (mesh_.changing())
+    {
+        calcDelta();
+    }
+}
+
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/PrandtlDelta/PrandtlDelta.H

@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::PrandtlDelta
+
+Description
+    Simple cube-root of cell volume delta used in LES models.
+
+SourceFiles
+    PrandtlDelta.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PrandtlDelta_H
+#define PrandtlDelta_H
+
+#include "LESdelta.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class PrandtlDelta Declaration
+\*---------------------------------------------------------------------------*/
+
+class PrandtlDelta
+:
+    public LESdelta
+{
+    // Private data
+
+        autoPtr<LESdelta> geometricDelta_;
+        scalar kappa_;
+        scalar Cdelta_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct and assignment
+        PrandtlDelta(const PrandtlDelta&);
+        void operator=(const PrandtlDelta&);
+
+        // Calculate the delta values
+        void calcDelta();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Prandtl");
+
+
+    // Constructors
+
+        //- Construct from name, mesh and IOdictionary
+        PrandtlDelta
+        (
+            const word& name,
+            const fvMesh& mesh,
+            const dictionary&
+        );
+
+
+    //- Destructor
+    virtual ~PrandtlDelta()
+    {}
+
+
+    // Member Functions
+
+        //- Read the LESdelta dictionary
+        virtual void read(const dictionary&);
+
+        // Correct values
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.C

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "cubeRootVolDelta.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(cubeRootVolDelta, 0);
+    addToRunTimeSelectionTable(LESdelta, cubeRootVolDelta, dictionary);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::cubeRootVolDelta::calcDelta()
+{
+    label nD = mesh().nGeometricD();
+
+    if (nD == 3)
+    {
+        delta_.internalField() = deltaCoeff_*pow(mesh().V(), 1.0/3.0);
+    }
+    else if (nD == 2)
+    {
+        WarningIn("cubeRootVolDelta::calcDelta()")
+            << "Case is 2D, LES is not strictly applicable\n"
+            << endl;
+
+        const Vector<label>& directions = mesh().geometricD();
+
+        scalar thickness = 0.0;
+        for (direction dir=0; dir<directions.nComponents; dir++)
+        {
+            if (directions[dir] == -1)
+            {
+                thickness = mesh().bounds().span()[dir];
+                break;
+            }
+        }
+
+        delta_.internalField() = deltaCoeff_*sqrt(mesh().V()/thickness);
+    }
+    else
+    {
+        FatalErrorIn("cubeRootVolDelta::calcDelta()")
+            << "Case is not 3D or 2D, LES is not applicable"
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cubeRootVolDelta::cubeRootVolDelta
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    LESdelta(name, mesh),
+    deltaCoeff_
+    (
+        dict.subDict(type() + "Coeffs").lookupOrDefault<scalar>("deltaCoeff", 1)
+    )
+{
+    calcDelta();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::cubeRootVolDelta::read(const dictionary& dict)
+{
+    dict.subDict(type() + "Coeffs").readIfPresent<scalar>
+    (
+        "deltaCoeff",
+        deltaCoeff_
+    );
+
+    calcDelta();
+}
+
+
+void Foam::cubeRootVolDelta::correct()
+{
+    if (mesh_.changing())
+    {
+        calcDelta();
+    }
+}
+
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/cubeRootVolDelta/cubeRootVolDelta.H

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::cubeRootVolDelta
+
+Description
+    Simple cube-root of cell volume delta used in LES models.
+
+SourceFiles
+    cubeRootVolDelta.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cubeRootVolDelta_H
+#define cubeRootVolDelta_H
+
+#include "LESdelta.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class cubeRootVolDelta Declaration
+\*---------------------------------------------------------------------------*/
+
+class cubeRootVolDelta
+:
+    public LESdelta
+{
+    // Private data
+
+        scalar deltaCoeff_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct and assignment
+        cubeRootVolDelta(const cubeRootVolDelta&);
+        void operator=(const cubeRootVolDelta&);
+
+        // Calculate the delta values
+        void calcDelta();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("cubeRootVol");
+
+
+    // Constructors
+
+        //- Construct from name, mesh and IOdictionary
+        cubeRootVolDelta
+        (
+            const word& name,
+            const fvMesh& mesh,
+            const dictionary&
+        );
+
+
+    //- Destructor
+    virtual ~cubeRootVolDelta()
+    {}
+
+
+    // Member Functions
+
+        //- Read the LESdelta dictionary
+        virtual void read(const dictionary&);
+
+        // Correct values
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C

@@ -0,0 +1,145 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "maxDeltaxyz.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(maxDeltaxyz, 0);
+    addToRunTimeSelectionTable(LESdelta, maxDeltaxyz, dictionary);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::maxDeltaxyz::calcDelta()
+{
+    label nD = mesh().nGeometricD();
+
+    tmp<volScalarField> hmax
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "hmax",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("zrero", dimLength, 0.0)
+        )
+    );
+
+    const cellList& cells = mesh().cells();
+
+    forAll(cells,cellI)
+    {
+        scalar deltaMaxTmp = 0.0;
+        const labelList& cFaces = mesh().cells()[cellI];
+        const point& centrevector = mesh().cellCentres()[cellI];
+
+        forAll(cFaces, cFaceI)
+        {
+            label faceI = cFaces[cFaceI];
+            const point& facevector = mesh().faceCentres()[faceI];
+            scalar tmp = mag(facevector - centrevector);
+            if (tmp > deltaMaxTmp)
+            {
+                deltaMaxTmp = tmp;
+            }
+        }
+        hmax()[cellI] = deltaCoeff_*deltaMaxTmp;
+    }
+
+    if (nD == 3)
+    {
+        delta_.internalField() = hmax();
+    }
+    else if (nD == 2)
+    {
+        WarningIn("maxDeltaxyz::calcDelta()")
+            << "Case is 2D, LES is not strictly applicable\n"
+            << endl;
+
+        delta_.internalField() = hmax();
+    }
+    else
+    {
+        FatalErrorIn("maxDeltaxyz::calcDelta()")
+            << "Case is not 3D or 2D, LES is not applicable"
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::maxDeltaxyz::maxDeltaxyz
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    LESdelta(name, mesh),
+    deltaCoeff_
+    (
+        dict.subDict(type() + "Coeffs").lookupOrDefault<scalar>("deltaCoeff", 1)
+    )
+{
+    calcDelta();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::maxDeltaxyz::read(const dictionary& dict)
+{
+    dict.subDict(type() + "Coeffs").readIfPresent<scalar>
+    (
+        "deltaCoeff",
+        deltaCoeff_
+    );
+
+    calcDelta();
+}
+
+
+void Foam::maxDeltaxyz::correct()
+{
+    if (mesh_.changing())
+    {
+        calcDelta();
+    }
+}
+
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.H

@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::maxDeltaxyz
+
+Description
+    maxDeltaxyz takes the maximum of the three dimensions per cell:
+    max(hx, hy, hz). Valid for structures hexahedral cells only.
+
+
+SourceFiles
+    maxDeltaxyz.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef maxDeltaxyzDelta_H
+#define maxDeltaxyzDelta_H
+
+#include "LESdelta.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class maxDeltaxyz Declaration
+\*---------------------------------------------------------------------------*/
+
+class maxDeltaxyz
+:
+    public LESdelta
+{
+    // Private data
+
+        scalar deltaCoeff_; //
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct and assignment
+        maxDeltaxyz(const maxDeltaxyz&);
+        void operator=(const maxDeltaxyz&);
+
+        // Calculate the delta values
+        void calcDelta();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("maxDeltaxyz");
+
+
+    // Constructors
+
+        //- Construct from name, mesh and IOdictionary
+        maxDeltaxyz
+        (
+            const word& name,
+            const fvMesh& mesh,
+            const dictionary&
+        );
+
+
+    //- Destructor
+    virtual ~maxDeltaxyz()
+    {}
+
+
+    // Member Functions
+
+        //- Read the LESdelta dictionary
+        virtual void read(const dictionary&);
+
+        // Correct values
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.C

@@ -0,0 +1,192 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "smoothDelta.H"
+#include "addToRunTimeSelectionTable.H"
+#include "FaceCellWave.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(smoothDelta, 0);
+    addToRunTimeSelectionTable(LESdelta, smoothDelta, dictionary);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Fill changedFaces (with face labels) and changedFacesInfo (with delta)
+// This is the initial set of faces from which to start the waves.
+// Since there might be lots of places with delta jumps we can follow various
+// strategies for this initial 'seed'.
+// - start from single cell/face and let FaceCellWave pick up all others
+//   from there. might be quite a few waves before everything settles.
+// - start from all faces. Lots of initial transfers.
+// We do something inbetween:
+// - start from all faces where there is a jump. Since we cannot easily
+//   determine this across coupled patches (cyclic, processor) introduce
+//   all faces of these and let FaceCellWave sort it out.
+void Foam::smoothDelta::setChangedFaces
+(
+    const polyMesh& mesh,
+    const volScalarField& delta,
+    DynamicList<label>& changedFaces,
+    DynamicList<deltaData>& changedFacesInfo
+)
+{
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        scalar ownDelta = delta[mesh.faceOwner()[faceI]];
+
+        scalar neiDelta = delta[mesh.faceNeighbour()[faceI]];
+
+        // Check if owner delta much larger than neighbour delta or vice versa
+
+        if (ownDelta > maxDeltaRatio_ * neiDelta)
+        {
+            changedFaces.append(faceI);
+            changedFacesInfo.append(deltaData(ownDelta));
+        }
+        else if (neiDelta > maxDeltaRatio_ * ownDelta)
+        {
+            changedFaces.append(faceI);
+            changedFacesInfo.append(deltaData(neiDelta));
+        }
+    }
+
+    // Insert all faces of coupled patches no matter what. Let FaceCellWave
+    // sort it out.
+    forAll(mesh.boundaryMesh(), patchI)
+    {
+        const polyPatch& patch = mesh.boundaryMesh()[patchI];
+
+        if (patch.coupled())
+        {
+            forAll(patch, patchFaceI)
+            {
+                label meshFaceI = patch.start() + patchFaceI;
+
+                scalar ownDelta = delta[mesh.faceOwner()[meshFaceI]];
+
+                changedFaces.append(meshFaceI);
+                changedFacesInfo.append(deltaData(ownDelta));
+            }
+        }
+    }
+
+    changedFaces.shrink();
+    changedFacesInfo.shrink();
+}
+
+
+void Foam::smoothDelta::calcDelta()
+{
+    const volScalarField& geometricDelta = geometricDelta_();
+
+    // Fill changed faces with info
+    DynamicList<label> changedFaces(mesh_.nFaces()/100 + 100);
+    DynamicList<deltaData> changedFacesInfo(changedFaces.size());
+
+    setChangedFaces(mesh_, geometricDelta, changedFaces, changedFacesInfo);
+
+    // Set initial field on cells.
+    List<deltaData> cellDeltaData(mesh_.nCells());
+
+    forAll(geometricDelta, cellI)
+    {
+        cellDeltaData[cellI] = geometricDelta[cellI];
+    }
+
+    // Set initial field on faces.
+    List<deltaData> faceDeltaData(mesh_.nFaces());
+
+
+    // Propagate information over whole domain.
+    FaceCellWave<deltaData, scalar> deltaCalc
+    (
+        mesh_,
+        changedFaces,
+        changedFacesInfo,
+        faceDeltaData,
+        cellDeltaData,
+        mesh_.globalData().nTotalCells()+1,  // max iterations
+        maxDeltaRatio_
+    );
+
+    forAll(delta_, cellI)
+    {
+        delta_[cellI] = cellDeltaData[cellI].delta();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::smoothDelta::smoothDelta
+(
+    const word& name,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    LESdelta(name, mesh),
+    geometricDelta_
+    (
+        LESdelta::New("geometricDelta", mesh, dict.subDict(type() + "Coeffs"))
+    ),
+    maxDeltaRatio_
+    (
+        readScalar(dict.subDict(type() + "Coeffs").lookup("maxDeltaRatio"))
+    )
+{
+    calcDelta();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::smoothDelta::read(const dictionary& dict)
+{
+    const dictionary& coeffsDict(dict.subDict(type() + "Coeffs"));
+
+    geometricDelta_().read(coeffsDict);
+    coeffsDict.lookup("maxDeltaRatio") >> maxDeltaRatio_;
+    calcDelta();
+}
+
+
+void Foam::smoothDelta::correct()
+{
+    geometricDelta_().correct();
+
+    if (mesh_.changing())
+    {
+        calcDelta();
+    }
+}
+
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.H

@@ -0,0 +1,290 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::smoothDelta
+
+Description
+    Smoothed delta which takes a given simple geometric delta and applies
+    smoothing to it such that the ratio of deltas between two cells is no
+    larger than a specified amount, typically 1.15.
+
+SourceFiles
+    smoothDelta.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef smoothDelta_H
+#define smoothDelta_H
+
+#include "LESdelta.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class smoothDelta Declaration
+\*---------------------------------------------------------------------------*/
+
+class smoothDelta
+:
+    public LESdelta
+{
+public:
+
+    //- Public member class used by mesh-wave to propagate the delta-ratio
+    class deltaData
+    {
+        scalar delta_;
+
+        // Private Member Functions
+
+            //- Update. Gets information from neighbouring face/cell and
+            //  uses this to update itself (if necessary) and return true.
+            template<class TrackingData>
+            inline bool update
+            (
+                const deltaData& w2,
+                const scalar scale,
+                const scalar tol,
+                TrackingData& td
+            );
+
+
+    public:
+
+        // Constructors
+
+            //- Construct null
+            inline deltaData();
+
+            //- Construct from delta value
+            inline deltaData(const scalar delta);
+
+
+        // Member Functions
+
+            // Access
+
+            scalar delta() const
+            {
+                return delta_;
+            }
+
+
+            // Needed by FaceCellWave
+
+                //- Check whether origin has been changed at all or
+                //  still contains original (invalid) value.
+                template<class TrackingData>
+                inline bool valid(TrackingData& td) const;
+
+                //- Check for identical geometrical data.
+                //  Used for cyclics checking.
+                template<class TrackingData>
+                inline bool sameGeometry
+                (
+                    const polyMesh&,
+                    const deltaData&,
+                    const scalar,
+                    TrackingData& td
+                ) const;
+
+                //- Convert any absolute coordinates into relative to
+                //  (patch)face centre
+                template<class TrackingData>
+                inline void leaveDomain
+                (
+                    const polyMesh&,
+                    const polyPatch&,
+                    const label patchFaceI,
+                    const point& faceCentre,
+                    TrackingData& td
+                );
+
+                //- Reverse of leaveDomain
+                template<class TrackingData>
+                inline void enterDomain
+                (
+                    const polyMesh&,
+                    const polyPatch&,
+                    const label patchFaceI,
+                    const point& faceCentre,
+                    TrackingData& td
+                );
+
+                //- Apply rotation matrix to any coordinates
+                template<class TrackingData>
+                inline void transform
+                (
+                    const polyMesh&,
+                    const tensor&,
+                    TrackingData& td
+                );
+
+                //- Influence of neighbouring face.
+                template<class TrackingData>
+                inline bool updateCell
+                (
+                    const polyMesh&,
+                    const label thisCellI,
+                    const label neighbourFaceI,
+                    const deltaData& neighbourInfo,
+                    const scalar tol,
+                    TrackingData& td
+                );
+
+                //- Influence of neighbouring cell.
+                template<class TrackingData>
+                inline bool updateFace
+                (
+                    const polyMesh&,
+                    const label thisFaceI,
+                    const label neighbourCellI,
+                    const deltaData& neighbourInfo,
+                    const scalar tol,
+                    TrackingData& td
+                );
+
+                //- Influence of different value on same face.
+                template<class TrackingData>
+                inline bool updateFace
+                (
+                    const polyMesh&,
+                    const label thisFaceI,
+                    const deltaData& neighbourInfo,
+                    const scalar tol,
+                    TrackingData& td
+                );
+
+                //- Same (like operator==)
+                template<class TrackingData>
+                inline bool equal(const deltaData&, TrackingData& td) const;
+
+            // Member Operators
+
+                // Needed for List IO
+                inline bool operator==(const deltaData&) const;
+
+                inline bool operator!=(const deltaData&) const;
+
+            // IOstream Operators
+
+                friend Ostream& operator<<
+                (
+                    Ostream& os,
+                    const deltaData& wDist
+                )
+                {
+                    return os << wDist.delta_;
+                }
+
+                friend Istream& operator>>(Istream& is, deltaData& wDist)
+                {
+                    return is >> wDist.delta_;
+                }
+    };
+
+
+private:
+
+    // Private data
+
+        autoPtr<LESdelta> geometricDelta_;
+        scalar maxDeltaRatio_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct and assignment
+        smoothDelta(const smoothDelta&);
+        void operator=(const smoothDelta&);
+
+        // Calculate the delta values
+        void calcDelta();
+
+
+        void setChangedFaces
+        (
+            const polyMesh& mesh,
+            const volScalarField& delta,
+            DynamicList<label>& changedFaces,
+            DynamicList<deltaData>& changedFacesInfo
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("smooth");
+
+
+    // Constructors
+
+        //- Construct from name, mesh and IOdictionary
+        smoothDelta
+        (
+            const word& name,
+            const fvMesh& mesh,
+            const dictionary&
+        );
+
+
+    //- Destructor
+    virtual ~smoothDelta()
+    {}
+
+
+    // Member Functions
+
+        //- Read the LESdelta dictionary
+        virtual void read(const dictionary&);
+
+        // Correct values
+        virtual void correct();
+};
+
+
+//- Data associated with deltaData type are contiguous
+template<>
+inline bool contiguous<smoothDelta::deltaData>()
+{
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "smoothDeltaDeltaDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDeltaDeltaDataI.H

@@ -0,0 +1,227 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Update this with w2 if applicable
+template<class TrackingData>
+inline bool smoothDelta::deltaData::update
+(
+    const smoothDelta::deltaData& w2,
+    const scalar scale,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    if (!valid(td) || (delta_ < VSMALL))
+    {
+        // My delta not set. Take over neighbour.
+        delta_ = w2.delta()/scale;
+
+        // Something changed. Let caller know.
+        return true;
+    }
+    else if (w2.delta() > (1 + tol)*scale*delta_)
+    {
+        // Neighbour is too big for me. Up my delta.
+        delta_ = w2.delta()/scale;
+
+        // Something changed. Let caller know.
+        return true;
+    }
+    else
+    {
+        // Neighbour is not too big for me or change is too small
+        // Nothing changed.
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Null constructor
+inline smoothDelta::deltaData::deltaData()
+:
+    delta_(-GREAT)
+{}
+
+
+// Construct from components
+inline smoothDelta::deltaData::deltaData(const scalar delta)
+:
+    delta_(delta)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class TrackingData>
+inline bool smoothDelta::deltaData::valid(TrackingData& td) const
+{
+    return delta_ > -SMALL;
+}
+
+
+// Checks for cyclic faces
+template<class TrackingData>
+inline bool smoothDelta::deltaData::sameGeometry
+(
+    const polyMesh&,
+    const deltaData&,
+    const scalar,
+    TrackingData& td
+) const
+{
+    return true;
+}
+
+
+template<class TrackingData>
+inline void smoothDelta::deltaData::leaveDomain
+(
+    const polyMesh&,
+    const polyPatch&,
+    const label,
+    const point&,
+    TrackingData& td
+)
+{}
+
+
+template<class TrackingData>
+inline void smoothDelta::deltaData::transform
+(
+    const polyMesh&,
+    const tensor&,
+    TrackingData& td
+)
+{}
+
+
+// Update absolute geometric quantities.
+template<class TrackingData>
+inline void smoothDelta::deltaData::enterDomain
+(
+    const polyMesh&,
+    const polyPatch&,
+    const label,
+    const point&,
+    TrackingData& td
+)
+{}
+
+
+// Update this (cellI) with face information.
+template<class TrackingData>
+inline bool smoothDelta::deltaData::updateCell
+(
+    const polyMesh&,
+    const label,
+    const label,
+    const deltaData& neighbourWallInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    // Take over info from face if more than deltaRatio larger.
+    return update(neighbourWallInfo, td, tol, td);
+}
+
+
+// Update this (face) with cell information.
+template<class TrackingData>
+inline bool smoothDelta::deltaData::updateFace
+(
+    const polyMesh&,
+    const label,
+    const label,
+    const deltaData& neighbourWallInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    // Take over information from cell without any scaling (scale = 1.0)
+    return update(neighbourWallInfo, 1.0, tol, td);
+}
+
+
+// Update this (face) with coupled face information.
+template<class TrackingData>
+inline bool smoothDelta::deltaData::updateFace
+(
+    const polyMesh&,
+    const label,
+    const deltaData& neighbourWallInfo,
+    const scalar tol,
+    TrackingData& td
+)
+{
+    // Take over information from coupled face without any scaling (scale = 1.0)
+    return update(neighbourWallInfo, 1.0, tol, td);
+}
+
+
+template<class TrackingData>
+inline bool smoothDelta::deltaData::equal
+(
+    const deltaData& rhs,
+    TrackingData& td
+) const
+{
+    return operator==(rhs);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool smoothDelta::deltaData::operator==
+(
+    const deltaData& rhs
+) const
+{
+    return delta_ == rhs.delta();
+}
+
+
+inline bool smoothDelta::deltaData::operator!=
+(
+    const deltaData& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESfilters/LESfilter/LESfilter.C

@@ -0,0 +1,67 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "LESfilter.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(LESfilter, 0);
+    defineRunTimeSelectionTable(LESfilter, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::LESfilter> Foam::LESfilter::New
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+{
+    const word filterType(dict.lookup("filter"));
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(filterType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "LESfilter::New(const fvMesh&, const dictionary&)"
+        )   << "Unknown LESfilter type "
+            << filterType << nl << nl
+            << "Valid LESfilter types are :" << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<LESfilter>(cstrIter()(mesh, dict));
+}
+
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESfilters/LESfilter/LESfilter.H

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESfilter
+
+Description
+    Abstract class for LES filters
+
+SourceFiles
+    LESfilter.C
+    newFilter.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LESfilter_H
+#define LESfilter_H
+
+#include "volFields.H"
+#include "typeInfo.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class LESfilter Declaration
+\*---------------------------------------------------------------------------*/
+
+class LESfilter
+{
+    // Private data
+
+        const fvMesh& mesh_;
+
+
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        LESfilter(const LESfilter&);
+        void operator=(const LESfilter&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("LESfilter");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            LESfilter,
+            dictionary,
+            (
+                const fvMesh& mesh,
+                const dictionary& LESfilterDict
+            ),
+            (mesh, LESfilterDict)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        LESfilter(const fvMesh& mesh)
+        :
+            mesh_(mesh)
+        {}
+
+
+    // Selectors
+
+        //- Return a reference to the selected LES filter
+        static autoPtr<LESfilter> New
+        (
+            const fvMesh&,
+            const dictionary&
+        );
+
+
+    //- Destructor
+    virtual ~LESfilter()
+    {}
+
+
+    // Member Functions
+
+        //- Return mesh reference
+        const fvMesh& mesh() const
+        {
+            return mesh_;
+        }
+
+        //- Read the LESfilter dictionary
+        virtual void read(const dictionary&) = 0;
+
+
+    // Member Operators
+
+        virtual tmp<volScalarField> operator()
+        (
+            const tmp<volScalarField>&
+        ) const = 0;
+
+        virtual tmp<volVectorField> operator()
+        (
+            const tmp<volVectorField>&
+        ) const = 0;
+
+        virtual tmp<volSymmTensorField> operator()
+        (
+            const tmp<volSymmTensorField>&
+        ) const = 0;
+
+        virtual tmp<volTensorField> operator()
+        (
+            const tmp<volTensorField>&
+        ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESfilters/anisotropicFilter/anisotropicFilter.C

@@ -0,0 +1,239 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "anisotropicFilter.H"
+#include "addToRunTimeSelectionTable.H"
+#include "zeroGradientFvPatchFields.H"
+#include "wallFvPatch.H"
+#include "fvc.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(anisotropicFilter, 0);
+    addToRunTimeSelectionTable(LESfilter, anisotropicFilter, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::anisotropicFilter::anisotropicFilter
+(
+    const fvMesh& mesh,
+    scalar widthCoeff
+)
+:
+    LESfilter(mesh),
+    widthCoeff_(widthCoeff),
+    coeff_
+    (
+        IOobject
+        (
+            "anisotropicFilterCoeff",
+            mesh.time().timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedVector("zero", dimLength*dimLength, vector::zero),
+        calculatedFvPatchVectorField::typeName
+    )
+{
+    for (direction d=0; d<vector::nComponents; d++)
+    {
+        coeff_.internalField().replace
+        (
+            d,
+            (1/widthCoeff_)*
+            sqr
+            (
+                2.0*mesh.V()
+               /fvc::surfaceSum(mag(mesh.Sf().component(d)))().internalField()
+            )
+        );
+    }
+}
+
+
+Foam::anisotropicFilter::anisotropicFilter
+(
+    const fvMesh& mesh,
+    const dictionary& bd
+)
+:
+    LESfilter(mesh),
+    widthCoeff_(readScalar(bd.subDict(type() + "Coeffs").lookup("widthCoeff"))),
+    coeff_
+    (
+        IOobject
+        (
+            "anisotropicFilterCoeff",
+            mesh.time().timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedVector("zero", dimLength*dimLength, vector::zero),
+        calculatedFvPatchScalarField::typeName
+    )
+{
+    for (direction d=0; d<vector::nComponents; d++)
+    {
+        coeff_.internalField().replace
+        (
+            d,
+            (1/widthCoeff_)*
+            sqr
+            (
+                2.0*mesh.V()
+               /fvc::surfaceSum(mag(mesh.Sf().component(d)))().internalField()
+            )
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::anisotropicFilter::read(const dictionary& bd)
+{
+    bd.subDict(type() + "Coeffs").lookup("widthCoeff") >> widthCoeff_;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::anisotropicFilter::operator()
+(
+    const tmp<volScalarField>& unFilteredField
+) const
+{
+    tmp<volScalarField> tmpFilteredField =
+        unFilteredField
+      + (
+           coeff_
+         & fvc::surfaceIntegrate
+           (
+               mesh().Sf()
+              *fvc::snGrad(unFilteredField())
+           )
+        );
+
+    unFilteredField.clear();
+
+    return tmpFilteredField;
+}
+
+
+Foam::tmp<Foam::volVectorField> Foam::anisotropicFilter::operator()
+(
+    const tmp<volVectorField>& unFilteredField
+) const
+{
+    tmp<volVectorField> tmpFilteredField =
+        unFilteredField
+      + (
+           coeff_
+         & fvc::surfaceIntegrate
+           (
+               mesh().Sf()
+              *fvc::snGrad(unFilteredField())
+           )
+        );
+
+    unFilteredField.clear();
+
+    return tmpFilteredField;
+}
+
+
+Foam::tmp<Foam::volSymmTensorField> Foam::anisotropicFilter::operator()
+(
+    const tmp<volSymmTensorField>& unFilteredField
+) const
+{
+    tmp<volSymmTensorField> tmpFilteredField
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                "anisotropicFilteredSymmTensorField",
+                mesh().time().timeName(),
+                mesh()
+            ),
+            mesh(),
+            unFilteredField().dimensions()
+        )
+    );
+
+    for (direction d=0; d<symmTensor::nComponents; d++)
+    {
+        tmpFilteredField().replace
+        (
+            d, anisotropicFilter::operator()(unFilteredField().component(d))
+        );
+    }
+
+    unFilteredField.clear();
+
+    return tmpFilteredField;
+}
+
+
+Foam::tmp<Foam::volTensorField> Foam::anisotropicFilter::operator()
+(
+    const tmp<volTensorField>& unFilteredField
+) const
+{
+    tmp<volTensorField> tmpFilteredField
+    (
+        new volTensorField
+        (
+            IOobject
+            (
+                "anisotropicFilteredTensorField",
+                mesh().time().timeName(),
+                mesh()
+            ),
+            mesh(),
+            unFilteredField().dimensions()
+        )
+    );
+
+    for (direction d=0; d<tensor::nComponents; d++)
+    {
+        tmpFilteredField().replace
+        (
+            d, anisotropicFilter::operator()(unFilteredField().component(d))
+        );
+    }
+
+    unFilteredField.clear();
+
+    return tmpFilteredField;
+}
+
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESfilters/anisotropicFilter/anisotropicFilter.H

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::anisotropicFilter
+
+Description
+    anisotropic filter
+
+    \verbatim
+    Kernel                 as filter          as Test filter with ratio 2
+    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    Box filter:            g = delta2/24  ->  g = delta2/6
+    Spherical box filter:  g = delta2/64  ->  g = delta2/16
+    Gaussian filter:       g = delta2/24  ->  g = delta2/6
+    \endverbatim
+
+SourceFiles
+    anisotropicFilter.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef anisotropicFilter_H
+#define anisotropicFilter_H
+
+#include "LESfilter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class anisotropicFilter Declaration
+\*---------------------------------------------------------------------------*/
+
+class anisotropicFilter
+:
+    public LESfilter
+{
+    // Private data
+
+        scalar widthCoeff_;
+        volVectorField coeff_;
+
+
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        anisotropicFilter(const anisotropicFilter&);
+        void operator=(const anisotropicFilter&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("anisotropic");
+
+    // Constructors
+
+        //- Construct from components
+        anisotropicFilter(const fvMesh& mesh, scalar widthCoeff);
+
+        //- Construct from IOdictionary
+        anisotropicFilter(const fvMesh& mesh, const dictionary&);
+
+
+    //- Destructor
+    virtual ~anisotropicFilter()
+    {}
+
+
+    // Member Functions
+
+        //- Read the LESfilter dictionary
+        virtual void read(const dictionary&);
+
+
+    // Member Operators
+
+        virtual tmp<volScalarField> operator()
+        (
+            const tmp<volScalarField>&
+        ) const;
+
+        virtual tmp<volVectorField> operator()
+        (
+            const tmp<volVectorField>&
+        ) const;
+
+        virtual tmp<volSymmTensorField> operator()
+        (
+            const tmp<volSymmTensorField>&
+        ) const;
+
+        virtual tmp<volTensorField> operator()
+        (
+            const tmp<volTensorField>&
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESfilters/laplaceFilter/laplaceFilter.C

@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "laplaceFilter.H"
+#include "addToRunTimeSelectionTable.H"
+#include "calculatedFvPatchFields.H"
+#include "fvm.H"
+#include "fvc.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(laplaceFilter, 0);
+    addToRunTimeSelectionTable(LESfilter, laplaceFilter, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::laplaceFilter::laplaceFilter(const fvMesh& mesh, scalar widthCoeff)
+:
+    LESfilter(mesh),
+    widthCoeff_(widthCoeff),
+    coeff_
+    (
+        IOobject
+        (
+            "laplaceFilterCoeff",
+            mesh.time().timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("zero", dimLength*dimLength, 0),
+        calculatedFvPatchScalarField::typeName
+    )
+{
+    coeff_.dimensionedInternalField() = pow(mesh.V(), 2.0/3.0)/widthCoeff_;
+}
+
+
+Foam::laplaceFilter::laplaceFilter(const fvMesh& mesh, const dictionary& bd)
+:
+    LESfilter(mesh),
+    widthCoeff_(readScalar(bd.subDict(type() + "Coeffs").lookup("widthCoeff"))),
+    coeff_
+    (
+        IOobject
+        (
+            "laplaceFilterCoeff",
+            mesh.time().timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("zero", dimLength*dimLength, 0),
+        calculatedFvPatchScalarField::typeName
+    )
+{
+    coeff_.dimensionedInternalField() = pow(mesh.V(), 2.0/3.0)/widthCoeff_;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::laplaceFilter::read(const dictionary& bd)
+{
+    bd.subDict(type() + "Coeffs").lookup("widthCoeff") >> widthCoeff_;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::laplaceFilter::operator()
+(
+    const tmp<volScalarField>& unFilteredField
+) const
+{
+    tmp<volScalarField> filteredField =
+        unFilteredField() + fvc::laplacian(coeff_, unFilteredField());
+
+    unFilteredField.clear();
+
+    return filteredField;
+}
+
+
+Foam::tmp<Foam::volVectorField> Foam::laplaceFilter::operator()
+(
+    const tmp<volVectorField>& unFilteredField
+) const
+{
+    tmp<volVectorField> filteredField =
+        unFilteredField() + fvc::laplacian(coeff_, unFilteredField());
+
+    unFilteredField.clear();
+
+    return filteredField;
+}
+
+
+Foam::tmp<Foam::volSymmTensorField> Foam::laplaceFilter::operator()
+(
+    const tmp<volSymmTensorField>& unFilteredField
+) const
+{
+    tmp<volSymmTensorField> filteredField =
+        unFilteredField() + fvc::laplacian(coeff_, unFilteredField());
+
+    unFilteredField.clear();
+
+    return filteredField;
+}
+
+
+Foam::tmp<Foam::volTensorField> Foam::laplaceFilter::operator()
+(
+    const tmp<volTensorField>& unFilteredField
+) const
+{
+    tmp<volTensorField> filteredField =
+        unFilteredField() + fvc::laplacian(coeff_, unFilteredField());
+
+    unFilteredField.clear();
+
+    return filteredField;
+}
+
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESfilters/laplaceFilter/laplaceFilter.H

@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::laplaceFilter
+
+Description
+    Laplace filter for LES
+
+    \verbatim
+    Kernel                 as filter          as Test filter with ratio 2
+    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    Box filter:            g = delta2/24  ->  g = delta2/6
+    Spherical box filter:  g = delta2/64  ->  g = delta2/16
+    Gaussian filter:       g = delta2/24  ->  g = delta2/6
+    \endverbatim
+
+SourceFiles
+    laplaceFilter.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef laplaceFilter_H
+#define laplaceFilter_H
+
+#include "LESfilter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class laplaceFilter Declaration
+\*---------------------------------------------------------------------------*/
+
+class laplaceFilter
+:
+    public LESfilter
+{
+    // Private data
+
+        scalar widthCoeff_;
+        volScalarField coeff_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct and assignment
+        laplaceFilter(const laplaceFilter&);
+        void operator=(const laplaceFilter&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("laplace");
+
+
+    // Constructors
+
+        //- Construct from components
+        laplaceFilter(const fvMesh& mesh, scalar widthCoeff);
+
+        //- Construct from IOdictionary
+        laplaceFilter(const fvMesh& mesh, const dictionary&);
+
+
+    //- Destructor
+    virtual ~laplaceFilter()
+    {}
+
+
+    // Member Functions
+
+        //- Read the LESfilter dictionary
+        virtual void read(const dictionary&);
+
+
+    // Member Operators
+
+        virtual tmp<volScalarField> operator()
+        (
+            const tmp<volScalarField>&
+        ) const;
+
+        virtual tmp<volVectorField> operator()
+        (
+            const tmp<volVectorField>&
+        ) const;
+
+        virtual tmp<volSymmTensorField> operator()
+        (
+            const tmp<volSymmTensorField>&
+        ) const;
+
+        virtual tmp<volTensorField> operator()
+        (
+            const tmp<volTensorField>&
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESfilters/simpleFilter/simpleFilter.C

@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "simpleFilter.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvc.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(simpleFilter, 0);
+    addToRunTimeSelectionTable(LESfilter, simpleFilter, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::simpleFilter::simpleFilter
+(
+    const fvMesh& mesh
+)
+:
+    LESfilter(mesh)
+{}
+
+
+Foam::simpleFilter::simpleFilter(const fvMesh& mesh, const dictionary&)
+:
+    LESfilter(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::simpleFilter::read(const dictionary&)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::simpleFilter::operator()
+(
+    const tmp<volScalarField>& unFilteredField
+) const
+{
+    tmp<volScalarField> filteredField = fvc::surfaceSum
+    (
+        mesh().magSf()*fvc::interpolate(unFilteredField)
+    )/fvc::surfaceSum(mesh().magSf());
+
+    unFilteredField.clear();
+
+    return filteredField;
+}
+
+
+Foam::tmp<Foam::volVectorField> Foam::simpleFilter::operator()
+(
+    const tmp<volVectorField>& unFilteredField
+) const
+{
+    tmp<volVectorField> filteredField = fvc::surfaceSum
+    (
+        mesh().magSf()*fvc::interpolate(unFilteredField)
+    )/fvc::surfaceSum(mesh().magSf());
+
+    unFilteredField.clear();
+
+    return filteredField;
+}
+
+
+Foam::tmp<Foam::volSymmTensorField> Foam::simpleFilter::operator()
+(
+    const tmp<volSymmTensorField>& unFilteredField
+) const
+{
+    tmp<volSymmTensorField> filteredField = fvc::surfaceSum
+    (
+        mesh().magSf()*fvc::interpolate(unFilteredField)
+    )/fvc::surfaceSum(mesh().magSf());
+
+    unFilteredField.clear();
+
+    return filteredField;
+}
+
+
+Foam::tmp<Foam::volTensorField> Foam::simpleFilter::operator()
+(
+    const tmp<volTensorField>& unFilteredField
+) const
+{
+    tmp<volTensorField> filteredField = fvc::surfaceSum
+    (
+        mesh().magSf()*fvc::interpolate(unFilteredField)
+    )/fvc::surfaceSum(mesh().magSf());
+
+    unFilteredField.clear();
+
+    return filteredField;
+}
+
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESfilters/simpleFilter/simpleFilter.H

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::simpleFilter
+
+Description
+    Simple top-hat filter used in dynamic LES models.
+
+    Implemented as a surface integral of the face interpolate of the field.
+
+SourceFiles
+    simpleFilter.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef simpleFilter_H
+#define simpleFilter_H
+
+#include "LESfilter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class simpleFilter Declaration
+\*---------------------------------------------------------------------------*/
+
+class simpleFilter
+:
+    public LESfilter
+
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct and assignment
+        simpleFilter(const simpleFilter&);
+        void operator=(const simpleFilter&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("simple");
+
+
+    // Constructors
+
+        //- Construct from components
+        simpleFilter(const fvMesh& mesh);
+
+        //- Construct from IOdictionary
+        simpleFilter(const fvMesh& mesh, const dictionary&);
+
+
+    //- Destructor
+    virtual ~simpleFilter()
+    {}
+
+
+    // Member Functions
+
+        //- Read the LESfilter dictionary
+        virtual void read(const dictionary&);
+
+
+    // Member Operators
+
+        virtual tmp<volScalarField> operator()
+        (
+            const tmp<volScalarField>&
+        ) const;
+
+        virtual tmp<volVectorField> operator()
+        (
+            const tmp<volVectorField>&
+        ) const;
+
+        virtual tmp<volSymmTensorField> operator()
+        (
+            const tmp<volSymmTensorField>&
+        ) const;
+
+        virtual tmp<volTensorField> operator()
+        (
+            const tmp<volTensorField>&
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/Smagorinsky/Smagorinsky.C

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Smagorinsky.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+Smagorinsky<BasicTurbulenceModel>::Smagorinsky
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName,
+    const word& type
+)
+:
+    LESeddyViscosity<BasicTurbulenceModel>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    Ck_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Ck",
+            this->coeffDict_,
+            0.094
+        )
+    )
+{
+    if (type == typeName)
+    {
+        correctNut();
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool Smagorinsky<BasicTurbulenceModel>::read()
+{
+    if (LESeddyViscosity<BasicTurbulenceModel>::read())
+    {
+        Ck_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> Smagorinsky<BasicTurbulenceModel>::k
+(
+    const tmp<volTensorField>& gradU
+) const
+{
+    volSymmTensorField D(symm(gradU));
+
+    volScalarField a(this->Ce_/this->delta());
+    volScalarField b((2.0/3.0)*tr(D));
+    volScalarField c(2*Ck_*this->delta()*(dev(D) && D));
+
+    return sqr((-b + sqrt(sqr(b) + 4*a*c))/(2*a));
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> Smagorinsky<BasicTurbulenceModel>::epsilon() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("epsilon", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            this->Ce_*k()*sqrt(k())/this->delta()
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+void Smagorinsky<BasicTurbulenceModel>::correctNut()
+{
+    volScalarField k(this->k(fvc::grad(this->U_)));
+
+    this->nut_ = Ck_*this->delta()*sqrt(k);
+    this->nut_.correctBoundaryConditions();
+}
+
+
+template<class BasicTurbulenceModel>
+void Smagorinsky<BasicTurbulenceModel>::correct()
+{
+    LESeddyViscosity<BasicTurbulenceModel>::correct();
+    correctNut();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/Smagorinsky/Smagorinsky.H

@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESModels::Smagorinsky
+
+Group
+    grpLESTurbulence
+
+Description
+    The Smagorinsky SGS model.
+
+    Algebraic eddy viscosity SGS model founded on the assumption that
+    local equilibrium prevails.
+    Thus,
+    \verbatim
+
+        B = 2/3*k*I - 2*nuSgs*dev(D)
+
+    where
+
+        D = symm(grad(U));
+        k from D:B + Ce*k^3/2/delta = 0
+        nuSgs = Ck*sqrt(k)*delta
+    \endverbatim
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        SmagorinskyCoeffs
+        {
+            Ck                  0.094;
+            Ce                  1.048;
+        }
+    \endverbatim
+
+SourceFiles
+    Smagorinsky.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Smagorinsky_H
+#define Smagorinsky_H
+
+#include "LESModel.H"
+#include "LESeddyViscosity.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Smagorinsky Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class Smagorinsky
+:
+    public LESeddyViscosity<BasicTurbulenceModel>
+{
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        Smagorinsky(const Smagorinsky&);
+        Smagorinsky& operator=(const Smagorinsky&);
+
+
+protected:
+
+    // Protected data
+
+        dimensionedScalar Ck_;
+
+
+    // Protected Member Functions
+
+        //- Return SGS kinetic energy
+        //  calculated from the given velocity gradient
+        tmp<volScalarField> k(const tmp<volTensorField>& gradU) const;
+
+        virtual void correctNut();
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("Smagorinsky");
+
+
+    // Constructors
+
+        //- Construct from components
+        Smagorinsky
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~Smagorinsky()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return SGS kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return k(fvc::grad(this->U_));
+        }
+
+        //- Return sub-grid disipation rate
+        virtual tmp<volScalarField> epsilon() const;
+
+        //- Correct Eddy-Viscosity and related properties
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Smagorinsky.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/eddyViscosity/LESeddyViscosity.C

@@ -0,0 +1,118 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "LESeddyViscosity.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+LESeddyViscosity<BasicTurbulenceModel>::LESeddyViscosity
+(
+    const word& type,
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName
+)
+:
+    eddyViscosity<LESModel<BasicTurbulenceModel> >
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    Ce_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Ce",
+            this->coeffDict_,
+            1.048
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool LESeddyViscosity<BasicTurbulenceModel>::read()
+{
+    if (eddyViscosity<LESModel<BasicTurbulenceModel> >::read())
+    {
+        Ce_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> LESeddyViscosity<BasicTurbulenceModel>::epsilon() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("epsilon", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            Ce_*this->k()*sqrt(this->k())/this->delta()
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/eddyViscosity/LESeddyViscosity.H

@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESModels::LESeddyViscosity
+
+Group
+    grpLESTurbulence
+
+Description
+    Eddy viscosity LES SGS model base class
+
+SourceFiles
+    LESeddyViscosity.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LESeddyViscosity_H
+#define LESeddyViscosity_H
+
+#include "LESModel.H"
+#include "eddyViscosity.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class LESeddyViscosity Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class LESeddyViscosity
+:
+    public eddyViscosity<LESModel<BasicTurbulenceModel> >
+{
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        LESeddyViscosity(const LESeddyViscosity&);
+        LESeddyViscosity& operator=(const LESeddyViscosity&);
+
+
+protected:
+
+    // Protected data
+
+        dimensionedScalar Ce_;
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    // Constructors
+
+        //- Construct from components
+        LESeddyViscosity
+        (
+            const word& type,
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName
+        );
+
+
+    //- Destructor
+    virtual ~LESeddyViscosity()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return sub-grid disipation rate
+        virtual tmp<volScalarField> epsilon() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "LESeddyViscosity.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/kEqn/kEqn.C

@@ -0,0 +1,205 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "kEqn.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+kEqn<BasicTurbulenceModel>::kEqn
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName,
+    const word& type
+)
+:
+    LESeddyViscosity<BasicTurbulenceModel>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    k_
+    (
+        IOobject
+        (
+            IOobject::groupName("k", this->U_.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+
+    Ck_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Ck",
+            this->coeffDict_,
+            0.094
+        )
+    )
+{
+    if (type == typeName)
+    {
+        correctNut();
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool kEqn<BasicTurbulenceModel>::read()
+{
+    if (LESeddyViscosity<BasicTurbulenceModel>::read())
+    {
+        Ck_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kEqn<BasicTurbulenceModel>::epsilon() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("epsilon", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            this->Ce_*k()*sqrt(k())/this->delta()
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+void kEqn<BasicTurbulenceModel>::correctNut()
+{
+    this->nut_ = Ck_*sqrt(k_)*this->delta();
+    this->nut_.correctBoundaryConditions();
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<fvScalarMatrix> kEqn<BasicTurbulenceModel>::kSource() const
+{
+    return tmp<fvScalarMatrix>
+    (
+        new fvScalarMatrix
+        (
+            k_,
+            dimVolume*this->rho_.dimensions()*k_.dimensions()
+            /dimTime
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+void kEqn<BasicTurbulenceModel>::correct()
+{
+    if (!this->turbulence_)
+    {
+        return;
+    }
+
+    // Local references
+    const alphaField& alpha = this->alpha_;
+    const rhoField& rho = this->rho_;
+    const surfaceScalarField& alphaPhi = this->alphaPhi_;
+    const surfaceScalarField& phi = this->phi_;
+    const volVectorField& U = this->U_;
+    volScalarField& nut = this->nut_;
+
+    LESeddyViscosity<BasicTurbulenceModel>::correct();
+
+    volScalarField divU(fvc::div(fvc::absolute(phi/fvc::interpolate(rho), U)));
+
+    tmp<volTensorField> tgradU(fvc::grad(U));
+    volScalarField G(this->GName(), nut*(tgradU() && dev(twoSymm(tgradU()))));
+    tgradU.clear();
+
+    tmp<fvScalarMatrix> kEqn
+    (
+        fvm::ddt(alpha, rho, k_)
+      + fvm::div(alphaPhi, k_)
+      - fvm::Sp(fvc::ddt(alpha, rho) + fvc::div(alphaPhi), k_)
+      - fvm::laplacian(alpha*rho*DkEff(), k_)
+     ==
+        alpha*rho*G
+      - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
+      - fvm::Sp(this->Ce_*alpha*rho*sqrt(k_)/this->delta(), k_)
+      + kSource()
+    );
+
+    kEqn().relax();
+    solve(kEqn);
+    bound(k_, this->kMin_);
+
+    correctNut();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/kEqn/kEqn.H

@@ -0,0 +1,184 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESModels::kEqn
+
+Group
+    grpLESTurbulence
+
+Description
+    One Equation Eddy Viscosity Model
+
+    Eddy viscosity SGS model using a modeled balance equation to simulate the
+    behaviour of k, hence,
+    \verbatim
+        d/dt(rho*k) + div(rho*U*k) - div(rho*nuEff*grad(k))
+        =
+        -rho*D:B - Ce*rho*k^(3/2)/delta
+
+    and
+
+        B = 2/3*k*I - 2*nuSgs*dev(D)
+
+    where
+
+        D = symm(grad(U));
+        nuSgs = Ck*sqrt(k)*delta
+        nuEff = nuSgs + nu
+    \endverbatim
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        NicenoKEqnCoeffs
+        {
+            Ck                  0.094;
+            Ce                  1.048;
+        }
+    \endverbatim
+
+SourceFiles
+    kEqn.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kEqn_H
+#define kEqn_H
+
+#include "LESModel.H"
+#include "LESeddyViscosity.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class kEqn Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class kEqn
+:
+    public LESeddyViscosity<BasicTurbulenceModel>
+{
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        kEqn(const kEqn&);
+        kEqn& operator=(const kEqn&);
+
+
+protected:
+
+    // Protected data
+
+        volScalarField k_;
+
+        dimensionedScalar Ck_;
+
+
+    // Protected Member Functions
+
+        virtual void correctNut();
+        virtual tmp<fvScalarMatrix> kSource() const;
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("kEqn");
+
+
+    // Constructors
+
+        //- Constructor from components
+        kEqn
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~kEqn()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return SGS kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return k_;
+        }
+
+        //- Return sub-grid disipation rate
+        virtual tmp<volScalarField> epsilon() const;
+
+        //- Return the effective diffusivity for k
+        tmp<volScalarField> DkEff() const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField("DkEff", this->nut_ + this->nu())
+            );
+        }
+
+        //- Correct Eddy-Viscosity and related properties
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "kEqn.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/Make/files

@@ -1,5 +1,21 @@
 turbulenceModel.C
 
+LESdelta = LES/LESdeltas
+
+$(LESdelta)/LESdelta/LESdelta.C
+$(LESdelta)/cubeRootVolDelta/cubeRootVolDelta.C
+$(LESdelta)/PrandtlDelta/PrandtlDelta.C
+$(LESdelta)/smoothDelta/smoothDelta.C
+$(LESdelta)/maxDeltaxyz/maxDeltaxyz.C
+
+LESfilters = LES/LESfilters
+
+$(LESfilters)/LESfilter/LESfilter.C
+$(LESfilters)/simpleFilter/simpleFilter.C
+$(LESfilters)/laplaceFilter/laplaceFilter.C
+$(LESfilters)/anisotropicFilter/anisotropicFilter.C
+
+
 derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C
 derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchField.C
 
@@ -44,6 +60,4 @@ RAS/derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLeng
 RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
 RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
 
-/* backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C */
-
 LIB = $(FOAM_LIBBIN)/libturbulenceModels

src/TurbulenceModels/turbulenceModels/RAS/RASModel/RASModel.C

@@ -68,14 +68,39 @@ Foam::RASModel<BasicTurbulenceModel>::RASModel
     printCoeffs_(RASDict_.lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(RASDict_.subOrEmptyDict(type + "Coeffs")),
 
-    kMin_("kMin", sqr(dimVelocity), SMALL),
+    kMin_
-    epsilonMin_("epsilonMin", kMin_.dimensions()/dimTime, SMALL),
+    (
-    omegaMin_("omegaMin", dimless/dimTime, SMALL)
+        dimensioned<scalar>::lookupOrAddToDict
-{
+        (
-    kMin_.readIfPresent(RASDict_);
+            "kMin",
-    epsilonMin_.readIfPresent(RASDict_);
+            RASDict_,
-    omegaMin_.readIfPresent(RASDict_);
+            SMALL,
+            sqr(dimVelocity)
+        )
+    ),
 
+    epsilonMin_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "epsilonMin",
+            RASDict_,
+            SMALL,
+            kMin_.dimensions()/dimTime
+        )
+    ),
+
+    omegaMin_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "omegaMin",
+            RASDict_,
+            SMALL,
+            dimless/dimTime
+        )
+    )
+{
     // Force the construction of the mesh deltaCoeffs which may be needed
     // for the construction of the derived models and BCs
     this->mesh_.deltaCoeffs();

src/TurbulenceModels/turbulenceModels/RAS/RASModel/RASModel.H

@@ -42,14 +42,6 @@ SourceFiles
 #define RASModel_H
 
 #include "TurbulenceModel.H"
-#include "volFields.H"
-#include "surfaceFields.H"
-#include "fvm.H"
-#include "fvc.H"
-#include "fvMatrices.H"
-#include "IOdictionary.H"
-#include "Switch.H"
-#include "bound.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -178,6 +170,10 @@ public:
 
     // Member Functions
 
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+
         // Access
 
             //- Return the lower allowable limit for k (default: SMALL)
@@ -244,9 +240,6 @@ public:
 
         //- Solve the turbulence equations and correct the turbulence viscosity
         virtual void correct();
-
-        //- Read RASProperties dictionary
-        virtual bool read();
 };