Merge remote branch 'OpenCFD/olesenm' into olesenm

src/turbulenceModels/compressible/LES/DeardorffDiffStress/DeardorffDiffStress.C

@@ -106,7 +106,7 @@ void DeardorffDiffStress::correct(const tmp<volTensorField>& tgradU)
 
     volScalarField K = 0.5*tr(B_);
 
-    solve
+    tmp<fvSymmTensorMatrix> BEqn
     (
         fvm::ddt(rho(), B_)
       + fvm::div(phi(), B_)
@@ -118,6 +118,8 @@ void DeardorffDiffStress::correct(const tmp<volTensorField>& tgradU)
       - (2*ce_ - 0.667*cm_)*I*rho()*epsilon()
     );
 
+    BEqn().relax();
+    BEqn().solve();
 
     // Bounding the component kinetic energies
 

src/turbulenceModels/compressible/LES/GenEddyVisc/GenEddyVisc.C

@@ -108,7 +108,9 @@ GenEddyVisc::GenEddyVisc
         ),
         mesh_
     )
-{}
+{
+    bound(k_, kMin_);
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

src/turbulenceModels/compressible/LES/LESModel/LESModel.C

@@ -79,7 +79,7 @@ LESModel::LESModel
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(subOrEmptyDict(type + "Coeffs")),
 
-    kMin_("kMin", dimVelocity*dimVelocity, SMALL),
+    kMin_("kMin", sqr(dimVelocity), SMALL),
 
     delta_(LESdelta::New("delta", U.mesh(), *this))
 {
@@ -131,8 +131,12 @@ autoPtr<LESModel> LESModel::New
     {
         FatalErrorIn
         (
-            "LESModel::New(const volVectorField& U, const "
+            "LESModel::New"
-            "surfaceScalarField& phi, const basicThermo&)"
+            "("
+                "const volVectorField&, "
+                "const surfaceScalarField&, "
+                "const basicThermo&"
+            ")"
         )   << "Unknown LESModel type " << modelName
             << endl << endl
             << "Valid LESModel types are :" << endl
@@ -167,7 +171,7 @@ bool LESModel::read()
             coeffDict_ <<= *dictPtr;
         }
 
-        readIfPresent("kMin", kMin_);
+        kMin_.readIfPresent(*this);
 
         delta_().read(*this);
 

src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.C

@@ -302,7 +302,7 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& tgradU)
     volScalarField Stilda =
         fv3()*::sqrt(2.0)*mag(skew(gradU)) + fv2()*nuTilda_/sqr(kappa_*dTilda_);
 
-    solve
+    tmp<fvScalarMatrix> nuTildaEqn
     (
         fvm::ddt(rho(), nuTilda_)
       + fvm::div(phi(), nuTilda_)
@@ -318,6 +318,9 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& tgradU)
       - fvm::Sp(rho()*Cw1_*fw(Stilda)*nuTilda_/sqr(dTilda_), nuTilda_)
     );
 
+    nuTildaEqn().relax();
+    nuTildaEqn().solve();
+
     bound(nuTilda_, dimensionedScalar("zero", nuTilda_.dimensions(), 0.0));
     nuTilda_.correctBoundaryConditions();
 

src/turbulenceModels/compressible/LES/dynOneEqEddy/dynOneEqEddy.C

@@ -118,7 +118,7 @@ void dynOneEqEddy::correct(const tmp<volTensorField>& tgradU)
     volScalarField divU = fvc::div(phi()/fvc::interpolate(rho()));
     volScalarField G = 2*muSgs_*(gradU && D);
 
-    solve
+    tmp<fvScalarMatrix> kEqn
     (
         fvm::ddt(rho(), k_)
       + fvm::div(phi(), k_)
@@ -129,9 +129,8 @@ void dynOneEqEddy::correct(const tmp<volTensorField>& tgradU)
       - fvm::Sp(ce_(D)*rho()*sqrt(k_)/delta(), k_)
     );
 
-    //FIXME: why not this?
+    kEqn().relax();
-    // kEqn.relax();
+    kEqn().solve();
-    // kEqn.solve();
 
     bound(k_, kMin_);
 

src/turbulenceModels/compressible/LES/lowReOneEqEddy/lowReOneEqEddy.C

