From e98214b0ada18616662c40450a5273845db27da0 Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Tue, 27 Jun 2017 15:47:53 -0700
Subject: [PATCH] BUG: Updating rho at the end of the pEq for rhoPimpleDyMFoam
 and limiting is needed

---
 .../rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H           | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index ed802ba6d2c..d1dcc87b410 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -1,4 +1,11 @@
-rho = thermo.rho();
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -92,9 +99,11 @@ K = 0.5*magSqr(U);
 if (pressureControl.limit(p))
 {
     p.correctBoundaryConditions();
-    rho = thermo.rho();
 }
 
+thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
+rho = thermo.rho();
+
 {
     rhoUf = fvc::interpolate(rho*U);
     surfaceVectorField n(mesh.Sf()/mesh.magSf());
-- 
GitLab