diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
index e8c2ce6d54..4671347b66 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@@ -1,10 +1,10 @@
 Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<psiThermo> pThermo
+autoPtr<fluidThermo> pThermo
 (
-    psiThermo::New(mesh)
+    fluidThermo::New(mesh)
 );
-psiThermo& thermo = pThermo();
+fluidThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
 volScalarField rho
@@ -38,35 +38,10 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-
-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, simple.dict(), pRefCell, pRefValue);
+pressureControl pressureControl(p, rho, simple.dict());
 
 mesh.setFluxRequired(p.name());
 
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        simple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        simple.dict(),
-        dimDensity,
-        0
-    )
-);
-
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index e46f2a6691..06fa890c05 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -30,7 +30,11 @@
             // Relax the pressure equation to ensure diagonal-dominance
             pEqn.relax();
 
-            pEqn.setReference(pRefCell, pRefValue);
+            pEqn.setReference
+            (
+                pressureControl.refCell(),
+                pressureControl.refValue()
+            );
 
             pEqn.solve();
 
@@ -60,7 +64,11 @@
                 fvOptions(psi, p, rho.name())
             );
 
-            pEqn.setReference(pRefCell, pRefValue);
+            pEqn.setReference
+            (
+                pressureControl.refCell(),
+                pressureControl.refValue()
+            );
 
             pEqn.solve();
 
@@ -81,6 +89,8 @@
     U.correctBoundaryConditions();
     fvOptions.correct(U);
 
+    pressureControl.limit(p);
+
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
@@ -89,16 +99,12 @@
             /fvc::domainIntegrate(psi);
     }
 
+    p.correctBoundaryConditions();
+
     rho = thermo.rho();
-    rho = max(rho, rhoMin);
-    rho = min(rho, rhoMax);
 
     if (!simple.transonic())
     {
         rho.relax();
     }
-
-    Info<< "rho max/min : "
-        << max(rho).value() << " "
-        << min(rho).value() << endl;
 }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
index 67b962f805..0dbde82960 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
@@ -40,7 +40,11 @@ if (simple.transonic())
         // Relax the pressure equation to maintain diagonal dominance
         pEqn.relax();
 
-        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
 
         pEqn.solve();
 
@@ -75,7 +79,11 @@ else
             fvOptions(psi, p, rho.name())
         );
 
-        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
 
         pEqn.solve();
 
@@ -97,6 +105,8 @@ U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 
+pressureControl.limit(p);
+
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
 if (closedVolume)
@@ -105,14 +115,12 @@ if (closedVolume)
         /fvc::domainIntegrate(psi);
 }
 
+p.correctBoundaryConditions();
+
 // Recalculate density from the relaxed pressure
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
 
 if (!simple.transonic())
 {
     rho.relax();
 }
-
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H
index 5daef0c2de..4671347b66 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H
@@ -1,10 +1,10 @@
 Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<rhoThermo> pThermo
+autoPtr<fluidThermo> pThermo
 (
-    rhoThermo::New(mesh)
+    fluidThermo::New(mesh)
 );
-rhoThermo& thermo = pThermo();
+fluidThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
 volScalarField rho
@@ -38,46 +38,10 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-
-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, simple.dict(), pRefCell, pRefValue);
+pressureControl pressureControl(p, rho, simple.dict());
 
 mesh.setFluxRequired(p.name());
 
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        simple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        simple.dict(),
-        dimDensity,
-        0
-    )
-);
-
-dimensionedScalar pMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "pMin",
-        simple.dict(),
-        dimPressure,
-        10
-    )
-);
-
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H
index ca8725b499..4167f80e59 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H
@@ -48,7 +48,11 @@
 
         fvScalarMatrix& pEqn = tpEqn.ref();
 
-        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
 
         pEqn.solve();
 
@@ -75,6 +79,8 @@
     U.correctBoundaryConditions();
     fvOptions.correct(U);
 
+    pressureControl.limit(p);
+
     // For closed-volume cases adjust the pressure and density levels
     // to obey overall mass continuity
     if (closedVolume)
@@ -83,13 +89,6 @@
             /fvc::domainIntegrate(psi);
     }
 
-    p = max(p, pMin);
-
     rho = thermo.rho();
-    rho = max(rho, rhoMin);
-    rho = min(rho, rhoMax);
     rho.relax();
-    Info<< "rho max/min : "
-        << max(rho).value() << " "
-        << min(rho).value() << endl;
 }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