@@ -106,7 +106,7 @@ void lowReOneEqEddy::correct(const tmp<volTensorField>& tgradU)
     volScalarField divU = fvc::div(phi()/fvc::interpolate(rho()));
     volScalarField G = 2*muSgs_*(gradU && dev(symm(gradU)));
 
-    solve
+    tmp<fvScalarMatrix> kEqn
     (
         fvm::ddt(rho(), k_)
       + fvm::div(phi(), k_)
@@ -117,6 +117,9 @@ void lowReOneEqEddy::correct(const tmp<volTensorField>& tgradU)
       - fvm::Sp(ce_*rho()*sqrt(k_)/delta(), k_)
     );
 
+    kEqn().relax();
+    kEqn().solve();
+
     bound(k_, kMin_);
 
     updateSubGridScaleFields();

src/turbulenceModels/compressible/LES/oneEqEddy/oneEqEddy.C

@@ -93,7 +93,7 @@ void oneEqEddy::correct(const tmp<volTensorField>& tgradU)
     volScalarField divU = fvc::div(phi()/fvc::interpolate(rho()));
     volScalarField G = 2*muSgs_*(gradU && dev(symm(gradU)));
 
-    fvScalarMatrix kEqn
+    tmp<fvScalarMatrix> kEqn
     (
         fvm::ddt(rho(), k_)
       + fvm::div(phi(), k_)
@@ -104,8 +104,8 @@ void oneEqEddy::correct(const tmp<volTensorField>& tgradU)
       - fvm::Sp(ce_*rho()*sqrt(k_)/delta(), k_)
     );
 
-    kEqn.relax();
+    kEqn().relax();
-    kEqn.solve();
+    kEqn().solve();
 
     bound(k_, kMin_);
 

src/turbulenceModels/compressible/RAS/LRR/LRR.C

@@ -229,6 +229,7 @@ LRR::LRR
             << exit(FatalError);
     }
 
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     mut_ = Cmu_*rho_*sqr(k_)/epsilon_;

src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C

@@ -258,6 +258,7 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             << exit(FatalError);
     }
 
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     mut_ = Cmu_*rho_*sqr(k_)/epsilon_;

src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C

