From 2c6b405043e8a3133629868bf20c819854ed27aa Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Sat, 23 Apr 2016 10:04:39 +0100
Subject: [PATCH] fireFoam: Added optional hydrostatic initialization of the
 pressure and density

Also added the new prghTotalHydrostaticPressure p_rgh BC which uses the
hydrostatic pressure field as the reference state for the far-field
which provides much more accurate entrainment is large open domains
typical of many fire simulations.

The hydrostatic field solution is controlled by the optional entries in
the fvSolution.PIMPLE dictionary, e.g.

    hydrostaticInitialization yes;
    nHydrostaticCorrectors 5;

and the solver must also be specified for the hydrostatic p_rgh field
ph_rgh e.g.

    ph_rgh
    {
        $p_rgh;
    }

Suitable boundary conditions for ph_rgh cannot always be derived from
those for p_rgh and so the ph_rgh is read to provide them.

To avoid accuracy issues with IO, restart and post-processing the p_rgh
and ph_rgh the option to specify a suitable reference pressure is
provided via the optional pRef file in the constant directory, e.g.

    dimensions      [1 -1 -2 0 0 0 0];
    value           101325;

which is used in the relationship between p_rgh and p:

    p = p_rgh + rho*gh + pRef;

Note that if pRef is specified all pressure BC specifications in the
p_rgh and ph_rgh files are relative to the reference to avoid round-off
errors.

For examples of suitable BCs for p_rgh and ph_rgh for a range of
fireFoam cases please study the tutorials in
tutorials/combustion/fireFoam/les which have all been updated.

Henry G. Weller
CFD Direct Ltd.
---
 .../combustion/fireFoam/createFields.H        | 101 ++++-----
 .../solvers/combustion/fireFoam/fireFoam.C    |   1 -
 .../solvers/combustion/fireFoam/pEqn.H        |   9 +-
 .../solvers/combustion/fireFoam/phrghEqn.H    |  62 ++++++
 src/finiteVolume/Make/files                   |   1 +
 .../cfdTools/general/include/readpRef.H       |  13 ++
 .../basic/fixedValue/fixedValueFvPatchField.H |   2 +-
 ...talHydrostaticPressureFvPatchScalarField.C | 168 +++++++++++++++
 ...talHydrostaticPressureFvPatchScalarField.H | 202 ++++++++++++++++++
 .../prghTotalPressureFvPatchScalarField.H     |  20 +-
 .../flameSpreadWaterSuppressionPanel/0/p_rgh  |  20 +-
 .../flameSpreadWaterSuppressionPanel/Allclean |   2 +-
 .../flameSpreadWaterSuppressionPanel/Allrun   |   4 +-
 .../constant/pRef                             |  21 ++
 .../system/fvSolution                         |   8 +
 .../les/oppositeBurningPanels/0/p_rgh         |   9 +-
 .../les/oppositeBurningPanels/Allclean        |   1 +
 .../fireFoam/les/oppositeBurningPanels/Allrun |   2 +
 .../les/oppositeBurningPanels/constant/pRef   |  21 ++
 .../oppositeBurningPanels/system/fvSolution   |   9 +
 .../fireFoam/les/smallPoolFire2D/0/p_rgh      |   9 +-
 .../fireFoam/les/smallPoolFire2D/Allclean     |  11 +
 .../fireFoam/les/smallPoolFire2D/Allrun       |   2 +
 .../les/smallPoolFire2D/constant/pRef         |  21 ++
 .../les/smallPoolFire2D/system/fvSolution     |   8 +
 .../fireFoam/les/smallPoolFire3D/0/p_rgh      |   8 +-
 .../fireFoam/les/smallPoolFire3D/Allclean     |  11 +
 .../fireFoam/les/smallPoolFire3D/Allrun       |   3 +
 .../les/smallPoolFire3D/constant/pRef         |  21 ++
 .../les/smallPoolFire3D/system/fvSolution     |   8 +
 30 files changed, 673 insertions(+), 105 deletions(-)
 create mode 100644 applications/solvers/combustion/fireFoam/phrghEqn.H
 create mode 100644 src/finiteVolume/cfdTools/general/include/readpRef.H
 create mode 100644 src/finiteVolume/fields/fvPatchFields/derived/prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C
 create mode 100644 src/finiteVolume/fields/fvPatchFields/derived/prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.H
 create mode 100644 tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/constant/pRef
 create mode 100644 tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/pRef
 create mode 100755 tutorials/combustion/fireFoam/les/smallPoolFire2D/Allclean
 create mode 100644 tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/pRef
 create mode 100755 tutorials/combustion/fireFoam/les/smallPoolFire3D/Allclean
 create mode 100644 tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/pRef