index 6cca40e2c0..29839274e8 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,11 +31,12 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "rhoThermo.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
+#include "simpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "IOporosityModelList.H"
-#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
index 4a3cdf8e01..200e620eef 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,9 +30,10 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "psiThermo.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "simpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/etc/templates/compressibleInflowOutflow/system/fvSchemes b/etc/templates/compressibleInflowOutflow/system/fvSchemes
index 9986456cf7..ad7ca94ca6 100644
--- a/etc/templates/compressibleInflowOutflow/system/fvSchemes
+++ b/etc/templates/compressibleInflowOutflow/system/fvSchemes
@@ -44,9 +44,9 @@ divSchemes
     div(phi,e)      $turbulence;
     div(phi,h)      $turbulence;
     div(phi,K)      $turbulence;
-    div(phi,Ekp)    $turbulence; 
+    div(phi,Ekp)    $turbulence;
 
-    div(phid,p)     bounded Gauss upwind;
+    div(phid,p)     Gauss upwind;
     div((phi|interpolate(rho)),p)  bounded Gauss upwind;
 
     div(((rho*nuEff)*dev2(T(grad(U)))))    Gauss linear;
diff --git a/etc/templates/compressibleInflowOutflow/system/fvSolution b/etc/templates/compressibleInflowOutflow/system/fvSolution
index 80579478d3..d20f03588b 100644
--- a/etc/templates/compressibleInflowOutflow/system/fvSolution
+++ b/etc/templates/compressibleInflowOutflow/system/fvSolution
@@ -43,8 +43,8 @@ SIMPLE
     }
 
     nNonOrthogonalCorrectors 0;
-    rhoMin          rhoMin [ 1 -3 0 0 0 ] 0.1;
-    rhoMax          rhoMax [ 1 -3 0 0 0 ] 1.5;
+    pMinFactor      0.1;
+    pMaxFactor      1.5;
 }
 
 relaxationFactors
@@ -56,7 +56,7 @@ relaxationFactors
     }
     equations
     {
-     	U               0.7;
+        U               0.7;
         "(e|h)"         0.7;
         "(k|epsilon|omega)" 0.7;
     }
diff --git a/src/TurbulenceModels/compressible/CompressibleTurbulenceModel/CompressibleTurbulenceModel.H b/src/TurbulenceModels/compressible/CompressibleTurbulenceModel/CompressibleTurbulenceModel.H
index 763303447e..07b32950d7 100644
--- a/src/TurbulenceModels/compressible/CompressibleTurbulenceModel/CompressibleTurbulenceModel.H
+++ b/src/TurbulenceModels/compressible/CompressibleTurbulenceModel/CompressibleTurbulenceModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,6 +115,19 @@ public:
             return this->transport_.mu(patchi);
         }
 
+        //- Return the laminar viscosity
+        virtual tmp<volScalarField> nu() const
+        {
+            return this->transport_.mu()/this->rho_;
+        }
+
+        //- Return the laminar viscosity on patchi
+        virtual tmp<scalarField> nu(const label patchi) const
+        {
+            return
+                this->transport_.mu(patchi)/this->rho_.boundaryField()[patchi];
+        }
+
         //- Return the turbulence dynamic viscosity
         virtual tmp<volScalarField> mut() const
         {
diff --git a/src/TurbulenceModels/turbulenceModels/TurbulenceModel/TurbulenceModel.H b/src/TurbulenceModels/turbulenceModels/TurbulenceModel/TurbulenceModel.H
index bff236d486..61e0ac3d7f 100644
--- a/src/TurbulenceModels/turbulenceModels/TurbulenceModel/TurbulenceModel.H
+++ b/src/TurbulenceModels/turbulenceModels/TurbulenceModel/TurbulenceModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -158,13 +158,13 @@ public:
         }
 
         //- Return the laminar viscosity
-        tmp<volScalarField> nu() const
+        virtual tmp<volScalarField> nu() const
         {
             return transport_.nu();
         }
 
         //- Return the laminar viscosity on patchi
-        tmp<scalarField> nu(const label patchi) const
+        virtual tmp<scalarField> nu(const label patchi) const
         {
             return transport_.nu(patchi);
         }
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 2d7fbce130..8f9122b58e 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -402,6 +402,7 @@ $(general)/constrainHbyA/constrainHbyA.C
 $(general)/adjustPhi/adjustPhi.C
 $(general)/bound/bound.C
 $(general)/CorrectPhi/correctUphiBCs.C