@@ -187,6 +187,7 @@ LaunderSharmaKE::LaunderSharmaKE
         autoCreateAlphat("alphat", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     mut_ = rho_*Cmu_*fMu()*sqr(k_)/epsilon_;

src/turbulenceModels/compressible/RAS/RASModel/RASModel.C

@@ -81,7 +81,7 @@ RASModel::RASModel
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(subOrEmptyDict(type + "Coeffs")),
 
-    kMin_("kMin", dimVelocity*dimVelocity, SMALL),
+    kMin_("kMin", sqr(dimVelocity), SMALL),
     epsilonMin_("epsilonMin", kMin_.dimensions()/dimTime, SMALL),
     omegaMin_("omegaMin", dimless/dimTime, SMALL),
 
@@ -133,9 +133,13 @@ autoPtr<RASModel> RASModel::New
     {
         FatalErrorIn
         (
-            "RASModel::New(const volScalarField&, "
+            "RASModel::New"
-            "const volVectorField&, const surfaceScalarField&, "
+            "("
-            "basicThermo&)"
+                "const volScalarField&, "
+                "const volVectorField&, "
+                "const surfaceScalarField&, "
+                "basicThermo&"
+            ")"
         )   << "Unknown RASModel type " << modelName
             << endl << endl
             << "Valid RASModel types are :" << endl

src/turbulenceModels/compressible/RAS/RNGkEpsilon/RNGkEpsilon.C

@@ -186,6 +186,7 @@ RNGkEpsilon::RNGkEpsilon
         autoCreateAlphat("alphat", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     mut_ = Cmu_*rho_*sqr(k_)/epsilon_;

src/turbulenceModels/compressible/RAS/kEpsilon/kEpsilon.C

@@ -168,6 +168,7 @@ kEpsilon::kEpsilon
         autoCreateAlphat("alphat", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     mut_ = Cmu_*rho_*sqr(k_)/epsilon_;

src/turbulenceModels/compressible/RAS/kOmegaSST/kOmegaSST.C

@@ -258,6 +258,7 @@ kOmegaSST::kOmegaSST
         autoCreateAlphat("alphat", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(omega_, omegaMin_);
 
     mut_ =

src/turbulenceModels/incompressible/LES/DeardorffDiffStress/DeardorffDiffStress.C

@@ -102,7 +102,7 @@ void DeardorffDiffStress::correct(const tmp<volTensorField>& tgradU)
     volScalarField K = 0.5*tr(B_);
     volScalarField Epsilon = 2*nuEff()*magSqr(D);
 
-    fvSymmTensorMatrix BEqn
+    tmp<fvSymmTensorMatrix> BEqn
     (
         fvm::ddt(B_)
       + fvm::div(phi(), B_)
@@ -114,8 +114,8 @@ void DeardorffDiffStress::correct(const tmp<volTensorField>& tgradU)
       - (2*ce_ - 0.667*cm_)*I*Epsilon
     );
 
-    BEqn.relax();
+    BEqn().relax();
-    BEqn.solve();
+    BEqn().solve();
 
     // Bounding the component kinetic energies
 

src/turbulenceModels/incompressible/LES/LESModel/LESModel.C

@@ -78,7 +78,7 @@ LESModel::LESModel
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(subOrEmptyDict(type + "Coeffs")),
 
-    kMin_("kMin", dimVelocity*dimVelocity, SMALL),
+    kMin_("kMin", sqr(dimVelocity), SMALL),
     delta_(LESdelta::New("delta", U.mesh(), *this))
 {
     readIfPresent("kMin", kMin_);
@@ -128,8 +128,12 @@ autoPtr<LESModel> LESModel::New
     {
         FatalErrorIn
         (
-            "LESModel::New(const volVectorField& U, const "
+            "LESModel::New"
-            "surfaceScalarField& phi, transportModel&)"
+            "("
+                "const volVectorField&, "
+                "const surfaceScalarField& ,"
+                "transportModel&"
+            ")"
         )   << "Unknown LESModel type " << modelName
             << endl << endl
             << "Valid LESModel types are :" << endl
@@ -167,7 +171,7 @@ bool LESModel::read()
 
         delta_().read(*this);
 
-        readIfPresent("kMin", kMin_);
+        kMin_.readIfPresent(*this);
 
         return true;
     }

src/turbulenceModels/incompressible/LES/LRRDiffStress/LRRDiffStress.C

@@ -90,7 +90,10 @@ LRRDiffStress::LRRDiffStress
         )
     )
 {
-    updateSubGridScaleFields(0.5*tr(B_));
+    volScalarField K = 0.5*tr(B_);
+    bound(K, kMin_);
+
+    updateSubGridScaleFields(K);
 
     printCoeffs();
 }
@@ -111,7 +114,7 @@ void LRRDiffStress::correct(const tmp<volTensorField>& tgradU)
     volScalarField K = 0.5*tr(B_);
     volScalarField Epsilon = 2*nuEff()*magSqr(D);
 
-    fvSymmTensorMatrix BEqn
+    tmp<fvSymmTensorMatrix> BEqn
     (
         fvm::ddt(B_)
       + fvm::div(phi(), B_)
@@ -124,8 +127,8 @@ void LRRDiffStress::correct(const tmp<volTensorField>& tgradU)
       - (0.667 - 2*c1_)*I*pow(K, 1.5)/delta()
     );
 
-    BEqn.relax();
+    BEqn().relax();
-    BEqn.solve();
+    BEqn().solve();
 
     // Bounding the component kinetic energies
 

src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C

@@ -295,7 +295,7 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
     const volScalarField dTilda = this->dTilda(S);
     const volScalarField STilda = this->STilda(S, dTilda);
 
-    fvScalarMatrix nuTildaEqn
+    tmp<fvScalarMatrix> nuTildaEqn
     (
         fvm::ddt(nuTilda_)
       + fvm::div(phi(), nuTilda_)
@@ -311,8 +311,8 @@ void SpalartAllmaras::correct(const tmp<volTensorField>& gradU)
       - fvm::Sp(Cw1_*fw(STilda, dTilda)*nuTilda_/sqr(dTilda), nuTilda_)
     );
 
-    nuTildaEqn.relax();
+    nuTildaEqn().relax();
-    nuTildaEqn.solve();
+    nuTildaEqn().solve();
 
     bound(nuTilda_, dimensionedScalar("zero", nuTilda_.dimensions(), 0.0));
     nuTilda_.correctBoundaryConditions();

src/turbulenceModels/incompressible/LES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C

@@ -176,7 +176,6 @@ SpalartAllmarasIDDES::SpalartAllmarasIDDES
             1.63
         )
     )
-
 {}
 
 

src/turbulenceModels/incompressible/LES/dynOneEqEddy/dynOneEqEddy.C

@@ -127,6 +127,8 @@ dynOneEqEddy::dynOneEqEddy
     filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
     filter_(filterPtr_())
 {
+    bound(k_, kMin_);
+
     updateSubGridScaleFields(symm(fvc::grad(U)));
 
     printCoeffs();
@@ -145,7 +147,7 @@ void dynOneEqEddy::correct(const tmp<volTensorField>& tgradU)
 
     volScalarField P = 2.0*nuSgs_*magSqr(D);
 
-    fvScalarMatrix kEqn
+    tmp<fvScalarMatrix> kEqn
     (
        fvm::ddt(k_)
      + fvm::div(phi(), k_)
@@ -155,8 +157,8 @@ void dynOneEqEddy::correct(const tmp<volTensorField>& tgradU)
      - fvm::Sp(ce(D)*sqrt(k_)/delta(), k_)
     );
 
-    kEqn.relax();
+    kEqn().relax();
-    kEqn.solve();
+    kEqn().solve();
 
     bound(k_, kMin_);
 

src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.C

@@ -118,6 +118,8 @@ dynSmagorinsky::dynSmagorinsky
     filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
     filter_(filterPtr_())
 {
+    bound(k_,  kMin_);
+
     updateSubGridScaleFields(dev(symm(fvc::grad(U))));
 
     printCoeffs();
@@ -133,6 +135,7 @@ void dynSmagorinsky::correct(const tmp<volTensorField>& gradU)
     volSymmTensorField D = dev(symm(gradU));
 
     k_ = cI(D)*sqr(delta())*magSqr(D);
+    bound(k_,  kMin_);
 
     updateSubGridScaleFields(D);
 }

src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C

@@ -323,6 +323,7 @@ kOmegaSSTSAS::kOmegaSSTSAS
         mesh_
     )
 {
+    bound(k_, kMin_);
     bound(omega_, omegaMin_);
 
     updateSubGridScaleFields(magSqr(symm(fvc::grad(U))));

src/turbulenceModels/incompressible/LES/locDynOneEqEddy/locDynOneEqEddy.C

@@ -119,6 +119,8 @@ locDynOneEqEddy::locDynOneEqEddy
     filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
     filter_(filterPtr_())
 {
+    bound(k_, kMin_);
+
     volScalarField KK = 0.5*(filter_(magSqr(U)) - magSqr(filter_(U)));
     updateSubGridScaleFields(symm(fvc::grad(U)), KK);
 
@@ -139,7 +141,7 @@ void locDynOneEqEddy::correct(const tmp<volTensorField>& gradU)
 
     volScalarField P = 2.0*nuSgs_*magSqr(D);
 
-    fvScalarMatrix kEqn
+    tmp<fvScalarMatrix> kEqn
     (
        fvm::ddt(k_)
      + fvm::div(phi(), k_)
@@ -149,8 +151,8 @@ void locDynOneEqEddy::correct(const tmp<volTensorField>& gradU)
      - fvm::Sp(ce(D, KK)*sqrt(k_)/delta(), k_)
     );
 
-    kEqn.relax();
+    kEqn().relax();
-    kEqn.solve();
+    kEqn().solve();
 
     bound(k_, kMin_);
 

src/turbulenceModels/incompressible/LES/oneEqEddy/oneEqEddy.C

@@ -86,6 +86,8 @@ oneEqEddy::oneEqEddy
         )
     )
 {
+    bound(k_, kMin_);
+
     updateSubGridScaleFields();
 
     printCoeffs();
@@ -100,7 +102,7 @@ void oneEqEddy::correct(const tmp<volTensorField>& gradU)
 
     volScalarField G = 2.0*nuSgs_*magSqr(symm(gradU));
 
-    fvScalarMatrix kEqn
+    tmp<fvScalarMatrix> kEqn
     (
        fvm::ddt(k_)
      + fvm::div(phi(), k_)
@@ -110,8 +112,8 @@ void oneEqEddy::correct(const tmp<volTensorField>& gradU)
      - fvm::Sp(ce_*sqrt(k_)/delta(), k_)
     );
 
-    kEqn.relax();
+    kEqn().relax();
-    kEqn.solve();
+    kEqn().solve();
 
     bound(k_, kMin_);
 

src/turbulenceModels/incompressible/LES/spectEddyVisc/spectEddyVisc.C

@@ -72,7 +72,6 @@ spectEddyVisc::spectEddyVisc
     LESModel(typeName, U, phi, transport),
     GenEddyVisc(U, phi, transport),
 
-
     cB_
     (
         dimensioned<scalar>::lookupOrAddToDict

src/turbulenceModels/incompressible/RAS/LRR/LRR.C

@@ -198,6 +198,7 @@ LRR::LRR
             << exit(FatalError);
     }
 
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*sqr(k_)/epsilon_;

src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C

@@ -140,6 +140,7 @@ LamBremhorstKE::LamBremhorstKE
         autoCreateLowReNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*fMu_*sqr(k_)/epsilon_;

src/turbulenceModels/incompressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C

@@ -224,6 +224,7 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
         autoCreateNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*sqr(k_)/epsilon_;

src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C

@@ -146,6 +146,7 @@ LaunderSharmaKE::LaunderSharmaKE
         autoCreateLowReNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilonTilda_, epsilonMin_);
 
     nut_ = Cmu_*fMu()*sqr(k_)/epsilonTilda_;

src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C

@@ -170,11 +170,17 @@ LienCubicKE::LienCubicKE
         autoCreateEpsilon("epsilon", mesh_)
     ),
 
-    //FIXME - epsilon is not bounded
-
     gradU_(fvc::grad(U)),
-    eta_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))),
+    eta_
-    ksi_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))),
+    (
+        k_/(epsilon_ + epsilonMin_)
+       *sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))
+    ),
+    ksi_
+    (
+        k_/(epsilon_ + epsilonMin_)
+       *sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))
+    ),
     Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
     fEta_(A2_ + pow(eta_, 3.0)),
 