diff --git a/applications/solvers/combustion/fireFoam/createFields.H b/applications/solvers/combustion/fireFoam/createFields.H
index c5fcc21c9..ce621258c 100644
--- a/applications/solvers/combustion/fireFoam/createFields.H
+++ b/applications/solvers/combustion/fireFoam/createFields.H
@@ -54,6 +54,9 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
+#include "createMRF.H"
+
+
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
@@ -69,6 +72,55 @@ autoPtr<compressible::turbulenceModel> turbulence
 // Set the turbulence into the combustion model
 combustion->setTurbulence(turbulence());
 
+
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+#include "readpRef.H"
+
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+mesh.setFluxRequired(p_rgh.name());
+
+#include "phrghEqn.H"
+
+
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+forAll(Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(thermo.he());
+
+IOdictionary additionalControlsDict
+(
+    IOobject
+    (
+        "additionalControls",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE
+    )
+);
+
+Switch solvePrimaryRegion
+(
+    additionalControlsDict.lookup("solvePrimaryRegion")
+);
+
 volScalarField dQ
 (
     IOobject
@@ -99,52 +151,3 @@ volScalarField dpdt
 
 Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
-
-
-#include "readGravitationalAcceleration.H"
-#include "readhRef.H"
-#include "gh.H"
-
-
-volScalarField p_rgh
-(
-    IOobject
-    (
-        "p_rgh",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-// Force p_rgh to be consistent with p
-p_rgh = p - rho*gh;
-
-mesh.setFluxRequired(p_rgh.name());
-
-multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
-
-forAll(Y, i)
-{
-    fields.add(Y[i]);
-}
-fields.add(thermo.he());
-
-IOdictionary additionalControlsDict
-(
-    IOobject
-    (
-        "additionalControls",
-        runTime.constant(),
-        mesh,
-        IOobject::MUST_READ_IF_MODIFIED,
-        IOobject::NO_WRITE
-    )
-);
-
-Switch solvePrimaryRegion
-(
-    additionalControlsDict.lookup("solvePrimaryRegion")
-);
diff --git a/applications/solvers/combustion/fireFoam/fireFoam.C b/applications/solvers/combustion/fireFoam/fireFoam.C
index fd8e2a435..3e1689f65 100644
--- a/applications/solvers/combustion/fireFoam/fireFoam.C
+++ b/applications/solvers/combustion/fireFoam/fireFoam.C
@@ -54,7 +54,6 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
-    #include "createMRF.H"
     #include "createFvOptions.H"
     #include "createClouds.H"
     #include "createSurfaceFilmModel.H"
diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
index 7b5249d57..0d42787ff 100644
--- a/applications/solvers/combustion/fireFoam/pEqn.H
+++ b/applications/solvers/combustion/fireFoam/pEqn.H
@@ -4,7 +4,7 @@ volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
-surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig("phig", -rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
 surfaceScalarField phiHbyA
 (
@@ -25,9 +25,10 @@ while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix p_rghEqn
     (
-        fvc::ddt(psi, rho)*gh
+        fvm::ddt(psi, p_rgh)
+      + fvc::ddt(psi, rho)*gh
+      + fvc::ddt(psi)*pRef
       + fvc::div(phiHbyA)
-      + fvm::ddt(psi, p_rgh)
       - fvm::laplacian(rhorAUf, p_rgh)
      ==
         parcels.Srho()
@@ -46,7 +47,7 @@ while (pimple.correctNonOrthogonal())
     }
 }
 
-p = p_rgh + rho*gh;
+p = p_rgh + rho*gh + pRef;
 
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
diff --git a/applications/solvers/combustion/fireFoam/phrghEqn.H b/applications/solvers/combustion/fireFoam/phrghEqn.H
new file mode 100644
index 000000000..3c5175f41
--- /dev/null
+++ b/applications/solvers/combustion/fireFoam/phrghEqn.H
@@ -0,0 +1,62 @@
+if (pimple.dict().lookupOrDefault<bool>("hydrostaticInitialization", false))
+{
+    volScalarField& ph_rgh = regIOobject::store
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "ph_rgh",
+                "0",
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh
+        )
+    );
+
+    if (equal(runTime.value(), 0))
+    {
+        p = ph_rgh + rho*gh + pRef;
+        thermo.correct();
+        rho = thermo.rho();
+
+        label nCorr
+        (
+            pimple.dict().lookupOrDefault<label>("nHydrostaticCorrectors", 5)
+        );
+
+        for (label i=0; i<nCorr; i++)
+        {
+            surfaceScalarField rhof("rhof", fvc::interpolate(rho));
+
+            surfaceScalarField phig
+            (
+                "phig",
+               -rhof*ghf*fvc::snGrad(rho)*mesh.magSf()
+            );
+
+            // Update the pressure BCs to ensure flux consistency
+            constrainPressure(ph_rgh, rho, U, phig, rhof);
+
+            fvScalarMatrix ph_rghEqn
+            (
+                fvm::laplacian(rhof, ph_rgh) == fvc::div(phig)
+            );
+
+            ph_rghEqn.solve();
+
+            p = ph_rgh + rho*gh + pRef;
+            thermo.correct();
+            rho = thermo.rho();
+
+            Info<< "Hydrostatic pressure variation "
+                << (max(ph_rgh) - min(ph_rgh)).value() << endl;
+        }
+
+        ph_rgh.write();
+
+        p_rgh = ph_rgh;
+    }
+}
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 87996b789..987729b1a 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -202,6 +202,7 @@ $(derivedFvPatchFields)/waveSurfacePressure/waveSurfacePressureFvPatchScalarFiel
 $(derivedFvPatchFields)/interstitialInletVelocity/interstitialInletVelocityFvPatchVectorField.C
 $(derivedFvPatchFields)/prghPressure/prghPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/prghTotalPressure/prghTotalPressureFvPatchScalarField.C
+$(derivedFvPatchFields)/prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/fixedProfile/fixedProfileFvPatchFields.C
 
 fvsPatchFields = fields/fvsPatchFields
diff --git a/src/finiteVolume/cfdTools/general/include/readpRef.H b/src/finiteVolume/cfdTools/general/include/readpRef.H
new file mode 100644
index 000000000..78ae5886e
--- /dev/null
+++ b/src/finiteVolume/cfdTools/general/include/readpRef.H
@@ -0,0 +1,13 @@
+    Info<< "\nReading pRef" << endl;
+    uniformDimensionedScalarField pRef
+    (
+        IOobject
+        (
+            "pRef",
+            runTime.constant(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        dimensionedScalar("pRef", dimPressure, 0)
+    );
diff --git a/src/finiteVolume/fields/fvPatchFields/basic/fixedValue/fixedValueFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/basic/fixedValue/fixedValueFvPatchField.H
index aef9052b5..fdf3ea51e 100644
--- a/src/finiteVolume/fields/fvPatchFields/basic/fixedValue/fixedValueFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/basic/fixedValue/fixedValueFvPatchField.H
@@ -43,7 +43,7 @@ Description
     myPatch
     {
         type            fixedValue;
-        value           uniform 0;  // example for scalar field usage
+        value           uniform 0;  // Example for scalar field usage
     }
     \endverbatim
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C
new file mode 100644
index 000000000..f0eb72c62
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.C
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 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 "prghTotalHydrostaticPressureFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "uniformDimensionedFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::prghTotalHydrostaticPressureFvPatchScalarField::
+prghTotalHydrostaticPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    UName_("U"),
+    phiName_("phi"),
+    rhoName_("rho"),
+    ph_rghName_("ph_rgh")
+{}
+
+
+Foam::prghTotalHydrostaticPressureFvPatchScalarField::
+prghTotalHydrostaticPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict),
+    UName_(dict.lookupOrDefault<word>("U", "U")),
+    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    ph_rghName_(dict.lookupOrDefault<word>("ph_rgh", "ph_rgh"))
+{}
+
+
+Foam::prghTotalHydrostaticPressureFvPatchScalarField::
+prghTotalHydrostaticPressureFvPatchScalarField
+(
+    const prghTotalHydrostaticPressureFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    UName_(ptf.UName_),
+    phiName_(ptf.phiName_),
+    rhoName_(ptf.rhoName_),
+    ph_rghName_(ptf.ph_rghName_)
+{}
+
+
+Foam::prghTotalHydrostaticPressureFvPatchScalarField::
+prghTotalHydrostaticPressureFvPatchScalarField
+(
+    const prghTotalHydrostaticPressureFvPatchScalarField& ptf
+)
+:
+    fixedValueFvPatchScalarField(ptf),
+    UName_(ptf.UName_),
+    phiName_(ptf.phiName_),
+    rhoName_(ptf.rhoName_),
+    ph_rghName_(ptf.ph_rghName_)
+{}
+
+
+Foam::prghTotalHydrostaticPressureFvPatchScalarField::
+prghTotalHydrostaticPressureFvPatchScalarField
+(
+    const prghTotalHydrostaticPressureFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(ptf, iF),
+    UName_(ptf.UName_),
+    phiName_(ptf.phiName_),
+    rhoName_(ptf.rhoName_),
+    ph_rghName_(ptf.ph_rghName_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::prghTotalHydrostaticPressureFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const scalarField& rhop =
+        patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+
+    const scalarField& ph_rghp =
+        patch().lookupPatchField<volScalarField, scalar>(ph_rghName_);
+
+    const scalarField& phip =
+        patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);
+
+    const vectorField& Up =
+        patch().lookupPatchField<volVectorField, vector>(UName_);
+
+    operator==
+    (
+        ph_rghp
+      - 0.5*rhop*(1.0 - pos(phip))*magSqr(Up)
+    );
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::prghTotalHydrostaticPressureFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
+    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    writeEntryIfDifferent<word>(os, "ph_rghName", "ph_rghName", ph_rghName_);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        prghTotalHydrostaticPressureFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.H
new file mode 100644
index 000000000..a37e645f1
--- /dev/null
+++ b/src/finiteVolume/fields/fvPatchFields/derived/prghTotalHydrostaticPressure/prghTotalHydrostaticPressureFvPatchScalarField.H
@@ -0,0 +1,202 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 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::prghTotalHydrostaticPressureFvPatchScalarField
+
+Group
+    grpGenericBoundaryConditions
+
+Description
+    This boundary condition provides static pressure condition for p_rgh,
+    calculated as:
+
+        \f[
+            p_rgh = ph_rgh - 0.5 \rho |U|^2
+        \f]
+
+    where
+    \vartable
+        p_rgh   | Pressure - \rho g.(h - hRef) [Pa]
+        ph_rgh  | Hydrostatic pressure - \rho g.(h - hRef) [Pa]
+        h       | Height in the opposite direction to gravity
+        hRef    | Reference height in the opposite direction to gravity
+        \rho    | Density
+        g       | Acceleration due to gravity [m/s^2]
+    \endtable
+
+    \heading Patch usage
+
+    \table
+        Property     | Description             | Required    | Default value
+        U            | Velocity field name     | no          | U
+        phi          | Flux field name         | no          | phi
+        rho          | Density field name      | no          | rho
+        ph_rgh       | ph_rgh field name       | no          | ph_rgh
+        value        | Patch face values       | yes         |
+    \endtable
+
+    Example of the boundary condition specification:
+    \verbatim
+    myPatch
+    {
+        type            prghTotalHydrostaticPressure;
+        value           uniform 0;
+    }
+    \endverbatim
+
+SeeAlso
+    Foam::fixedValueFvPatchScalarField
+    Foam::prghTotalPressureFvPatchScalarField
+
+SourceFiles
+    prghTotalHydrostaticPressureFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef prghTotalHydrostaticPressureFvPatchScalarField_H
+#define prghTotalHydrostaticPressureFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+         Class prghTotalHydrostaticPressureFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class prghTotalHydrostaticPressureFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+
+protected:
+
+    // Protected data
+
+        //- Name of the velocity field
+        word UName_;
+
+        //- Name of the flux transporting the field
+        word phiName_;
+
+        //- Name of density field
+        word rhoName_;
+
+        //- Name of hydrostatic pressure field
+        word ph_rghName_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("prghTotalHydrostaticPressure");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        prghTotalHydrostaticPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        prghTotalHydrostaticPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  prghTotalHydrostaticPressureFvPatchScalarField onto a new patch
+        prghTotalHydrostaticPressureFvPatchScalarField
+        (
+            const prghTotalHydrostaticPressureFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        prghTotalHydrostaticPressureFvPatchScalarField
+        (
+            const prghTotalHydrostaticPressureFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField >
+            (
+                new prghTotalHydrostaticPressureFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        prghTotalHydrostaticPressureFvPatchScalarField
+        (
+            const prghTotalHydrostaticPressureFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new prghTotalHydrostaticPressureFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/prghTotalPressure/prghTotalPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/prghTotalPressure/prghTotalPressureFvPatchScalarField.H
index 1723e2091..120133223 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/prghTotalPressure/prghTotalPressureFvPatchScalarField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/prghTotalPressure/prghTotalPressureFvPatchScalarField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,7 +32,7 @@ Description
     calculated as:
 
         \f[
-            p_rgh = p - \rho g (h - hRef)
+            p_rgh = p - \rho g.(h - hRef)
             p = p0 - 0.5 \rho |U|^2
         \f]
 
@@ -43,18 +43,18 @@ Description
         p0      | Total pressure [Pa]
         h       | Height in the opposite direction to gravity
         hRef    | Reference height in the opposite direction to gravity
-        \rho    | density
-        g       | acceleration due to gravity [m/s^2]
+        \rho    | Density
+        g       | Acceleration due to gravity [m/s^2]
     \endtable
 
     \heading Patch usage
 
     \table
         Property     | Description             | Required    | Default value
-        U            | velocity field name     | no          | U
-        phi          | flux field name         | no          | phi
-        rho          | density field name      | no          | rho
-        p0           | total pressure          | yes         |
+        U            | Velocity field name     | no          | U
+        phi          | Flux field name         | no          | phi
+        rho          | Density field name      | no          | rho
+        p0           | Total pressure          | yes         |
     \endtable
 
     Example of the boundary condition specification:
@@ -62,11 +62,7 @@ Description
     myPatch
     {
         type            prghTotalPressure;
-        U               U;
-        phi             phi;
-        rho             rho;
         p0              uniform 0;
-        value           uniform 0; // optional initial value
     }
     \endverbatim
 
diff --git a/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/0/p_rgh b/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/0/p_rgh
index d95b85f7b..dbef4a438 100644
--- a/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/0/p_rgh
+++ b/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/0/p_rgh
@@ -17,32 +17,20 @@ FoamFile
 
 dimensions      [1 -1 -2 0 0 0 0];
 
-internalField   uniform 101325;
+internalField   uniform 0;
 
 boundaryField
 {
-    outlet
+    "(side|outlet)"
     {
-        type            fixedFluxPressure;
+        type            prghTotalHydrostaticPressure;
+        p0              $internalField;
         value           $internalField;
     }
 
     ground
     {
         type            fixedFluxPressure;
-
-    }
-
-    side
-    {
-        type            totalPressure;
-        U               U;
-        phi             phi;
-        rho             rho;
-        psi             none;
-        gamma           1.4;
-        p0              $internalField;
-        value           $internalField;
     }
 
     inlet
diff --git a/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/Allclean b/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/Allclean
index 8f51891b2..615ace08b 100755
--- a/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/Allclean
+++ b/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/Allclean
@@ -6,9 +6,9 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 cleanCase
 
+rm -f 0/ph_rgh
 rm -rf constant/filmRegion/polyMesh
 rm -rf constant/pyrolysisRegion/polyMesh
-
 rm -rf system/pyrolysisRegion/filmRegion
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/Allrun b/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/Allrun
index 6101e9637..4b568eed0 100755
--- a/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/Allrun
+++ b/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/Allrun
@@ -33,9 +33,11 @@ rm log.createPatch
 # create actual patches
 runApplication createPatch -region filmRegion -overwrite
 
+cp 0/ph_rgh.orig 0/ph_rgh
+
 # Run
 runApplication `getApplication`
 
 paraFoam -touchAll
 
-#------------------------------------------------------------------------------
\ No newline at end of file
+#------------------------------------------------------------------------------
diff --git a/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/constant/pRef b/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/constant/pRef
new file mode 100644
index 000000000..2f1d6789f
--- /dev/null
+++ b/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/constant/pRef
@@ -0,0 +1,21 @@
+/*--------------------------------*- 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       uniformDimensionedScalarField;
+    location    "constant";
+    object      pRef;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+value           101325;
+
+// ************************************************************************* //
diff --git a/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/system/fvSolution b/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/system/fvSolution
index 070876d36..c41caf799 100644
--- a/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/system/fvSolution
+++ b/tutorials/combustion/fireFoam/les/flameSpreadWaterSuppressionPanel/system/fvSolution
@@ -58,6 +58,11 @@ solvers
         mergeLevels      1;
     };
 
+    ph_rgh
+    {
+        $p_rgh;
+    }
+
     "(U|Yi|h|k)"
     {
         solver          PBiCG;
@@ -93,6 +98,9 @@ PIMPLE
     nOuterCorrectors  1;
     nCorrectors       2;
     nNonOrthogonalCorrectors 0;
+
+    hydrostaticInitialization yes;
+    nHydrostaticCorrectors 5;
 }
 
 relaxationFactors
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/p_rgh b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/p_rgh
index 8beb74136..a8142fbd7 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/p_rgh
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/0/p_rgh
@@ -17,7 +17,7 @@ FoamFile
 
 dimensions      [1 -1 -2 0 0 0 0];
 
-internalField   uniform 101325;
+internalField   uniform 0;
 
 boundaryField
 {
@@ -35,12 +35,7 @@ boundaryField
 
     "(top|sides)"
     {
-        type            totalPressure;
-        U               U;
-        phi             phi;
-        rho             rho;
-        psi             none;
-        gamma           1.4;
+        type            prghTotalHydrostaticPressure;
         p0              $internalField;
         value           $internalField;
     }
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/Allclean b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/Allclean
index df17cd588..ed4bd14d2 100755
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/Allclean
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/Allclean
@@ -6,6 +6,7 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 cleanCase
 
+rm -f 0/ph_rgh
 rm -rf constant/panelRegion/polyMesh
 rm  -f 0/polyMesh/cellMap
 rm  -f constant/polyMesh/boundary
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/Allrun b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/Allrun
index befcb5a01..2ee23ef99 100755
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/Allrun
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/Allrun
@@ -26,6 +26,8 @@ runApplication -s fBurner \
 
 runApplication createPatch -overwrite
 
+cp 0/ph_rgh.orig 0/ph_rgh
+
 runApplication -s master \
     decomposePar -force
 
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/pRef b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/pRef
new file mode 100644
index 000000000..2f1d6789f
--- /dev/null
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/pRef
@@ -0,0 +1,21 @@
+/*--------------------------------*- 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       uniformDimensionedScalarField;
+    location    "constant";
+    object      pRef;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+value           101325;
+
+// ************************************************************************* //
diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSolution b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSolution
index 5691a3dbd..0182eed8b 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSolution
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSolution
@@ -43,6 +43,12 @@ solvers
         relTol              0;
     };
 
+    ph_rgh
+    {
+        $p_rgh;
+    }
+
+
     "(U|Yi|k|h|omega)"
     {
         solver          PBiCG;
@@ -88,6 +94,9 @@ PIMPLE
     nOuterCorrectors  1;
     nCorrectors       2;
     nNonOrthogonalCorrectors 0;
+
+    hydrostaticInitialization yes;
+    nHydrostaticCorrectors 5;
 }
 
 relaxationFactors
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p_rgh b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p_rgh
index 5868c5276..eb19cbc78 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p_rgh
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p_rgh
@@ -17,19 +17,14 @@ FoamFile
 
 dimensions      [1 -1 -2 0 0 0 0];
 
-internalField   uniform 101325;
+internalField   uniform 0;
 
 boundaryField
 {
     "(outlet|sides)"
     {
-        type            totalPressure;
+        type            prghTotalHydrostaticPressure;
         p0              $internalField;
-        U               U;
-        phi             phi;
-        rho             none;
-        psi             none;
-        gamma           1;
         value           $internalField;
     }
 
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/Allclean b/tutorials/combustion/fireFoam/les/smallPoolFire2D/Allclean
new file mode 100755
index 000000000..333e5d8d2
--- /dev/null
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/Allclean
@@ -0,0 +1,11 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+
+rm -f 0/ph_rgh
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/Allrun b/tutorials/combustion/fireFoam/les/smallPoolFire2D/Allrun
index 11ae4b0d7..5bfda3a6b 100755
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/Allrun
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/Allrun
@@ -9,6 +9,8 @@ runApplication blockMesh
 runApplication topoSet
 runApplication createPatch -overwrite
 
+cp 0/ph_rgh.orig 0/ph_rgh
+
 # Run
 runApplication `getApplication`
 
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/pRef b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/pRef
new file mode 100644
index 000000000..2f1d6789f
--- /dev/null
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/pRef
@@ -0,0 +1,21 @@
+/*--------------------------------*- 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       uniformDimensionedScalarField;
+    location    "constant";
+    object      pRef;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+value           101325;
+
+// ************************************************************************* //
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSolution b/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSolution
index a7e9319e2..674390637 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSolution
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSolution
@@ -52,6 +52,11 @@ solvers
         relTol              0;
     };
 
+    ph_rgh
+    {
+        $p_rgh;
+    }
+
 
     "(U|Yi|k|h|omega)"
     {
@@ -101,6 +106,9 @@ PIMPLE
     nOuterCorrectors  1;
     nCorrectors       2;
     nNonOrthogonalCorrectors 0;
+
+    hydrostaticInitialization yes;
+    nHydrostaticCorrectors 5;
 }
 
 relaxationFactors
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/p_rgh b/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/p_rgh
index 15a149ba6..746686864 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/p_rgh
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/0/p_rgh
@@ -17,18 +17,14 @@ FoamFile
 
 dimensions      [1 -1 -2 0 0 0 0];
 
-internalField   uniform 101325;
+internalField   uniform 0;
 
 boundaryField
 {
     "(outlet|sides)"
     {
-        type            totalPressure;
+        type            prghTotalHydrostaticPressure;
         p0              $internalField;
-        U               U;
-        phi             phi;
-        rho             rho;
-        psi             none;
         gamma           1;
         value           $internalField;
     }
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/Allclean b/tutorials/combustion/fireFoam/les/smallPoolFire3D/Allclean
new file mode 100755
index 000000000..333e5d8d2
--- /dev/null
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/Allclean
@@ -0,0 +1,11 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+
+rm -f 0/ph_rgh
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/Allrun b/tutorials/combustion/fireFoam/les/smallPoolFire3D/Allrun
index 81460ea38..d05756391 100755
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/Allrun
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/Allrun
@@ -9,6 +9,9 @@ application=`getApplication`
 runApplication blockMesh
 runApplication topoSet
 runApplication createPatch -overwrite
+
+cp 0/ph_rgh.orig 0/ph_rgh
+
 runApplication decomposePar -force
 
 # Run
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/pRef b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/pRef
new file mode 100644
index 000000000..2f1d6789f
--- /dev/null
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/pRef
@@ -0,0 +1,21 @@
+/*--------------------------------*- 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       uniformDimensionedScalarField;
+    location    "constant";
+    object      pRef;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+value           101325;
+
+// ************************************************************************* //
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/fvSolution b/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/fvSolution
index 574c3cfd7..a80537572 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/fvSolution
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/system/fvSolution
@@ -44,6 +44,11 @@ solvers
         relTol              0;
     };
 
+    ph_rgh
+    {
+        $p_rgh;
+    }
+
 
     "(U|Yi|k|h)"
     {
@@ -92,6 +97,9 @@ PIMPLE
     nOuterCorrectors  1;
     nCorrectors       2;
     nNonOrthogonalCorrectors 0;
+
+    hydrostaticInitialization yes;
+    nHydrostaticCorrectors 5;
 }
 
 relaxationFactors