diff --git a/applications/solvers/combustion/reactingFoam/buoyantReactingFoam/buoyantReactingFoam.C b/applications/solvers/combustion/reactingFoam/buoyantReactingFoam/buoyantReactingFoam.C
index 4fbda271ef..9c95e2cdf6 100644
--- a/applications/solvers/combustion/reactingFoam/buoyantReactingFoam/buoyantReactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/buoyantReactingFoam/buoyantReactingFoam.C
@@ -43,6 +43,7 @@ Description
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
 #include "pressureReference.H"
+#include "hydrostaticInitialisation.H"
 #include "CorrectPhi.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
diff --git a/applications/solvers/combustion/reactingFoam/buoyantReactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/buoyantReactingFoam/createFields.H
index 870f98bfa9..1ddaff17c5 100644
--- a/applications/solvers/combustion/reactingFoam/buoyantReactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/buoyantReactingFoam/createFields.H
@@ -96,7 +96,18 @@ pressureReference pressureReference
 
 mesh.setFluxRequired(p_rgh.name());
 
-#include "hydrostaticInitialisation.H"
+hydrostaticInitialisation
+(
+    p_rgh,
+    p,
+    rho,
+    U,
+    gh,
+    ghf,
+    pRef,
+    thermo,
+    pimple.dict()
+);
 
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index 4dbee77899..642f3cd898 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -41,6 +41,7 @@ Description
 #include "fluidThermophysicalTransportModel.H"
 #include "pimpleControl.H"
 #include "pressureReference.H"
+#include "hydrostaticInitialisation.H"
 #include "CorrectPhi.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
index 29b0d30a4f..8f04e6367c 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
@@ -76,9 +76,6 @@ volScalarField p_rgh
     mesh
 );
 