@@ -228,6 +234,7 @@ LienCubicKE::LienCubicKE
         )
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*sqr(k_)/epsilon_ + C5viscosity_;

src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C

@@ -220,14 +220,22 @@ LienCubicKELowRe::LienCubicKELowRe
     y_(mesh_),
 
     gradU_(fvc::grad(U)),
-    eta_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))),
+    eta_
-    ksi_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))),
+    (
+        k_/(epsilon_ + epsilonMin_)
+       *sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))
+    ),
+    ksi_
+    (
+        k_/(epsilon_ + epsilonMin_)
+       *sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))
+    ),
     Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
-    fEta_(A2_ + pow(eta_, 3.0)),
+    fEta_(A2_ + pow3(eta_)),
 
     C5viscosity_
     (
-        -2.0*pow(Cmu_, 3.0)*pow(k_, 4.0)/pow(epsilon_, 3.0)
+        -2.0*pow3(Cmu_)*pow4(k_)/pow3(epsilon_ + epsilonMin_)
        *(magSqr(gradU_ + gradU_.T()) - magSqr(gradU_ - gradU_.T()))
     ),
 
@@ -252,7 +260,7 @@ LienCubicKELowRe::LienCubicKELowRe
         symm
         (
             // quadratic terms
-            pow(k_, 3.0)/sqr(epsilon_)
+            pow3(k_)/sqr(epsilon_ + epsilonMin_)
            *(
                 Ctau1_/fEta_
                *(
@@ -263,8 +271,8 @@ LienCubicKELowRe::LienCubicKELowRe
               + Ctau3_/fEta_*(gradU_.T() & gradU_)
             )
             // cubic term C4
-          - 20.0*pow(k_, 4.0)/pow(epsilon_, 3.0)
+          - 20.0*pow4(k_)/pow3(epsilon_ + epsilonMin_)
-           *pow(Cmu_, 3.0)
+           *pow3(Cmu_)
            *(
                 ((gradU_ & gradU_) & gradU_.T())
               + ((gradU_ & gradU_.T()) & gradU_.T())
@@ -280,6 +288,7 @@ LienCubicKELowRe::LienCubicKELowRe
         )
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_
@@ -430,8 +439,8 @@ void LienCubicKELowRe::correct()
 
     epsEqn().relax();
 
-#   include "LienCubicKELowReSetWallDissipation.H"
+    #include "LienCubicKELowReSetWallDissipation.H"
-#   include "wallDissipationI.H"
+    #include "wallDissipationI.H"
 
     solve(epsEqn);
     bound(epsilon_, epsilonMin_);

src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C

@@ -180,6 +180,7 @@ LienLeschzinerLowRe::LienLeschzinerLowRe
         autoCreateLowReNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*(scalar(1) - exp(-Am_*yStar_))

src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C

@@ -144,7 +144,6 @@ NonlinearKEShih::NonlinearKEShih
         )
     ),
 