+$(general)/pressureControl/pressureControl.C
 
 solutionControl = $(general)/solutionControl
 $(solutionControl)/solutionControl/solutionControl.C
diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
index 0d24e83a7b..136fdda595 100644
--- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
+++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,7 +27,7 @@ License
 
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
-void Foam::setRefCell
+bool Foam::setRefCell
 (
     const volScalarField& field,
     const volScalarField& fieldRef,
@@ -109,11 +109,17 @@ void Foam::setRefCell
         }
 
         refValue = readScalar(dict.lookup(refValueName));
+
+        return true;
+    }
+    else
+    {
+        return false;
     }
 }
 
 
-void Foam::setRefCell
+bool Foam::setRefCell
 (
     const volScalarField& field,
     const dictionary& dict,
@@ -122,7 +128,7 @@ void Foam::setRefCell
     const bool forceReference
 )
 {
-    setRefCell(field, field, dict, refCelli, refValue, forceReference);
+    return setRefCell(field, field, dict, refCelli, refValue, forceReference);
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H
index 8a9d3d100a..1532c718e7 100644
--- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H
+++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H
@@ -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-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,8 +46,8 @@ namespace Foam
 
 //- If the field fieldRef needs referencing find the reference cell nearest
 //  (in index) to the given cell looked-up for field, but which is not on a
-//  cyclic, symmetry or processor patch.
-void setRefCell
+//  cyclic, symmetry or processor patch and return true, otherwise return false.
+bool setRefCell
 (
     const volScalarField& field,
     const volScalarField& fieldRef,
@@ -59,8 +59,8 @@ void setRefCell
 
 //- If the field needs referencing find the reference cell nearest
 //  (in index) to the given cell looked-up for field, but which is not on a
-//  cyclic, symmetry or processor patch.
-void setRefCell
+//  cyclic, symmetry or processor patch and return true, otherwise return false.
+bool setRefCell
 (
     const volScalarField& field,
     const dictionary& dict,
diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
new file mode 100644
index 0000000000..101126a5e8
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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 "pressureControl.H"
+#include "findRefCell.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::pressureControl::pressureControl
+(
+    const volScalarField& p,
+    const volScalarField& rho,
+    const dictionary& dict
+)
+:
+    refCell_(0),
+    refValue_(0),
+    pMax_("pMax", dimPressure, 0),
+    pMin_("pMin", dimPressure, GREAT)
+{
+    if (setRefCell(p, dict, refCell_, refValue_))
+    {
+        pMax_.value() = refValue_;
+        pMin_.value() = refValue_;
+    }
+
+    const volScalarField::Boundary& pbf = p.boundaryField();
+    const volScalarField::Boundary& rhobf = rho.boundaryField();
+
+    scalar rhoRefMax = -GREAT;
+    scalar rhoRefMin = GREAT;
+    bool rhoLimits = false;
+
+    forAll(pbf, patchi)
+    {
+        if (pbf[patchi].fixesValue())
+        {
+            rhoLimits = true;
+
+            pMax_.value() = max(pMax_.value(), max(pbf[patchi]));
+            pMin_.value() = min(pMin_.value(), min(pbf[patchi]));
+
+            rhoRefMax = max(rhoRefMax, max(rhobf[patchi]));
+            rhoRefMin = min(rhoRefMin, min(rhobf[patchi]));
+        }
+    }
+
+    if (dict.found("pMax"))
+    {
+        pMax_.value() = readScalar(dict.lookup("pMax"));
+    }
+    else if (dict.found("pMaxFactor"))
+    {
+        const scalar pMaxFactor(readScalar(dict.lookup("pMaxFactor")));
+        pMax_ *= pMaxFactor;
+    }
+    else if (dict.found("rhoMax"))
+    {
+        // For backward-compatibility infer the pMax from rhoMax
+
+        IOWarningInFunction(dict)
+            << "'rhoMax' specified rather than 'pMax' or 'pMaxFactor'" << nl
+            << "    This is supported for backward-compatibility but "
+               "'pMax' or 'pMaxFactor' are more reliable." << endl;
+
+        if (!rhoLimits)
+        {
+            FatalIOErrorInFunction(dict)
+                << "'rhoMax' specified rather than 'pMaxFactor'" << nl
+                << "    but the corresponding reference density cannot"
+                   " be evaluated from the boundary conditions." << nl
+                << "Please specify 'pMaxFactor' rather than 'rhoMax'"
+                << exit(FatalError);
+        }
+
+        dimensionedScalar rhoMax("rhoMax", dimDensity, dict);
+
+        pMax_ *= max(rhoMax.value()/rhoRefMax, 1);
+    }
+
+    if (dict.found("pMin"))
+    {
+        pMin_.value() = readScalar(dict.lookup("pMin"));
+    }
+    else if (dict.found("pMinFactor"))
+    {
+        const scalar pMinFactor(readScalar(dict.lookup("pMinFactor")));
+        pMin_ *= pMinFactor;
+    }
+    else if (dict.found("rhoMin"))
+    {
+        // For backward-compatibility infer the pMin from rhoMin
+
+        IOWarningInFunction(dict)
+            << "'rhoMin' specified rather than 'pMin' or 'pMinFactor'" << nl
+            << "    This is supported for backward-compatibility but"
+               "'pMin' or 'pMinFactor' are more reliable." << endl;
+
+        if (!rhoLimits)
+        {
+            FatalIOErrorInFunction(dict)
+                << "'rhoMin' specified rather than 'pMinFactor'" << nl
+                << "    but the corresponding reference density cannot"
+                   " be evaluated from the boundary conditions." << nl
+                << "Please specify 'pMinFactor' rather than 'rhoMin'"
+                << exit(FatalError);
+        }
+
+        dimensionedScalar rhoMin("rhoMin", dimDensity, dict);
+
+        pMin_ *= min(rhoMin.value()/rhoRefMin, 1);
+    }
+
+    Info<< "pressureControl" << nl
+        << "    pMax/pMin " << pMax_.value() << " " << pMin_.value()
+        << nl << endl;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::pressureControl::limit(volScalarField& p) const
+{
+    Info<< "pressureControl: p max/min "
+        << max(p).value() << " "
+        << min(p).value() << endl;
+
+    p = max(p, pMin_);
+    p = min(p, pMax_);
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
new file mode 100644
index 0000000000..19bd053d43
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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::pressureControl
+
+Description
+    Provides controls for the pressure reference is closed-volume simulations
+    and a general method for limiting the pressure during the startup of
+    steady-state simulations.
+
+SourceFiles
+    pressureControlI.H
+    pressureControl.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pressureControl_H
+#define pressureControl_H
+
+#include "dimensionedScalar.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class pressureControl Declaration
+\*---------------------------------------------------------------------------*/
+
+class pressureControl
+{
+    // Private data
+
+        //- Optional cell in which the reference pressure is set
+        label refCell_;
+
+        //- Optional pressure reference level
+        scalar refValue_;
+
+        //- Pressure lower-limit
+        dimensionedScalar pMax_;
+
+        //- Pressure upper-limit
+        dimensionedScalar pMin_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from the SIMPLE/PIMPLE sub-dictionary
+        pressureControl
+        (
+            const volScalarField& p,
+            const volScalarField& rho,
+            const dictionary& dict
+        );
+
+
+    // Member Functions
+
+        //- Return the cell in which the reference pressure is set
+        inline label refCell() const;
+
+        //- Return the pressure reference level
+        inline scalar refValue() const;
+
+        //- Limit the pressure
+        void limit(volScalarField& p) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pressureControlI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControlI.H b/src/finiteVolume/cfdTools/general/pressureControl/pressureControlI.H
new file mode 100644
index 0000000000..6548d43455
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControlI.H
@@ -0,0 +1,40 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::label Foam::pressureControl::refCell() const
+{
+    return refCell_;
+}
+
+
+inline Foam::scalar Foam::pressureControl::refValue() const
+{
+    return refValue_;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSolution b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSolution
index 9eee4c23fb..4b7b15b927 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSolution
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSolution
@@ -54,8 +54,8 @@ solvers
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.5;
-    rhoMax          1.5;
+    pMinFactor      0.5;
+    pMaxFactor      1.5;
 
     residualControl
     {
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSolution b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSolution
index 8699fd7095..421197bf13 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSolution
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSolution
@@ -48,8 +48,8 @@ SIMPLE
 {
     nUCorrectors        2;
     nNonOrthogonalCorrectors 0;
-    rhoMin              0.5;
-    rhoMax              2.0;
+    pMinFactor          0.5;
+    pMaxFactor          2.0;
 
     residualControl
     {
diff --git a/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSolution b/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSolution
index a054853c3e..25619f1f3f 100644
--- a/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSolution
+++ b/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSolution
@@ -54,8 +54,8 @@ solvers
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.4;
-    rhoMax          1.5;
+    pMinFactor      0.4;
+    pMaxFactor      1.5;
 
     residualControl
     {
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution
index 3855d4cb05..2ffe1fe4be 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution
@@ -39,8 +39,8 @@ solvers
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.1;
-    rhoMax          1.0;
+    pMinFactor      0.1;
+    pMaxFactor      1.5;
     transonic       yes;
     consistent      yes;
 
@@ -60,7 +60,6 @@ relaxationFactors
     fields
     {
         p               1;
-        rho             1;
     }
     equations
     {
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/T b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/T
new file mode 100644
index 0000000000..2df93ba37e
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/T
@@ -0,0 +1,42 @@
+/*--------------------------------*- 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       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            fixedValue;
+        value           uniform 350;
+    }
+
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/U b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/U
new file mode 100644
index 0000000000..99e3c1d635
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/U
@@ -0,0 +1,42 @@
+/*--------------------------------*- 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       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            noSlip;
+    }
+    inlet
+    {
+        type            flowRateInletVelocity;
+        massFlowRate    constant 5;
+        rhoInlet        1000;    // Guess for rho
+    }
+    outlet
+    {
+        type            inletOutlet;
+        value           uniform (0 0 0);
+        inletValue      uniform (0 0 0);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/alphat b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/alphat
new file mode 100644
index 0000000000..3a6625042c
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/alphat
@@ -0,0 +1,43 @@
+/*--------------------------------*- 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       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/epsilon b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/epsilon
new file mode 100644
index 0000000000..71b9c782d1
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/epsilon
@@ -0,0 +1,47 @@
+/*--------------------------------*- 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       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 200;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 200;
+    }
+    inlet
+    {
+        type            turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.005;
+        value           uniform 200;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 200;
+        value           uniform 200;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/k b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/k
new file mode 100644
index 0000000000..7ee52f8dbf
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/k
@@ -0,0 +1,44 @@
+/*--------------------------------*- 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       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            kqRWallFunction;
+        value           uniform 1;
+    }
+    inlet
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.05;
+        value           uniform 1;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/nut b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/nut
new file mode 100644
index 0000000000..ac8d562f40
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/nut
@@ -0,0 +1,45 @@
+/*--------------------------------*- 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       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/p b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/p
new file mode 100644
index 0000000000..5e0d23f65f
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/0/p
@@ -0,0 +1,39 @@
+/*--------------------------------*- 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       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            zeroGradient;
+    }
+    inlet
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/thermophysicalProperties b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/thermophysicalProperties
new file mode 100644
index 0000000000..426ab72116
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/thermophysicalProperties
@@ -0,0 +1,31 @@
+/*--------------------------------*- 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      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    properties      liquid;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    H2O;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/turbulenceProperties b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/turbulenceProperties
new file mode 100644
index 0000000000..cd2daf8229
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/turbulenceProperties
@@ -0,0 +1,30 @@
+/*--------------------------------*- 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      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RAS;
+
+RAS
+{
+    RASModel        kEpsilon;
+
+    turbulence      on;
+
+    printCoeffs     on;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/blockMeshDict b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/blockMeshDict
new file mode 100644
index 0000000000..9bec6a9547
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/blockMeshDict
@@ -0,0 +1,127 @@
+/*--------------------------------*- 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.001;
+
+vertices
+(
+    // front-plane: z = +25mm
+    // inlet region
+    (   -50  25   25)           // pt 0
+    (     0  25   25)           // pt 1
+    (   -50  75   25)           // pt 2
+    (     0  75   25)           // pt 3
+    // outlet region
+    (  -500 -75   25)           // pt 4
+    (     0 -75   25)           // pt 5
+    (  -500 -25   25)           // pt 6
+    (     0 -25   25)           // pt 7
+    // bend mid-points
+    (    25   0   25)           // pt 8
+    (    75   0   25)           // pt 9
+    // back-plane: z = -25mm
+    // inlet region
+    (   -50  25   -25)          // pt 0 + 10
+    (     0  25   -25)          // pt 1 + 10
+    (   -50  75   -25)          // pt 2 + 10
+    (     0  75   -25)          // pt 3 + 10
+    // outlet region
+    (  -500 -75   -25)          // pt 4 + 10
+    (     0 -75   -25)          // pt 5 + 10
+    (  -500 -25   -25)          // pt 7 + 10
+    (     0 -25   -25)          // pt 8 + 10
+    // bend mid-points
+    (    25   0   -25)          // pt 8 + 10
+    (    75   0   -25)          // pt 9 + 10
+);
+
+blocks
+(
+    hex (0 1 11 10  2 3 13 12) inlet  ( 20 20 20)  simpleGrading (1 1 1)
+    hex (4 5 15 14  6 7 17 16) outlet (200 20 20)  simpleGrading (1 1 1)
+
+    hex (1 8 18 11  3 9 19 13) bend1  ( 30 20 20)  simpleGrading (1 1 1)
+    hex (5 9 19 15  7 8 18 17) bend2  ( 30 20 20)  simpleGrading (1 1 1)
+);
+
+edges
+(
+   // block 2
+   arc  1  8  ( 17.678  17.678  25)
+   arc 11 18  ( 17.678  17.678 -25)
+   arc  3  9  ( 53.033  53.033  25)
+   arc 13 19  ( 53.033  53.033 -25)
+   // block 3
+   arc  7  8  ( 17.678  -17.678  25)
+   arc 17 18  ( 17.678  -17.678 -25)
+   arc  5  9  ( 53.033  -53.033  25)
+   arc 15 19  ( 53.033  -53.033 -25)
+);
+
+boundary
+(
+        // is there no way of defining all my 'defaultFaces' to be 'wall'?
+    Default_Boundary_Region
+    {
+        type wall;
+        faces
+        (
+            // block0
+            ( 0 1 3 2 )
+            ( 11 10 12 13 )
+            ( 0 10 11 1 )
+            ( 2 3 13 12 )
+            // block1
+            ( 4 5 7 6 )
+            ( 15 14 16 17 )
+            ( 4 14 15 5 )
+            ( 6 7 17 16 )
+            // block2
+            ( 1 8 9 3 )
+            ( 18 11 13 19 )
+            ( 3 9 19 13 )
+            ( 1 11 18 8 )
+            // block3
+            ( 5 9 8 7 )
+            ( 19 15 17 18 )
+            ( 5 15 19 9 )
+            ( 7 8 18 17 )
+        );
+    }
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 2 12 10)
+        );
+    }
+
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (4 6 16 14)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/controlDict b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/controlDict
new file mode 100644
index 0000000000..edddd719e6
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/controlDict
@@ -0,0 +1,51 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     rhoSimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         500;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   100;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+graphFormat     raw;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/decomposeParDict b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/decomposeParDict
new file mode 100644
index 0000000000..9e844a62ba
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/decomposeParDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- 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      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 8;
+
+method          hierarchical;
+
+simpleCoeffs
+{
+    n               (8 1 1);
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               (4 2 1);
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSchemes b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSchemes
new file mode 100644
index 0000000000..c68f548c95
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSchemes
@@ -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      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default             steadyState;
+}
+
+gradSchemes
+{
+    default             Gauss linear;
+
+    limited             cellLimited Gauss linear 1;
+}
+
+divSchemes
+{
+    default             none;
+
+    div(phi,U)          bounded Gauss linearUpwind limited;
+    div(phi,e)          bounded Gauss linearUpwind limited;
+    div(phi,epsilon)    bounded Gauss linearUpwind limited;
+    div(phi,k)          bounded Gauss linearUpwind limited;
+    div(phi,Ekp)        bounded Gauss linearUpwind limited;
+
+    div(((rho*nuEff)*dev2(T(grad(U)))))      Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSolution b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSolution
new file mode 100644
index 0000000000..3b30d4ebe5
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSolution
@@ -0,0 +1,71 @@
+/*--------------------------------*- 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          GAMG;
+        smoother        GaussSeidel;
+
+        tolerance       1e-7;
+        relTol          0.1;
+    }
+
+    "(U|e|k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+
+        tolerance       1e-7;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    pMinFactor      0.1;
+    pMaxFactor      1.5;
+
+    transonic       no;
+    consistent      no;
+
+    residualControl
+    {
+        p               1e-4;
+        U               1e-4;
+        "(k|omega|epsilon|e|h)" 1e-3;
+    }
+}
+
+relaxationFactors
+{
+    fields
+    {
+        p               0.3;
+    }
+    equations
+    {
+        U               0.7;
+        e               0.7;
+        k               0.7;
+        epsilon         0.7;
+    }
+}
+
+// ************************************************************************* //