From 3e03e3ba075ee59ec190b55cd5d5bc8cb535f4be Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Wed, 30 May 2012 15:23:56 +0100
Subject: [PATCH] MULES: corrected handling of density in implicitSolve

---
 .../fvMatrices/solvers/MULES/MULESTemplates.C | 44 ++++++++++++++++---
 1 file changed, 39 insertions(+), 5 deletions(-)

diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
index c16288ea35..c9dd621d3e 100644
--- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
+++ b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C
@@ -103,6 +103,29 @@ void Foam::MULES::explicitSolve
 }
 
 
+namespace Foam
+{
+namespace MULES
+{
+    template<class RhoType>
+    inline tmp<surfaceScalarField> interpolate(const RhoType& rho)
+    {
+        notImplemented
+        (
+            "tmp<surfaceScalarField> interpolate(const RhoType& rho)"
+        );
+        return tmp<surfaceScalarField>(NULL);
+    }
+
+    template<>
+    inline tmp<surfaceScalarField> interpolate(const volScalarField& rho)
+    {
+        return fvc::interpolate(rho);
+    }
+}
+}
+
+
 template<class RhoType, class SpType, class SuType>
 void Foam::MULES::implicitSolve
 (
@@ -143,10 +166,6 @@ void Foam::MULES::implicitSolve
     scalarField allCoLambda(mesh.nFaces());
 
     {
-        tmp<surfaceScalarField> Cof =
-            mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
-           *mag(phi)/mesh.magSf();
-
         slicedSurfaceScalarField CoLambda
         (
             IOobject
@@ -164,7 +183,22 @@ void Foam::MULES::implicitSolve
             false   // Use slices for the couples
         );
 
-        CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
+        if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
+        {
+            tmp<surfaceScalarField> Cof =
+                mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
+               *mag(phi/interpolate(rho))/mesh.magSf();
+
+            CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
+        }
+        else
+        {
+            tmp<surfaceScalarField> Cof =
+                mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
+               *mag(phi)/mesh.magSf();
+
+            CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
+        }
     }
 
     scalarField allLambda(allCoLambda);