-    //FIXME: should be named 'kappa_' or 'kappa'?
     kappa_
     (
         dimensioned<scalar>::lookupOrAddToDict
@@ -175,7 +174,6 @@ NonlinearKEShih::NonlinearKEShih
             IOobject::AUTO_WRITE
         ),
         mesh_
-        //FIXME: what about autoCreateK("k", mesh_)
     ),
 
     epsilon_
@@ -189,39 +187,30 @@ NonlinearKEShih::NonlinearKEShih
             IOobject::AUTO_WRITE
         ),
         mesh_
-        //FIXME: what about autoCreateK("epsilon", mesh_)
     ),
 
-    //FIXME: epsilon is not bounded
-
     gradU_(fvc::grad(U)),
-    eta_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))),
+    eta_
-    ksi_(k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))),
+    (
+        k_/(epsilon_ + epsilonMin_)
+       *sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))
+    ),
+    ksi_
+    (
+        k_/(epsilon_+ epsilonMin_)
+       *sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))
+    ),
     Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
     fEta_(A2_ + pow(eta_, 3.0)),
 
-    // FIXME: epsilon is not bounded
     nut_("nut", Cmu_*sqr(k_)/(epsilon_ + epsilonMin_)),
-    // FIXME: why not use the following?
-    // nut_
-    // (
-    //     IOobject
-    //     (
-    //         "nut",
-    //         runTime_.timeName(),
-    //         mesh_,
-    //         IOobject::NO_READ,
-    //         IOobject::AUTO_WRITE
-    //     ),
-    //     autoCreateNut("nut", mesh_)
-    // ),
 
     nonlinearStress_
     (
         "nonlinearStress",
         symm
         (
-            pow(k_, 3.0)/sqr(epsilon_)
+            pow3(k_)/sqr(epsilon_ + epsilonMin_)
            *(
                 Ctau1_/fEta_
                *(
@@ -234,11 +223,9 @@ NonlinearKEShih::NonlinearKEShih
         )
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
-    //FIXME: could use this
-    // nut_ = Cmu_*sqr(k_)/epsilon_;
-
     #include "wallNonlinearViscosityI.H"
 
     printCoeffs();

src/turbulenceModels/incompressible/RAS/RASModel/RASModel.C

@@ -80,7 +80,7 @@ RASModel::RASModel
     printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
     coeffDict_(subOrEmptyDict(type + "Coeffs")),
 
-    kMin_("kMin", dimVelocity*dimVelocity, SMALL),
+    kMin_("kMin", sqr(dimVelocity), SMALL),
     epsilonMin_("epsilonMin", kMin_.dimensions()/dimTime, SMALL),
     omegaMin_("omegaMin", dimless/dimTime, SMALL),
 
@@ -131,8 +131,12 @@ autoPtr<RASModel> RASModel::New
     {
         FatalErrorIn
         (
-            "RASModel::New(const volVectorField&, "
+            "RASModel::New"
-            "const surfaceScalarField&, transportModel&)"
+            "("
+                "const volVectorField&, "
+                "const surfaceScalarField&, "
+                "transportModel&"
+            ")"
         )   << "Unknown RASModel type " << modelName
             << endl << endl
             << "Valid RASModel types are :" << endl

src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C

@@ -155,6 +155,7 @@ RNGkEpsilon::RNGkEpsilon
         autoCreateNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*sqr(k_)/epsilon_;

src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C

@@ -128,6 +128,7 @@ kEpsilon::kEpsilon
         autoCreateNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
 
     nut_ = Cmu_*sqr(k_)/epsilon_;

src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C

@@ -137,6 +137,7 @@ kOmega::kOmega
         autoCreateNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(omega_, omegaMin_);
 
     nut_ = k_/omega_;

src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C

@@ -236,6 +236,7 @@ kOmegaSST::kOmegaSST
         autoCreateNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(omega_, omegaMin_);
 
     nut_ =

src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C

@@ -196,6 +196,7 @@ qZeta::qZeta
         autoCreateNut("nut", mesh_)
     )
 {
+    bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
     bound(q_, qMin_);
     bound(zeta_, zetaMin_);