-// Force p_rgh to be consistent with p
-p_rgh = p - rho*gh - pRef;
-
 pressureReference pressureReference
 (
     p,
@@ -89,6 +86,19 @@ pressureReference pressureReference
 
 mesh.setFluxRequired(p_rgh.name());
 
+hydrostaticInitialisation
+(
+    p_rgh,
+    p,
+    rho,
+    U,
+    gh,
+    ghf,
+    pRef,
+    thermo,
+    pimple.dict()
+);
+
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt
 (
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
index a5d7f26838..12b0ae5a92 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
@@ -35,6 +35,7 @@ Description
 #include "fluidThermophysicalTransportModel.H"
 #include "simpleControl.H"
 #include "pressureReference.H"
+#include "hydrostaticInitialisation.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
 
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
index fe3e113e03..893b50fcdc 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
@@ -74,9 +74,6 @@ volScalarField p_rgh
     mesh
 );
 
-// Force p_rgh to be consistent with p
-p_rgh = p - rho*gh - pRef;
-
 pressureReference pressureReference
 (
     p,
@@ -87,6 +84,19 @@ pressureReference pressureReference
 
 mesh.setFluxRequired(p_rgh.name());
 
+hydrostaticInitialisation
+(
+    p_rgh,
+    p,
+    rho,
+    U,
+    gh,
+    ghf,
+    pRef,
+    thermo,
+    simple.dict()
+);
+
 dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 
 #include "createMRF.H"
diff --git a/src/finiteVolume/cfdTools/compressible/hydrostaticInitialisation.H b/src/finiteVolume/cfdTools/compressible/hydrostaticInitialisation.H
deleted file mode 100644
index 6ea10a0568..0000000000
--- a/src/finiteVolume/cfdTools/compressible/hydrostaticInitialisation.H
+++ /dev/null
@@ -1,67 +0,0 @@
-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;
-    }
-}
-else
-{
-    // Force p_rgh to be consistent with p
-    p_rgh = p - rho*gh - pRef;
-}
diff --git a/src/thermophysicalModels/basic/fluidThermo/hydrostaticInitialisation.C b/src/thermophysicalModels/basic/fluidThermo/hydrostaticInitialisation.C
new file mode 100644
index 0000000000..5f1857245b
--- /dev/null
+++ b/src/thermophysicalModels/basic/fluidThermo/hydrostaticInitialisation.C
@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 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 "hydrostaticInitialisation.H"
+
+#include "fluidThermo.H"
+#include "fvmLaplacian.H"
+#include "fvcDiv.H"
+#include "fvcSnGrad.H"
+#include "surfaceInterpolate.H"
+#include "constrainPressure.H"
+#include "uniformDimensionedFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::hydrostaticInitialisation
+(
+    volScalarField& p_rgh,
+    volScalarField& p,
+    volScalarField& rho,
+    const volVectorField& U,
+    const volScalarField& gh,
+    const surfaceScalarField& ghf,
+    const uniformDimensionedScalarField& pRef,
+    fluidThermo& thermo,
+    const dictionary& dict
+)
+{
+    if (dict.lookupOrDefault<bool>("hydrostaticInitialisation", false))
+    {
+        const fvMesh& mesh = p_rgh.mesh();
+
+        volScalarField& ph_rgh = regIOobject::store
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "ph_rgh",
+                    "0",
+                    mesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh
+            )
+        );
+
+        if (!mesh.time().restart())
+        {
+            p = ph_rgh + rho*gh + pRef;
+            thermo.correct();
+            rho = thermo.rho();
+
+            label nCorr
+            (
+                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;
+        }
+    }
+    else
+    {
+        // Force p_rgh to be consistent with p
+        p_rgh = p - rho*gh - pRef;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/basic/fluidThermo/hydrostaticInitialisation.H b/src/thermophysicalModels/basic/fluidThermo/hydrostaticInitialisation.H
new file mode 100644
index 0000000000..ea88a7601c
--- /dev/null
+++ b/src/thermophysicalModels/basic/fluidThermo/hydrostaticInitialisation.H
@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 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::hydrostaticInitialisation
+
+Description
+    Optional hydrostatic initialisation of p_rgh and p by solving for and
+    caching the hydrostatic ph_rgh and updating the density such that
+
+        p = ph_rgh + rho*gh + pRef
+
+    This initialisation process is applied at the beginning of the run (not on
+    restart) if the \c hydrostaticInitialisation switch is set true in
+    fvSolution/PIMPLE or fvSolution/SIMPLE.  The calculation is iterative if the
+    density is a function of pressure and an optional number of iterations \c
+    nHydrostaticCorrectors may be specified which defaults to 5.
+
+SourceFiles
+    hydrostaticInitialisation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef hydrostaticInitialisation_H
+#define hydrostaticInitialisation_H
+
+#include "volFieldsFwd.H"
+#include "surfaceFieldsFwd.H"
+#include "uniformDimensionedFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    class fluidThermo;
+
+    void hydrostaticInitialisation
+    (
+        volScalarField& p_rgh,
+        volScalarField& p,
+        volScalarField& rho,
+        const volVectorField& U,
+        const volScalarField& gh,
+        const surfaceScalarField& ghf,
+        const uniformDimensionedScalarField& pRef,
+        fluidThermo& thermo,
+        const dictionary& dict
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/combustion/buoyantReactingFoam/RAS/smallPoolFire2D/system/fvSolution b/tutorials/combustion/buoyantReactingFoam/RAS/smallPoolFire2D/system/fvSolution
index c959476d57..2f5940f2e9 100644
--- a/tutorials/combustion/buoyantReactingFoam/RAS/smallPoolFire2D/system/fvSolution
+++ b/tutorials/combustion/buoyantReactingFoam/RAS/smallPoolFire2D/system/fvSolution
@@ -98,7 +98,7 @@ PIMPLE
     nCorrectors       2;
     nNonOrthogonalCorrectors 0;
 
-    hydrostaticInitialization yes;
+    hydrostaticInitialisation yes;
     nHydrostaticCorrectors 5;
 }
 
diff --git a/tutorials/combustion/buoyantReactingFoam/RAS/smallPoolFire3D/system/fvSolution b/tutorials/combustion/buoyantReactingFoam/RAS/smallPoolFire3D/system/fvSolution
index 7c884e7b02..b9ed9c37e8 100644
--- a/tutorials/combustion/buoyantReactingFoam/RAS/smallPoolFire3D/system/fvSolution
+++ b/tutorials/combustion/buoyantReactingFoam/RAS/smallPoolFire3D/system/fvSolution
@@ -89,7 +89,7 @@ PIMPLE
     nCorrectors       2;
     nNonOrthogonalCorrectors 0;
 
-    hydrostaticInitialization yes;
+    hydrostaticInitialisation yes;
     nHydrostaticCorrectors 5;
 }