From 870d71db1f57efc5b3edc5d687122c2c4867c651 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Wed, 28 Jan 2015 16:34:18 +0000
Subject: [PATCH] buoyantKEpsilon: filter sources depending on mag(g) so that
 it operates as kEpsilon if mag(g) = 0

---
 .../RAS/buoyantKEpsilon/buoyantKEpsilon.C     | 35 ++++++++++++++-----
 1 file changed, 27 insertions(+), 8 deletions(-)

diff --git a/src/TurbulenceModels/compressible/RAS/buoyantKEpsilon/buoyantKEpsilon.C b/src/TurbulenceModels/compressible/RAS/buoyantKEpsilon/buoyantKEpsilon.C
index 14b1b6ee6..2f7bdb161 100644
--- a/src/TurbulenceModels/compressible/RAS/buoyantKEpsilon/buoyantKEpsilon.C
+++ b/src/TurbulenceModels/compressible/RAS/buoyantKEpsilon/buoyantKEpsilon.C
@@ -116,7 +116,18 @@ template<class BasicTurbulenceModel>
 tmp<fvScalarMatrix>
 buoyantKEpsilon<BasicTurbulenceModel>::kSource() const
 {
-    return -fvm::SuSp(Gcoef(), this->k_);
+    const uniformDimensionedVectorField& g =
+        this->mesh_.objectRegistry::template
+        lookupObject<uniformDimensionedVectorField>("g");
+
+    if (mag(g.value()) > SMALL)
+    {
+        return -fvm::SuSp(Gcoef(), this->k_);
+    }
+    else
+    {
+        return kEpsilon<BasicTurbulenceModel>::kSource();
+    }
 }
 
 
@@ -128,15 +139,23 @@ buoyantKEpsilon<BasicTurbulenceModel>::epsilonSource() const
         this->mesh_.objectRegistry::template
         lookupObject<uniformDimensionedVectorField>("g");
 
-    vector gHat(g.value()/mag(g.value()));
+    if (mag(g.value()) > SMALL)
+    {
+        vector gHat(g.value()/mag(g.value()));
 
-    volScalarField v(gHat & this->U_);
-    volScalarField u
-    (
-        mag(this->U_ - gHat*v) + dimensionedScalar("SMALL", dimVelocity, SMALL)
-    );
+        volScalarField v(gHat & this->U_);
+        volScalarField u
+        (
+            mag(this->U_ - gHat*v)
+          + dimensionedScalar("SMALL", dimVelocity, SMALL)
+        );
 
-    return -fvm::SuSp(this->C1_*tanh(mag(v)/u)*Gcoef(), this->epsilon_);
+        return -fvm::SuSp(this->C1_*tanh(mag(v)/u)*Gcoef(), this->epsilon_);
+    }
+    else
+    {
+        return kEpsilon<BasicTurbulenceModel>::epsilonSource();
+    }
 }