diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
index c29d451091..865af56184 100644
--- a/applications/solvers/DNS/dnsFoam/dnsFoam.C
+++ b/applications/solvers/DNS/dnsFoam/dnsFoam.C
@@ -86,7 +86,7 @@ int main(int argc, char *argv[])
         for (int corr=1; corr<=1; corr++)
         {
             volScalarField rAU(1.0/UEqn.A());
-            surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+            surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
             volVectorField HbyA("HbyA", U);
             HbyA = rAU*UEqn.H();
 
@@ -94,12 +94,12 @@ int main(int argc, char *argv[])
             (
                 "phiHbyA",
                 (fvc::interpolate(HbyA) & mesh.Sf())
-              + Dp*fvc::ddtCorr(U, phi)
+              + rAUf*fvc::ddtCorr(U, phi)
             );
 
             fvScalarMatrix pEqn
             (
-                fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
+                fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
             );
 
             pEqn.solve();
diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
index 83a44ab7c2..995a59a737 100644
--- a/applications/solvers/combustion/XiFoam/pEqn.H
+++ b/applications/solvers/combustion/XiFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
@@ -14,7 +14,7 @@ if (pimple.transonic())
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )/fvc::interpolate(rho)
     );
 
@@ -26,7 +26,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
@@ -48,7 +48,7 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
 
@@ -60,7 +60,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index a92e7285ab..384ecbee59 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/applications/solvers/combustion/engineFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
@@ -15,7 +15,7 @@ if (pimple.transonic())
        *(
             (
                 (fvc::interpolate(rho*HbyA) & mesh.Sf())
-              //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+              //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
             )/fvc::interpolate(rho)
           - fvc::meshPhi(rho, U)
         )
@@ -51,7 +51,7 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+          //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
       - fvc::interpolate(rho)*fvc::meshPhi(rho, U)
     );
diff --git a/applications/solvers/combustion/fireFoam/fireFoam.C b/applications/solvers/combustion/fireFoam/fireFoam.C
index 14e46bef66..b8c48549a6 100644
--- a/applications/solvers/combustion/fireFoam/fireFoam.C
+++ b/applications/solvers/combustion/fireFoam/fireFoam.C
@@ -41,6 +41,7 @@ Description
 #include "psiCombustionModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
index c98ac61f3b..6652eab564 100644
--- a/applications/solvers/combustion/fireFoam/pEqn.H
+++ b/applications/solvers/combustion/fireFoam/pEqn.H
@@ -1,27 +1,35 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
-phi.boundaryField() =
-    fvc::interpolate(rho.boundaryField()*U.boundaryField())
-  & mesh.Sf().boundaryField();
 
-surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (
         (fvc::interpolate(rho*HbyA) & mesh.Sf())
-      + Dp*fvc::ddtCorr(rho, U, phi)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
     )
   + phig
 );
 
-fvOptions.makeRelative(phiHbyA);
+fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the fixedFluxPressure BCs to ensure flux consistency
+setSnGrad<fixedFluxPressureFvPatchScalarField>
+(
+    p_rgh.boundaryField(),
+    (
+        phiHbyA.boundaryField()
+      - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+       *rho.boundaryField()
+    )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+);
 
 while (pimple.correctNonOrthogonal())
 {
@@ -30,7 +38,7 @@ while (pimple.correctNonOrthogonal())
         fvc::ddt(psi, rho)*gh
       + fvc::div(phiHbyA)
       + fvm::ddt(psi, p_rgh)
-      - fvm::laplacian(Dp, p_rgh)
+      - fvm::laplacian(rhorAUf, p_rgh)
      ==
         parcels.Srho()
       + surfaceFilm.Srho()
@@ -44,7 +52,7 @@ while (pimple.correctNonOrthogonal())
     if (pimple.finalNonOrthogonalIter())
     {
         phi = phiHbyA + p_rghEqn.flux();
-        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/Dp);
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
         U.correctBoundaryConditions();
         fvOptions.correct(U);
     }
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index 0dc5d422db..f7a3004ac2 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
@@ -14,7 +14,7 @@ if (pimple.transonic())
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )/fvc::interpolate(rho)
     );
 
@@ -48,7 +48,7 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
 
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
index 8602725bdb..8894c7a9d7 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
@@ -6,25 +6,36 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
       + phig
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+           *rho.boundaryField()
+        )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+    );
+
     fvScalarMatrix p_rghDDtEqn
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
@@ -38,7 +49,7 @@
         fvScalarMatrix p_rghEqn
         (
             p_rghDDtEqn
-          - fvm::laplacian(Dp, p_rgh)
+          - fvm::laplacian(rhorAUf, p_rgh)
         );
 
         fvOptions.constrain(p_rghEqn);
@@ -55,7 +66,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
             K = 0.5*magSqr(U);
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C
index 19e9f96f51..fd4b6f58eb 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C
@@ -36,6 +36,7 @@ Description
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
index 3c92ce8119..f9b1bd65e4 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
@@ -6,7 +6,7 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -18,7 +18,7 @@
             "phiHbyA",
             (
                 (fvc::interpolate(rho*HbyA) & mesh.Sf())
-              + Dp*fvc::ddtCorr(rho, U, phi)
+              + rhorAUf*fvc::ddtCorr(rho, U, phi)
             )/fvc::interpolate(rho)
         );
 
@@ -39,7 +39,7 @@
             fvScalarMatrix pEqn
             (
                 pDDtEqn
-              - fvm::laplacian(Dp, p)
+              - fvm::laplacian(rhorAUf, p)
              ==
                 fvOptions(psi, p, rho.name())
             );
@@ -61,7 +61,7 @@
             "phiHbyA",
             (
                 (fvc::interpolate(rho*HbyA) & mesh.Sf())
-              + Dp*fvc::ddtCorr(rho, U, phi)
+              + rhorAUf*fvc::ddtCorr(rho, U, phi)
             )
         );
 
@@ -80,7 +80,7 @@
             fvScalarMatrix pEqn
             (
                 pDDtEqn
-              - fvm::laplacian(Dp, p)
+              - fvm::laplacian(rhorAUf, p)
             );
 
             fvOptions.constrain(pEqn);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 373bafcd3b..a4e5fa20da 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -4,7 +4,7 @@ rho = min(rho, rhoMax);
 rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
@@ -22,7 +22,7 @@ if (pimple.transonic())
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )/fvc::interpolate(rho)
     );
 
@@ -34,7 +34,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
           ==
             fvOptions(psi, p, rho.name())
         );
@@ -56,14 +56,12 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
-    volScalarField Dp("Dp", rho*rAU);
-
     while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
@@ -71,7 +69,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
           ==
             fvOptions(psi, p, rho.name())
         );
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index d2681663e6..1d70a0f354 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -4,7 +4,7 @@ rho = min(rho, rhoMax);
 rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
@@ -22,7 +22,7 @@ if (pimple.transonic())
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phiAbs)
+          + rhorAUf*fvc::ddtCorr(rho, U, phiAbs)
         )/fvc::interpolate(rho)
     );
 
@@ -34,7 +34,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
           ==
             fvOptions(psi, p, rho.name())
         );
@@ -55,7 +55,7 @@ else
     (
         "phiHbyA",
         (fvc::interpolate(rho*HbyA) & mesh.Sf())
-      + Dp*fvc::ddtCorr(rho, U, phiAbs)
+      + rhorAUf*fvc::ddtCorr(rho, U, phiAbs)
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
@@ -67,7 +67,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
           ==
             fvOptions(psi, p, rho.name())
         );
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
index f8837f0aa4..0503db4207 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimplecFoam/pEqn.H
@@ -35,7 +35,7 @@ if (pimple.transonic())
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField Dp("Dp", rho*rAtU);
+    volScalarField rhorAtU("rhorAtU", rho*rAtU);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -44,7 +44,7 @@ if (pimple.transonic())
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
           + fvc::div(phic)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAtU, p)
          ==
             fvOptions(psi, p, rho.name())
         );
@@ -75,7 +75,7 @@ else
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField Dp("Dp", rho*rAtU);
+    volScalarField rhorAtU("rhorAtU", rho*rAtU);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -83,7 +83,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAtU, p)
          ==
             fvOptions(psi, p, rho.name())
         );
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
index d61a80ffa8..c656b3c3d7 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H
@@ -26,7 +26,7 @@ if (simple.transonic())
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField Dp("Dp", rho*rAtU);
+    volScalarField rhorAtU("rhorAtU", rho*rAtU);
 
     while (simple.correctNonOrthogonal())
     {
@@ -34,7 +34,7 @@ if (simple.transonic())
         (
             fvm::div(phid, p)
           + fvc::div(phic)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAtU, p)
          ==
             fvOptions(psi, p, rho.name())
         );
@@ -67,14 +67,14 @@ else
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField Dp("Dp", rho*rAtU);
+    volScalarField rhorAtU("rhorAtU", rho*rAtU);
 
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAtU, p)
           ==
             fvOptions(psi, p, rho.name())
         );
diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H
index 2832bf421a..edde4b198e 100644
--- a/applications/solvers/compressible/sonicFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
@@ -12,7 +12,7 @@ surfaceScalarField phid
     fvc::interpolate(psi)*
     (
         (mesh.Sf() & fvc::interpolate(rho*HbyA))
-      + Dp*fvc::ddtCorr(rho, U, phi)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
     )/fvc::interpolate(rho)
 );
 
@@ -23,7 +23,7 @@ while (pimple.correctNonOrthogonal())
     (
         fvm::ddt(psi, p)
       + fvm::div(phid, p)
-      - fvm::laplacian(Dp, p)
+      - fvm::laplacian(rhorAUf, p)
     );
 
     pEqn.solve();
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
index 91becd8aa1..f8f99dafa5 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
@@ -13,7 +13,7 @@ surfaceScalarField phid
    *(
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+          //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )/fvc::interpolate(rho)
       - fvc::meshPhi(rho, U)
     )
@@ -25,7 +25,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     (
         fvm::ddt(psi, p)
       + fvm::div(phid, p)
-      - fvm::laplacian(Dp, p)
+      - fvm::laplacian(rhorAUf, p)
     );
 
     pEqn.solve();
diff --git a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
index 2a0b08074a..df3f3cdc61 100644
--- a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C
@@ -75,8 +75,12 @@ int main(int argc, char *argv[])
             // --- Pressure corrector loop
             while (pimple.correct())
             {
-                volScalarField rAU(1.0/UEqn.A());
-                surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+                volScalarField rAU("rAU", 1.0/UEqn.A());
+                surfaceScalarField rhorAUf
+                (
+                    "rhorAUf",
+                    fvc::interpolate(rho*rAU)
+                );
 
                 U = rAU*UEqn.H();
 
@@ -86,7 +90,7 @@ int main(int argc, char *argv[])
                     psi
                    *(
                        (fvc::interpolate(rho*U) & mesh.Sf())
-                     + Dp*fvc::ddtCorr(rho, U, phi)
+                     + rhorAUf*fvc::ddtCorr(rho, U, phi)
                     )/fvc::interpolate(rho)
                 );
 
@@ -97,7 +101,7 @@ int main(int argc, char *argv[])
                     fvm::ddt(psi, p)
                   + fvc::div(phi)
                   + fvm::div(phid, p)
-                  - fvm::laplacian(Dp, p)
+                  - fvm::laplacian(rhorAUf, p)
                 );
 
                 pEqn.solve();
diff --git a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
index 945aab6ba3..e9dedef5a9 100644
--- a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
+++ b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
@@ -93,7 +93,7 @@ int main(int argc, char *argv[])
             for (int corr=0; corr<nCorr; corr++)
             {
                 volScalarField rAU(1.0/UEqn.A());
-                surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+                surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
                 volVectorField HbyA("HbyA", U);
                 HbyA = rAU*UEqn.H();
@@ -102,14 +102,14 @@ int main(int argc, char *argv[])
                 (
                     "phiHbyA",
                     (fvc::interpolate(HbyA) & mesh.Sf())
-                  + Dp*fvc::ddtCorr(U, phi)
+                  + rAUf*fvc::ddtCorr(U, phi)
                 );
 
                 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
                 {
                     fvScalarMatrix pEqn
                     (
-                        fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
+                        fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
                     );
 
                     pEqn.setReference(pRefCell, pRefValue);
@@ -143,14 +143,14 @@ int main(int argc, char *argv[])
             BEqn.solve();
 
             volScalarField rAB(1.0/BEqn.A());
-            surfaceScalarField DpB("DpB", fvc::interpolate(rAB));
+            surfaceScalarField rABf("rABf", fvc::interpolate(rAB));
 
             phiB = (fvc::interpolate(B) & mesh.Sf())
-                + DpB*fvc::ddtCorr(B, phiB);
+                + rABf*fvc::ddtCorr(B, phiB);
 
             fvScalarMatrix pBEqn
             (
-                fvm::laplacian(DpB, pB) == fvc::div(phiB)
+                fvm::laplacian(rABf, pB) == fvc::div(phiB)
             );
             pBEqn.solve();
 
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
index 6168082e1d..43d8a4da53 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
@@ -51,6 +51,7 @@ Description
 #include "radiationModel.H"
 #include "fvIOoptionList.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
index 94299b84af..d77eb12b0e 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
@@ -1,25 +1,35 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rhok)*mesh.magSf());
+    surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
-      + Dp*fvc::ddtCorr(U, phi)
+      + rAUf*fvc::ddtCorr(U, phi)
       + phig
     );
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -36,7 +46,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
index b0d2f98aeb..483a723f7c 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
@@ -50,6 +50,7 @@ Description
 #include "RASModel.H"
 #include "fvIOoptionList.H"
 #include "simpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
index 53fec2e1b6..bae1bcef88 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn().A());
-    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn().H();
@@ -18,6 +18,16 @@
 
     phiHbyA += phig;
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index b1b6ceede2..9d62893f18 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -39,6 +39,7 @@ Description
 #include "radiationModel.H"
 #include "fvIOoptionList.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index ad98545b40..5f6d70e054 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -6,25 +6,36 @@
     thermo.rho() -= psi*p_rgh;
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
+    surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (
             (fvc::interpolate(rho*U) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rAUf*fvc::ddtCorr(rho, U, phi)
         )
       + phig
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+           *rho.boundaryField()
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
     fvScalarMatrix p_rghDDtEqn
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
@@ -38,7 +49,7 @@
         fvScalarMatrix p_rghEqn
         (
             p_rghDDtEqn
-          - fvm::laplacian(Dp, p_rgh)
+          - fvm::laplacian(rAUf, p_rgh)
         );
 
         fvOptions.constrain(p_rghEqn);
@@ -55,7 +66,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
             K = 0.5*magSqr(U);
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
index 67df5b61e0..2e52da8f53 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
@@ -36,6 +36,7 @@ Description
 #include "radiationModel.H"
 #include "simpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index 41c8900dc0..60b3c1a0cf 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -2,8 +2,8 @@
     rho = thermo.rho();
     rho.relax();
 
-    volScalarField rAU(1.0/UEqn().A());
-    surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
+    volScalarField rAU("rAU", 1.0/UEqn().A());
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn().H();
@@ -23,6 +23,17 @@
 
     phiHbyA += phig;
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+           *rho.boundaryField()
+        )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+    );
+
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
index b2cd9ac54e..6aae0cf5b6 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@@ -45,6 +45,7 @@ Description
 #include "radiationModel.H"
 #include "fvIOoptionList.H"
 #include "coordinateSystem.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
index a3615c4309..75720dd51c 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
@@ -38,6 +38,7 @@ Description
 #include "radiationModel.H"
 #include "fvIOoptionList.H"
 #include "coordinateSystem.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index f42ed7e39c..fc854c40a5 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@@ -4,8 +4,8 @@
     rho = min(rho, rhoMax[i]);
     rho.relax();
 
-    volScalarField rAU(1.0/UEqn().A());
-    surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
+    volScalarField rAU("rAU", 1.0/UEqn().A());
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn().H();
@@ -25,6 +25,17 @@
 
     phiHbyA += phig;
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+           *rho.boundaryField()
+        )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+    );
+
     dimensionedScalar compressibility = fvc::domainIntegrate(psi);
     bool compressible = (compressibility.value() > SMALL);
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index a712250f35..38ba49bd74 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -5,26 +5,37 @@
 
     rho = thermo.rho();
 
-    volScalarField rAU(1.0/UEqn().A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+    volScalarField rAU("rAU", 1.0/UEqn().A());
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn().H();
 
-    surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
       + phig
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+           *rho.boundaryField()
+        )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+    );
+
     {
         fvScalarMatrix p_rghDDtEqn
         (
@@ -41,7 +52,7 @@
             fvScalarMatrix p_rghEqn
             (
                 p_rghDDtEqn
-              - fvm::laplacian(Dp, p_rgh)
+              - fvm::laplacian(rhorAUf, p_rgh)
             );
 
             p_rghEqn.solve
@@ -63,7 +74,7 @@
             {
                 phi = phiHbyA + p_rghEqn.flux();
                 U = HbyA
-                  + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
+                  + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
                 U.correctBoundaryConditions();
                 fvOptions.correct(U);
                 K = 0.5*magSqr(U);
diff --git a/applications/solvers/incompressible/pimpleFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pEqn.H
index 6824d8ee0c..4799df09f7 100644
--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H
@@ -1,4 +1,4 @@
-surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
@@ -12,7 +12,7 @@ surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (fvc::interpolate(HbyA) & mesh.Sf())
-  + Dp*fvc::ddtCorr(U, phi)
+  + rAUf*fvc::ddtCorr(U, phi)
 );
 
 fvOptions.makeRelative(phiHbyA);
@@ -25,7 +25,7 @@ while (pimple.correctNonOrthogonal())
     // Pressure corrector
     fvScalarMatrix pEqn
     (
-        fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
+        fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
     );
 
     pEqn.setReference(pRefCell, pRefValue);
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createPcorrTypes.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createPcorrTypes.H
deleted file mode 100644
index a602fd4843..0000000000
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createPcorrTypes.H
+++ /dev/null
@@ -1,13 +0,0 @@
-    wordList pcorrTypes
-    (
-        p.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    for (label i=0; i<p.boundaryField().size(); i++)
-    {
-        if (p.boundaryField()[i].fixesValue())
-        {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
-        }
-    }
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
index 827daa4ce8..7ecfe1f662 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
@@ -1,4 +1,4 @@
-surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
@@ -12,7 +12,7 @@ surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (fvc::interpolate(HbyA) & mesh.Sf())
-  + Dp*fvc::ddtCorr(U, Uf)
+  + rAUf*fvc::ddtCorr(U, Uf)
 );
 
 if (p.needReference())
@@ -26,7 +26,7 @@ while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix pEqn
     (
-        fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
+        fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
     );
 
     pEqn.setReference(pRefCell, pRefValue);
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index 200ea289d8..fce4b6e64d 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
@@ -14,7 +14,7 @@ if (pimple.transonic())
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )/fvc::interpolate(rho)
     );
 
@@ -26,7 +26,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             coalParcels.Srho()
           + fvOptions(psi, p, rho.name())
@@ -49,7 +49,7 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
 
@@ -61,7 +61,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             coalParcels.Srho()
           + fvOptions(psi, p, rho.name())
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
index 1e504462b6..6652eab564 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H
@@ -1,25 +1,36 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 
-surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (
         (fvc::interpolate(rho*HbyA) & mesh.Sf())
-      + Dp*fvc::ddtCorr(rho, U, phi)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
     )
   + phig
 );
 
 fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+// Update the fixedFluxPressure BCs to ensure flux consistency
+setSnGrad<fixedFluxPressureFvPatchScalarField>
+(
+    p_rgh.boundaryField(),
+    (
+        phiHbyA.boundaryField()
+      - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+       *rho.boundaryField()
+    )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+);
+
 while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix p_rghEqn
@@ -27,7 +38,7 @@ while (pimple.correctNonOrthogonal())
         fvc::ddt(psi, rho)*gh
       + fvc::div(phiHbyA)
       + fvm::ddt(psi, p_rgh)
-      - fvm::laplacian(Dp, p_rgh)
+      - fvm::laplacian(rhorAUf, p_rgh)
      ==
         parcels.Srho()
       + surfaceFilm.Srho()
@@ -41,7 +52,7 @@ while (pimple.correctNonOrthogonal())
     if (pimple.finalNonOrthogonalIter())
     {
         phi = phiHbyA + p_rghEqn.flux();
-        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/Dp);
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
         U.correctBoundaryConditions();
         fvOptions.correct(U);
     }
diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
index 4067c125eb..27f272100a 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
@@ -39,6 +39,7 @@ Description
 #include "SLGThermo.H"
 #include "fvIOoptionList.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index 420bc1ef0c..56b22dba9c 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -6,7 +6,7 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -16,7 +16,7 @@
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
 
@@ -36,7 +36,7 @@
         fvScalarMatrix pEqn
         (
             pDDtEqn
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
         );
 
         fvOptions.constrain(pEqn);
diff --git a/applications/solvers/lagrangian/sprayFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/pEqn.H
index f62a17286c..26511b001d 100644
--- a/applications/solvers/lagrangian/sprayFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
@@ -14,7 +14,7 @@ if (pimple.transonic())
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )/fvc::interpolate(rho)
     );
 
@@ -26,7 +26,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             parcels.Srho()
           + fvOptions(psi, p, rho.name())
@@ -49,7 +49,7 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
 
@@ -61,7 +61,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             parcels.Srho()
           + fvOptions(psi, p, rho.name())
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
index b7de7e2f79..62c3e0376e 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
@@ -15,7 +15,7 @@ if (pimple.transonic())
        *(
             (
                 (fvc::interpolate(rho*HbyA) & mesh.Sf())
-              //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+              //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
             )/fvc::interpolate(rho)
           - fvc::meshPhi(rho, U)
         )
@@ -29,7 +29,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             parcels.Srho()
           + fvOptions(psi, p, rho.name())
@@ -52,7 +52,7 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(HbyA) & mesh.Sf())
-          //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+          //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
       - fvc::interpolate(rho)*fvc::meshPhi(rho, U)
     );
@@ -65,7 +65,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             parcels.Srho()
           + fvOptions(psi, p, rho.name())
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
index f869961a3c..e7f078d1a6 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
@@ -12,16 +12,16 @@
     surfaceScalarField rhof("rhof", fvc::interpolate(rho));
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
     phiv = (fvc::interpolate(HbyA) & mesh.Sf())
-         + Dp*fvc::ddtCorr(U, phivAbs);
+         + rhorAUf*fvc::ddtCorr(U, phivAbs);
     fvc::makeRelative(phiv, U);
 
-    surfaceScalarField phiGradp(Dp*mesh.magSf()*fvc::snGrad(p));
+    surfaceScalarField phiGradp(rhorAUf*mesh.magSf()*fvc::snGrad(p));
 
     phiv -= phiGradp/rhof;
 
@@ -35,7 +35,7 @@
           + psi*correction(fvm::ddt(p))
           + fvc::div(phiv, rho)
           + fvc::div(phiGradp)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
diff --git a/applications/solvers/multiphase/cavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
index e9d059aadb..210a30fd14 100644
--- a/applications/solvers/multiphase/cavitatingFoam/pEqn.H
+++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H
@@ -12,15 +12,15 @@
     surfaceScalarField rhof("rhof", fvc::interpolate(rho));
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
     phiv = (fvc::interpolate(HbyA) & mesh.Sf())
-         + Dp*fvc::ddtCorr(U, phiv);
+         + rhorAUf*fvc::ddtCorr(U, phiv);
 
-    surfaceScalarField phiGradp(Dp*mesh.magSf()*fvc::snGrad(p));
+    surfaceScalarField phiGradp(rhorAUf*mesh.magSf()*fvc::snGrad(p));
 
     phiv -= phiGradp/rhof;
 
@@ -32,7 +32,7 @@
           - (rhol0 + (psil - psiv)*pSat)*fvc::ddt(alphav) - pSat*fvc::ddt(psi)
           + fvc::div(phiv, rho)
           + fvc::div(phiGradp)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
index 7590e6ff45..a9c3365bdf 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C
@@ -46,6 +46,7 @@ Description
 #include "twoPhaseMixtureThermo.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,7 +62,7 @@ int main(int argc, char *argv[])
     #include "readControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
-    #include "createPcorrTypes.H"
+    #include "createPrghCorrTypes.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
index ed5faa7c86..030b997f48 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@@ -44,6 +44,7 @@ Description
 #include "twoPhaseMixtureThermo.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index adc380978c..dbd859b954 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -11,18 +11,27 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
     );
-    phi = phiHbyA;
 
     surfaceScalarField phig
     (
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
     );
 
     phiHbyA += phig;
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
     tmp<fvScalarMatrix> p_rghEqnComp1;
     tmp<fvScalarMatrix> p_rghEqnComp2;
 
@@ -70,7 +79,7 @@
         fvScalarMatrix p_rghEqnIncomp
         (
             fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p_rgh)
+          - fvm::laplacian(rAUf, p_rgh)
         );
 
         solve
@@ -97,7 +106,7 @@
             phi = phiHbyA + p_rghEqnIncomp.flux();
 
             U = HbyA
-              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/Dp);
+              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
index 3d1b93f429..d1e164e3cb 100644
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
@@ -46,6 +46,7 @@ Description
 #include "fvcSmooth.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
index b819068ddd..670ce727fd 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/MRFInterFoam.C
@@ -45,6 +45,7 @@ Description
 #include "IOMRFZoneList.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
index be0e723493..6c315c8eb8 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -13,23 +13,32 @@
     );
     adjustPhi(phiHbyA, U, p_rgh);
     mrfZones.makeRelative(phiHbyA);
-    phi = phiHbyA;
 
     surfaceScalarField phig
     (
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
     );
 
     phiHbyA += phig;
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - mrfZones.relative(mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -40,7 +49,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/createPcorrTypes.H b/applications/solvers/multiphase/interFoam/interDyMFoam/createPcorrTypes.H
deleted file mode 100644
index dfd4afb49b..0000000000
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/createPcorrTypes.H
+++ /dev/null
@@ -1,13 +0,0 @@
-    wordList pcorrTypes
-    (
-        p_rgh.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    for (label i=0; i<p_rgh.boundaryField().size(); i++)
-    {
-        if (p_rgh.boundaryField()[i].fixesValue())
-        {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
-        }
-    }
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/createUf.H b/applications/solvers/multiphase/interFoam/interDyMFoam/createUf.H
deleted file mode 100644
index aab92ce0e2..0000000000
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/createUf.H
+++ /dev/null
@@ -1,56 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 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/>.
-
-Global
-    createUf
-
-Description
-    Creates and initialises the velocity velocity field Uf.
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef createUf_H
-#define createUf_H
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-Info<< "Reading/calculating face velocity Uf\n" << endl;
-
-surfaceVectorField Uf
-(
-    IOobject
-    (
-        "Uf",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    linearInterpolate(U)
-);
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index 75e4c2f8d4..4c986ef2d8 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
@@ -41,6 +41,7 @@ Description
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,7 +57,7 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "createUf.H"
     #include "readTimeControls.H"
-    #include "createPcorrTypes.H"
+    #include "createPrghCorrTypes.H"
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
index 51190711c6..4fbd2666b4 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -19,23 +19,31 @@
         fvc::makeAbsolute(phiHbyA, U);
     }
 
-    surfaceScalarField phiAbs("phiAbs", phiHbyA);
-
     surfaceScalarField phig
     (
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
     );
 
     phiHbyA += phig;
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -46,7 +54,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C
index dccf3f8b28..fbcbd7c7e4 100644
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@@ -45,6 +45,7 @@ Description
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
index 1fb80a4b2c..1a10e45f3f 100644
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
@@ -38,6 +38,7 @@ Description
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H
index d8651bfa81..1194f97099 100644
--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -13,23 +13,32 @@
     );
 
     adjustPhi(phiHbyA, U, p_rgh);
-    phi = phiHbyA;
 
     surfaceScalarField phig
     (
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
     );
 
     phiHbyA += phig;
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -40,7 +49,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
diff --git a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
index 39c81c7a8d..4e6177ef9a 100644
--- a/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
+++ b/applications/solvers/multiphase/interFoam/porousInterFoam/porousInterFoam.C
@@ -47,6 +47,7 @@ Description
 #include "IOporosityModelList.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
index 1246ee203d..630b94534e 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
@@ -50,6 +50,7 @@ Description
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,9 +65,9 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
-    #include "../interFoam/interDyMFoam/createUf.H"
+    #include "createUf.H"
     #include "readTimeControls.H"
-    #include "../interFoam/interDyMFoam/createPcorrTypes.H"
+    #include "createPcorrTypes.H"
     #include "../interFoam/interDyMFoam/correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
index 1169de5c54..7863f5d748 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -19,18 +19,26 @@
         fvc::makeAbsolute(phiHbyA, U);
     }
 
-    surfaceScalarField phiAbs("phiAbs", phiHbyA);
-
     surfaceScalarField phig
     (
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
     );
 
     phiHbyA += phig;
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
     Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
     const volScalarField& vDotcP = vDotP[0]();
     const volScalarField& vDotvP = vDotP[1]();
@@ -39,7 +47,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvc::div(phiHbyA) - fvm::laplacian(Dp, p_rgh)
+            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
           - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
         );
 
@@ -51,7 +59,7 @@
         {
             phi = phiHbyA + p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
index bc6af055f0..3dd2b50fed 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C
@@ -48,6 +48,7 @@ Description
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
index 75f06420ee..6714fce91c 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -12,18 +12,27 @@
       + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
     );
     adjustPhi(phiHbyA, U, p_rgh);
-    phi = phiHbyA;
 
     surfaceScalarField phig
     (
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
     );
 
     phiHbyA += phig;
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
     Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
     const volScalarField& vDotcP = vDotP[0]();
     const volScalarField& vDotvP = vDotP[1]();
@@ -32,7 +41,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvc::div(phiHbyA) - fvm::laplacian(Dp, p_rgh)
+            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
           - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
         );
 
@@ -44,7 +53,7 @@
         {
             phi = phiHbyA + p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
index b509b53562..0c5c8a2a86 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam.C
@@ -35,12 +35,12 @@ Description
 #include "phaseModel.H"
 #include "dragModel.H"
 #include "heatTransferModel.H"
-#include "pimpleControl.H"
-
 #include "singlePhaseTransportModel.H"
 #include "LESModel.H"
-
+#include "pimpleControl.H"
 #include "IOMRFZoneList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index 00d41d3bb4..0f32807132 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -97,6 +97,8 @@
           + rAlphaAUfs[phasei]*fvc::ddtCorr(phase.U(), phase.phi())
         );
         mrfZones.makeRelative(phiHbyAs[phasei]);
+        mrfZones.makeRelative(phase.phi().oldTime());
+        mrfZones.makeRelative(phase.phi());
 
         phiHbyAs[phasei] +=
             rAlphaAUfs[phasei]
@@ -157,54 +159,58 @@
         phasei++;
     }
 
-    // Reset phi BCs
-    phi.boundaryField() = 0;
-
-    phasei = 0;
-    forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
-    {
-        phaseModel& phase = iter();
-
-        phase.phi().boundaryField() ==
-            (mesh.Sf().boundaryField() & phase.U().boundaryField());
-
-        mrfZones.makeRelative(phase.phi().oldTime());
-        mrfZones.makeRelative(phase.phi());
-
-        // Update phi BCs before pEqn
-        phi.boundaryField() +=
-            alphafs[phasei].boundaryField()*phase.phi().boundaryField();
-
-        phasei++;
-    }
-
-    surfaceScalarField Dp
+    surfaceScalarField rAUf
     (
         IOobject
         (
-            "Dp",
+            "rAUf",
             runTime.timeName(),
             mesh
         ),
         mesh,
-        dimensionedScalar("Dp", dimensionSet(-1, 3, 1, 0, 0), 0)
+        dimensionedScalar("rAUf", dimensionSet(-1, 3, 1, 0, 0), 0)
     );
 
     phasei = 0;
     forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     {
         phaseModel& phase = iter();
-        Dp += mag(alphafs[phasei]*rAlphaAUfs[phasei])/phase.rho();
+        rAUf += mag(alphafs[phasei]*rAlphaAUfs[phasei])/phase.rho();
 
         phasei++;
     }
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    {
+        surfaceScalarField::GeometricBoundaryField phib(phi.boundaryField());
+        phib = 0;
+        phasei = 0;
+        forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
+        {
+            phaseModel& phase = iter();
+
+            phib +=
+                alphafs[phasei].boundaryField()
+               *(mesh.Sf().boundaryField() & phase.U().boundaryField());
+
+            phasei++;
+        }
+
+        setSnGrad<fixedFluxPressureFvPatchScalarField>
+        (
+            p.boundaryField(),
+            (
+                phiHbyA.boundaryField() - mrfZones.relative(phib)
+            )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+        );
+    }
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqnIncomp
         (
             fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rAUf, p)
         );
 
         pEqnIncomp.setReference(pRefCell, pRefValue);
@@ -221,7 +227,7 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            surfaceScalarField mSfGradp(pEqnIncomp.flux()/Dp);
+            surfaceScalarField mSfGradp(pEqnIncomp.flux()/rAUf);
 
             phasei = 0;
             phi = dimensionedScalar("phi", phi.dimensions(), 0);
@@ -248,7 +254,7 @@
             // );
 
             p.relax();
-            mSfGradp = pEqnIncomp.flux()/Dp;
+            mSfGradp = pEqnIncomp.flux()/rAUf;
 
             U = dimensionedVector("U", dimVelocity, vector::zero);
 
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
index 1f4b4cdb01..bdcb810581 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/MRFMultiphaseInterFoam.C
@@ -37,6 +37,7 @@ Description
 #include "turbulenceModel.H"
 #include "IOMRFZoneList.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
index afa9d6058e..1286226676 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -13,20 +13,29 @@
     );
     adjustPhi(phiHbyA, U, p_rgh);
     mrfZones.makeRelative(phiHbyA);
-    phi = phiHbyA;
 
     surfaceScalarField phig
     (
-        - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
     );
 
     phiHbyA += phig;
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - mrfZones.relative(mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -37,7 +46,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
index 9f3fec616b..9bfeef4674 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
@@ -36,6 +36,7 @@ Description
 #include "multiphaseMixture.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
index d092aab3e2..a40c2be91b 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -12,23 +12,32 @@
       + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
     );
     adjustPhi(phiHbyA, U, p_rgh);
-    phi = phiHbyA;
 
     surfaceScalarField phig
     (
         (
             mixture.surfaceTensionForce()
           - ghf*fvc::snGrad(rho)
-        )*Dp*mesh.magSf()
+        )*rAUf*mesh.magSf()
     );
 
     phiHbyA += phig;
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -39,7 +48,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
index 9d5f21a94d..fac21e8b1c 100644
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
@@ -20,8 +20,15 @@ adjustPhi(phiHbyA, U, p_gh);
 
 fvOptions.makeRelative(phiHbyA);
 
-// Update the phi BCs from U before p BCs are updated
-phi.boundaryField() = mesh.Sf().boundaryField() & U.boundaryField();
+// Update the fixedFluxPressure BCs to ensure flux consistency
+setSnGrad<fixedFluxPressureFvPatchScalarField>
+(
+    p_gh.boundaryField(),
+    (
+        phiHbyA.boundaryField()
+      - (mesh.Sf().boundaryField() & U.boundaryField())
+    )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+);
 
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
diff --git a/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C b/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C
index ca35e81116..9b4caeed61 100644
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C
@@ -39,6 +39,7 @@ Description
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/settlingFoam/pEqn.H b/applications/solvers/multiphase/settlingFoam/pEqn.H
index 6bc6e9a6c1..13c3ddf7ff 100644
--- a/applications/solvers/multiphase/settlingFoam/pEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -10,23 +10,33 @@
         "phiHbyA",
         (
            (fvc::interpolate(rho*HbyA) & mesh.Sf())
-         + Dp*fvc::ddtCorr(rho, U, phi)
+         + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
-    phi = phiHbyA;
 
     surfaceScalarField phig
     (
-        - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
+        - ghf*fvc::snGrad(rho)*rhorAUf*mesh.magSf()
     );
 
     phiHbyA += phig;
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+           *rho.boundaryField()
+        )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(Dp, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
+            fvm::laplacian(rhorAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -37,7 +47,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/multiphase/settlingFoam/settlingFoam.C b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
index 57c531ba83..3b2454c8bf 100644
--- a/applications/solvers/multiphase/settlingFoam/settlingFoam.C
+++ b/applications/solvers/multiphase/settlingFoam/settlingFoam.C
@@ -38,6 +38,7 @@ Description
 #include "plasticViscosity.H"
 #include "yieldStress.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
index 94495f4f30..4a89f67d8f 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -12,20 +12,29 @@
       + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
     );
     adjustPhi(phiHbyA, U, p_rgh);
-    phi = phiHbyA;
 
     surfaceScalarField phig
     (
-        - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
     );
 
     phiHbyA += phig;
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -36,7 +45,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
         }
     }
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
index 8c51a402ee..247c12a8b5 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
@@ -37,6 +37,7 @@ Description
 #include "incompressibleTwoPhaseMixture.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
index 6485ee0323..9e99c5cc61 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
@@ -80,22 +80,12 @@
 
     surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
 
-    // Update phi BCs before pEqn
-    phi.boundaryField() =
-        mrfZones.relative
-        (
-            alpha1f.boundaryField()
-           *(mesh.Sf().boundaryField() & U1.boundaryField())
-          + alpha2f.boundaryField()
-           *(mesh.Sf().boundaryField() & U2.boundaryField())
-        );
-
     HbyA1 += (1.0/rho1)*rAU1*dragCoeff*U2;
     HbyA2 += (1.0/rho2)*rAU2*dragCoeff*U1;
 
-    surfaceScalarField Dp
+    surfaceScalarField rAUf
     (
-        "Dp",
+        "rAUf",
         mag
         (
             alpha1f*rAlphaAU1f/fvc::interpolate(rho1)
@@ -103,6 +93,22 @@
         )
     );
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - mrfZones.relative
+            (
+                alpha1f.boundaryField()
+               *(mesh.Sf().boundaryField() & U1.boundaryField())
+              + alpha2f.boundaryField()
+               *(mesh.Sf().boundaryField() & U2.boundaryField())
+            )
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
     tmp<fvScalarMatrix> pEqnComp1;
     tmp<fvScalarMatrix> pEqnComp2;
 
@@ -160,7 +166,7 @@
         fvScalarMatrix pEqnIncomp
         (
             fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rAUf, p)
         );
 
         solve
@@ -175,7 +181,7 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            surfaceScalarField mSfGradp(pEqnIncomp.flux()/Dp);
+            surfaceScalarField mSfGradp(pEqnIncomp.flux()/rAUf);
 
             phi1.boundaryField() ==
                 mrfZones.relative
@@ -200,7 +206,7 @@
             );
 
             p.relax();
-            mSfGradp = pEqnIncomp.flux()/Dp;
+            mSfGradp = pEqnIncomp.flux()/rAUf;
 
             U1 = HbyA1
               + fvc::reconstruct
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
index 7b8509337d..1e7cbafea7 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
@@ -40,6 +40,7 @@ Description
 #include "PhaseIncompressibleTurbulenceModel.H"
 #include "pimpleControl.H"
 #include "IOMRFZoneList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 2f0eaee518..6d8dabcec8 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -154,7 +154,6 @@ $(derivedFvPatchFields)/mappedFixedValue/mappedFixedValueFvPatchFields.C
 $(derivedFvPatchFields)/mappedFlowRate/mappedFlowRateFvPatchVectorField.C
 $(derivedFvPatchFields)/mappedVelocityFluxFixedValue/mappedVelocityFluxFixedValueFvPatchField.C
 $(derivedFvPatchFields)/movingWallVelocity/movingWallVelocityFvPatchVectorField.C
-$(derivedFvPatchFields)/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/oscillatingFixedValue/oscillatingFixedValueFvPatchFields.C
 $(derivedFvPatchFields)/outletInlet/outletInletFvPatchFields.C
 $(derivedFvPatchFields)/outletMappedUniformInlet/outletMappedUniformInletFvPatchFields.C
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/createRhoUf.H b/src/finiteVolume/cfdTools/compressible/createRhoUf.H
similarity index 93%
rename from applications/solvers/compressible/sonicFoam/sonicDyMFoam/createRhoUf.H
rename to src/finiteVolume/cfdTools/compressible/createRhoUf.H
index bb05d10f66..e91115ae75 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/createRhoUf.H
+++ b/src/finiteVolume/cfdTools/compressible/createRhoUf.H
@@ -22,15 +22,15 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Global
-    createUf
+    createRhoUf
 
 Description
-    Creates and initialises the velocity velocity field Uf.
+    Creates and initialises the velocity velocity field rhoUf.
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef createUf_H
-#define createUf_H
+#ifndef createRhoUf_H
+#define createRhoUf_H
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
index 70ce53da1e..8565193823 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
@@ -461,6 +461,28 @@ void Foam::MRFZone::makeRelative(volVectorField& U) const
 }
 
 
+void Foam::MRFZone::makeRelative(surfaceScalarField& phi) const
+{
+    makeRelativeRhoFlux(geometricOneField(), phi);
+}
+
+
+void Foam::MRFZone::makeRelative(FieldField<fvsPatchField, scalar>& phi) const
+{
+    makeRelativeRhoFlux(oneFieldField(), phi);
+}
+
+
+void Foam::MRFZone::makeRelative
+(
+    const surfaceScalarField& rho,
+    surfaceScalarField& phi
+) const
+{
+    makeRelativeRhoFlux(rho, phi);
+}
+
+
 void Foam::MRFZone::makeAbsolute(volVectorField& U) const
 {
     const volVectorField& C = mesh_.C();
@@ -499,28 +521,6 @@ void Foam::MRFZone::makeAbsolute(volVectorField& U) const
 }
 
 
-void Foam::MRFZone::makeRelative(surfaceScalarField& phi) const
-{
-    makeRelativeRhoFlux(geometricOneField(), phi);
-}
-
-
-void Foam::MRFZone::makeRelative(FieldField<fvsPatchField, scalar>& phi) const
-{
-    return makeRelativeRhoFlux(oneFieldField(), phi);
-}
-
-
-void Foam::MRFZone::makeRelative
-(
-    const surfaceScalarField& rho,
-    surfaceScalarField& phi
-) const
-{
-    makeRelativeRhoFlux(rho, phi);
-}
-
-
 void Foam::MRFZone::makeAbsolute(surfaceScalarField& phi) const
 {
     makeAbsoluteRhoFlux(geometricOneField(), phi);
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
index 97e59e2f6d..1396914122 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
@@ -218,9 +218,6 @@ public:
             //- Make the given absolute velocity relative within the MRF region
             void makeRelative(volVectorField& U) const;
 
-            //- Make the given relative velocity absolute within the MRF region
-            void makeAbsolute(volVectorField& U) const;
-
             //- Make the given absolute flux relative within the MRF region
             void makeRelative(surfaceScalarField& phi) const;
 
@@ -235,6 +232,9 @@ public:
                 surfaceScalarField& phi
             ) const;
 
+            //- Make the given relative velocity absolute within the MRF region
+            void makeAbsolute(volVectorField& U) const;
+
             //- Make the given relative flux absolute within the MRF region
             void makeAbsolute(surfaceScalarField& phi) const;
 
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C
index 8ab2aa1731..c56106f40c 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C
@@ -169,15 +169,6 @@ void Foam::MRFZoneList::makeRelative(volVectorField& U) const
 }
 
 
-void Foam::MRFZoneList::makeAbsolute(volVectorField& U) const
-{
-    forAll(*this, i)
-    {
-        operator[](i).makeAbsolute(U);
-    }
-}
-
-
 void Foam::MRFZoneList::makeRelative(surfaceScalarField& phi) const
 {
     forAll(*this, i)
@@ -193,7 +184,7 @@ Foam::MRFZoneList::relative
     const tmp<FieldField<fvsPatchField, scalar> >& phi
 ) const
 {
-    tmp<FieldField<fvsPatchField, scalar> > rphi(phi);
+    tmp<FieldField<fvsPatchField, scalar> > rphi(phi.ptr());
 
     forAll(*this, i)
     {
@@ -217,6 +208,15 @@ void Foam::MRFZoneList::makeRelative
 }
 
 
+void Foam::MRFZoneList::makeAbsolute(volVectorField& U) const
+{
+    forAll(*this, i)
+    {
+        operator[](i).makeAbsolute(U);
+    }
+}
+
+
 void Foam::MRFZoneList::makeAbsolute(surfaceScalarField& phi) const
 {
     forAll(*this, i)
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/createPcorrTypes.H b/src/finiteVolume/cfdTools/general/include/createPcorrTypes.H
similarity index 100%
rename from applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/createPcorrTypes.H
rename to src/finiteVolume/cfdTools/general/include/createPcorrTypes.H
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/createPcorrTypes.H b/src/finiteVolume/cfdTools/general/include/createPrghCorrTypes.H
similarity index 100%
rename from applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/createPcorrTypes.H
rename to src/finiteVolume/cfdTools/general/include/createPrghCorrTypes.H
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createUf.H b/src/finiteVolume/cfdTools/incompressible/createUf.H
similarity index 100%
rename from applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createUf.H
rename to src/finiteVolume/cfdTools/incompressible/createUf.H
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
index 12f59c1d11..0f21637953 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
@@ -38,11 +38,7 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
 )
 :
     fixedGradientFvPatchScalarField(p, iF),
-    phiHbyAName_("phiHbyA"),
-    phiName_("phi"),
-    rhoName_("rho"),
-    DpName_("Dp"),
-    adjoint_(false)
+    curTimeIndex_(-1)
 {}
 
 
@@ -55,11 +51,7 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
 )
 :
     fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
-    phiHbyAName_(ptf.phiHbyAName_),
-    phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_),
-    DpName_(ptf.DpName_),
-    adjoint_(ptf.adjoint_)
+    curTimeIndex_(-1)
 {}
 
 
@@ -71,11 +63,7 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
 )
 :
     fixedGradientFvPatchScalarField(p, iF),
-    phiHbyAName_(dict.lookupOrDefault<word>("phiHbyA", "phiHbyA")),
-    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
-    DpName_(dict.lookupOrDefault<word>("Dp", "Dp")),
-    adjoint_(dict.lookupOrDefault<Switch>("adjoint", false))
+    curTimeIndex_(-1)
 {
     if (dict.found("value") && dict.found("gradient"))
     {
@@ -99,11 +87,7 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
 )
 :
     fixedGradientFvPatchScalarField(wbppsf),
-    phiHbyAName_(wbppsf.phiHbyAName_),
-    phiName_(wbppsf.phiName_),
-    rhoName_(wbppsf.rhoName_),
-    DpName_(wbppsf.DpName_),
-    adjoint_(wbppsf.adjoint_)
+    curTimeIndex_(-1)
 {}
 
 
@@ -114,16 +98,29 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
 )
 :
     fixedGradientFvPatchScalarField(wbppsf, iF),
-    phiHbyAName_(wbppsf.phiHbyAName_),
-    phiName_(wbppsf.phiName_),
-    rhoName_(wbppsf.rhoName_),
-    DpName_(wbppsf.DpName_),
-    adjoint_(wbppsf.adjoint_)
+    curTimeIndex_(-1)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs
+(
+    const scalarField& snGradp
+)
+{
+    if (updated())
+    {
+        return;
+    }
+
+    curTimeIndex_ = this->db().time().timeIndex();
+
+    gradient() = snGradp;
+    fixedGradientFvPatchScalarField::updateCoeffs();
+}
+
+
 void Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs()
 {
     if (updated())
@@ -131,76 +128,19 @@ void Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const surfaceScalarField& phiHbyA =
-        db().lookupObject<surfaceScalarField>(phiHbyAName_);
-
-    const surfaceScalarField& phi =
-        db().lookupObject<surfaceScalarField>(phiName_);
-
-    fvsPatchField<scalar> phiHbyAp =
-        patch().patchField<surfaceScalarField, scalar>(phiHbyA);
-
-    fvsPatchField<scalar> phip =
-        patch().patchField<surfaceScalarField, scalar>(phi);
-
-    /*
-    if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
+    if (curTimeIndex_ != this->db().time().timeIndex())
     {
-        const fvPatchField<scalar>& rhop =
-            patch().lookupPatchField<volScalarField, scalar>(rhoName_);
-
-        phip /= rhop;
+        FatalErrorIn("fixedFluxPressureFvPatchScalarField::updateCoeffs()")
+            << "updateCoeffs(const scalarField& snGradp) MUST be called before"
+               " updateCoeffs() or evaluate() to set the boundary gradient."
+            << exit(FatalError);
     }
-
-    if (phiHbyA.dimensions() == dimDensity*dimVelocity*dimArea)
-    {
-        const fvPatchField<scalar>& rhop =
-            patch().lookupPatchField<volScalarField, scalar>(rhoName_);
-
-        phiHbyAp /= rhop;
-    }
-    */
-
-    const scalarField *DppPtr = NULL;
-
-    if (db().foundObject<volScalarField>(DpName_))
-    {
-        DppPtr =
-            &patch().lookupPatchField<volScalarField, scalar>(DpName_);
-    }
-    else if (db().foundObject<surfaceScalarField>(DpName_))
-    {
-        const surfaceScalarField& Dp =
-            db().lookupObject<surfaceScalarField>(DpName_);
-
-        DppPtr =
-            &patch().patchField<surfaceScalarField, scalar>(Dp);
-    }
-
-    if (adjoint_)
-    {
-        gradient() = (phip - phiHbyAp)/patch().magSf()/(*DppPtr);
-    }
-    else
-    {
-        gradient() = (phiHbyAp - phip)/patch().magSf()/(*DppPtr);
-    }
-
-    fixedGradientFvPatchScalarField::updateCoeffs();
 }
 
 
 void Foam::fixedFluxPressureFvPatchScalarField::write(Ostream& os) const
 {
     fixedGradientFvPatchScalarField::write(os);
-    writeEntryIfDifferent<word>(os, "phiHbyA", "phiHbyA", phiHbyAName_);
-    writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
-    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
-    writeEntryIfDifferent<word>(os, "Dp", "Dp", DpName_);
-    if (adjoint_)
-    {
-        os.writeKeyword("adjoint") << adjoint_ << token::END_STATEMENT << nl;
-    }
     writeEntry("value", os);
 }
 
@@ -216,4 +156,5 @@ namespace Foam
     );
 }
 
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.H
index e18ef462ea..2f43c6943f 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.H
@@ -28,44 +28,15 @@ Group
     grpInletBoundaryConditions grpWallBoundaryConditions
 
 Description
-    This boundary condition adjusts the pressure gradient such that the flux
-    on the boundary is that specified by the velocity boundary condition.
-
-    The predicted flux to be compensated by the pressure gradient is evaluated
-    as \f$(\phi - \phi_{H/A})\f$, both of which are looked-up from the database,
-    as is the pressure diffusivity used to calculate the gradient using:
-
-        \f[
-            \nabla(p) = \frac{\phi_{H/A} - \phi}{|Sf| D_p}
-        \f]
-
-    where
-
-    \vartable
-        phi     | flux
-        D_p     | pressure diffusivity
-        Sf      | patch face areas [m2]
-    \endvartable
-
-    \heading Patch usage
-
-    \table
-        Property     | Description             | Required    | Default value
-        phiHbyA      | name of predicted flux field | no | phiHbyA
-        phi          | name of flux field      | no |   phi
-        rho          | name of density field   | no |   rho
-        Dp           | name of pressure diffusivity field | no |   Dp
-    \endtable
+    This boundary condition sets the pressure gradient to the provided value
+    such that the flux on the boundary is that specified by the velocity
+    boundary condition.
 
     Example of the boundary condition specification:
     \verbatim
     myPatch
     {
         type            fixedFluxPressure;
-        phiHbyA         phiHbyA;
-        phi             phi;
-        rho             rho;
-        Dp              Dp;
     }
     \endverbatim
 
@@ -99,21 +70,8 @@ class fixedFluxPressureFvPatchScalarField
 {
     // Private data
 
-        //- Name of the predicted flux transporting the field
-        word phiHbyAName_;
-
-        //- Name of the flux transporting the field
-        word phiName_;
-
-        //- Name of the density field used to normalise the mass flux
-        //  if neccessary
-        word rhoName_;
-
-        //- Name of the pressure diffusivity field
-        word DpName_;
-
-        //- Is the pressure adjoint, i.e. has the opposite sign
-        Switch adjoint_;
+        //- Current time index (used for updating)
+        label curTimeIndex_;
 
 
 public:
@@ -186,7 +144,10 @@ public:
 
     // Member functions
 
-        //- Update the coefficients associated with the patch field
+        //- Update the patch pressure gradient field from the given snGradp
+        virtual void updateCoeffs(const scalarField& snGradp);
+
+        //- Update the patch pressure gradient field
         virtual void updateCoeffs();
 
         //- Write
@@ -200,6 +161,39 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "volFields.H"
+
+namespace Foam
+{
+    template<class GradBC>
+    inline void setSnGrad
+    (
+        volScalarField::GeometricBoundaryField& bf,
+        const FieldField<fvsPatchField, scalar>& snGrad
+    )
+    {
+        forAll(bf, patchi)
+        {
+            if (isA<GradBC>(bf[patchi]))
+            {
+                refCast<GradBC>(bf[patchi]).updateCoeffs(snGrad[patchi]);
+            }
+        }
+    }
+
+    template<class GradBC>
+    inline void setSnGrad
+    (
+        volScalarField::GeometricBoundaryField& bf,
+        const tmp<FieldField<fvsPatchField, scalar> >& tsnGrad
+    )
+    {
+        setSnGrad<GradBC>(bf, tsnGrad());
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C
deleted file mode 100644
index 89cf73e8ad..0000000000
--- a/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C
+++ /dev/null
@@ -1,197 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 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 "multiphaseFixedFluxPressureFvPatchScalarField.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "surfaceFields.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(p, iF),
-    phiHbyAName_("phiHbyA"),
-    phiName_("phi"),
-    rhoName_("rho"),
-    DpName_("Dp")
-{}
-
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
-    const multiphaseFixedFluxPressureFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
-    phiHbyAName_(ptf.phiHbyAName_),
-    phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_),
-    DpName_(ptf.DpName_)
-{}
-
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedGradientFvPatchScalarField(p, iF),
-    phiHbyAName_(dict.lookupOrDefault<word>("phiHbyA", "phiHbyA")),
-    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
-    DpName_(dict.lookupOrDefault<word>("Dp", "Dp"))
-{
-    if (dict.found("gradient"))
-    {
-        gradient() = scalarField("gradient", dict, p.size());
-        fixedGradientFvPatchScalarField::updateCoeffs();
-        fixedGradientFvPatchScalarField::evaluate();
-    }
-    else
-    {
-        fvPatchField<scalar>::operator=(patchInternalField());
-        gradient() = 0.0;
-    }
-}
-
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
-    const multiphaseFixedFluxPressureFvPatchScalarField& mfppsf
-)
-:
-    fixedGradientFvPatchScalarField(mfppsf),
-    phiHbyAName_(mfppsf.phiHbyAName_),
-    phiName_(mfppsf.phiName_),
-    rhoName_(mfppsf.rhoName_),
-    DpName_(mfppsf.DpName_)
-{}
-
-
-Foam::multiphaseFixedFluxPressureFvPatchScalarField::
-multiphaseFixedFluxPressureFvPatchScalarField
-(
-    const multiphaseFixedFluxPressureFvPatchScalarField& mfppsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(mfppsf, iF),
-    phiHbyAName_(mfppsf.phiHbyAName_),
-    phiName_(mfppsf.phiName_),
-    rhoName_(mfppsf.rhoName_),
-    DpName_(mfppsf.DpName_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::multiphaseFixedFluxPressureFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const surfaceScalarField& phiHbyA =
-        db().lookupObject<surfaceScalarField>(phiHbyAName_);
-
-    const surfaceScalarField& phi =
-        db().lookupObject<surfaceScalarField>(phiName_);
-
-    fvsPatchField<scalar> phiHbyAp =
-        patch().patchField<surfaceScalarField, scalar>(phiHbyA);
-
-    fvsPatchField<scalar> phip =
-        patch().patchField<surfaceScalarField, scalar>(phi);
-
-    /*
-    if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
-    {
-        const fvPatchField<scalar>& rhop =
-            patch().lookupPatchField<volScalarField, scalar>(rhoName_);
-
-        phip /= rhop;
-    }
-
-    if (phiHbyA.dimensions() == dimDensity*dimVelocity*dimArea)
-    {
-        const fvPatchField<scalar>& rhop =
-            patch().lookupPatchField<volScalarField, scalar>(rhoName_);
-
-        phiHbyAp /= rhop;
-    }
-    */
-
-    const fvsPatchField<scalar>& Dpp =
-        patch().lookupPatchField<surfaceScalarField, scalar>(DpName_);
-
-    gradient() = (phiHbyAp - phip)/patch().magSf()/Dpp;
-
-    fixedGradientFvPatchScalarField::updateCoeffs();
-}
-
-
-void Foam::multiphaseFixedFluxPressureFvPatchScalarField::write
-(
-    Ostream& os
-) const
-{
-    fvPatchScalarField::write(os);
-    writeEntryIfDifferent<word>(os, "phiHbyA", "phiHbyA", phiHbyAName_);
-    writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
-    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
-    writeEntryIfDifferent<word>(os, "Dp", "Dp", DpName_);
-    gradient().writeEntry("gradient", os);
-    writeEntry("value", os);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    makePatchTypeField
-    (
-        fvPatchScalarField,
-        multiphaseFixedFluxPressureFvPatchScalarField
-    );
-}
-
-// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H
deleted file mode 100644
index 18b9698c0e..0000000000
--- a/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H
+++ /dev/null
@@ -1,199 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 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::multiphaseFixedFluxPressureFvPatchScalarField
-
-Group
-    grpWallBoundaryConditions grpGenericBoundaryConditions
-
-Description
-    This boundary condition adjusts the pressure gradient such that the flux
-    on the boundary is that specified by the velocity boundary condition.
-
-    The predicted flux to be compensated by the pressure gradient is evaluated
-    as \f$(\phi - \phi_{H/A})\f$, both of which are looked-up from the database,
-    as is the pressure diffusivity Dp used to calculate the gradient using:
-
-        \f[
-            \nabla(p) = \frac{\phi_{H/A} - \phi}{|Sf| Dp}
-        \f]
-
-    where
-
-    \vartable
-        \phi    | flux
-        Dp      | pressure diffusivity
-        Sf      | patch face areas [m2]
-    \endvartable
-
-    \heading Patch usage
-
-    \table
-        Property     | Description             | Required    | Default value
-        phiHbyA      | name of predicted flux field | no | phiHbyA
-        phi          | name of flux field      | no |   phi
-        rho          | name of density field   | no |   rho
-        Dp           | name of pressure diffusivity field | no |   Dp
-    \endtable
-
-    Example of the boundary condition specification:
-    \verbatim
-    myPatch
-    {
-        type            multiphaseFixedFluxPressure;
-        phiHbyA         phiHbyA;
-        phi             phi;
-        rho             rho;
-        Dp              Dp;
-    }
-    \endverbatim
-
-SourceFiles
-    multiphaseFixedFluxPressureFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef multiphaseFixedFluxPressureFvPatchScalarFields_H
-#define multiphaseFixedFluxPressureFvPatchScalarFields_H
-
-#include "fvPatchFields.H"
-#include "fixedGradientFvPatchFields.H"
-#include "Switch.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-              Class multiphaseFixedFluxPressureFvPatch Declaration
-\*---------------------------------------------------------------------------*/
-
-class multiphaseFixedFluxPressureFvPatchScalarField
-:
-    public fixedGradientFvPatchScalarField
-{
-    // Private data
-
-        //- Name of the predicted flux transporting the field
-        word phiHbyAName_;
-
-        //- Name of the flux transporting the field
-        word phiName_;
-
-        //- Name of the density field used to normalise the mass flux
-        //  if neccessary
-        word rhoName_;
-
-        //- Name of the pressure diffusivity field
-        word DpName_;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("multiphaseFixedFluxPressure");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        multiphaseFixedFluxPressureFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        multiphaseFixedFluxPressureFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given
-        //  multiphaseFixedFluxPressureFvPatchScalarField onto a new patch
-        multiphaseFixedFluxPressureFvPatchScalarField
-        (
-            const multiphaseFixedFluxPressureFvPatchScalarField&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct as copy
-        multiphaseFixedFluxPressureFvPatchScalarField
-        (
-            const multiphaseFixedFluxPressureFvPatchScalarField&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new multiphaseFixedFluxPressureFvPatchScalarField(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        multiphaseFixedFluxPressureFvPatchScalarField
-        (
-            const multiphaseFixedFluxPressureFvPatchScalarField&,
-            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 multiphaseFixedFluxPressureFvPatchScalarField(*this, iF)
-            );
-        }
-
-
-    // Member 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/fvOptions/fvOptions/fvOption.C b/src/fvOptions/fvOptions/fvOption.C
index 0751b64c77..6a34f9ae07 100644
--- a/src/fvOptions/fvOptions/fvOption.C
+++ b/src/fvOptions/fvOptions/fvOption.C
@@ -27,6 +27,7 @@ License
 #include "fvMesh.H"
 #include "fvMatrices.H"
 #include "volFields.H"
+#include "fvsPatchFields.H"
 #include "ListOps.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -480,6 +481,15 @@ void Foam::fv::option::makeRelative(surfaceScalarField& phi) const
 }
 
 
+void Foam::fv::option::makeRelative
+(
+    FieldField<fvsPatchField, scalar>& phi
+) const
+{
+    // do nothing
+}
+
+
 void Foam::fv::option::makeRelative
 (
     const surfaceScalarField& rho,
diff --git a/src/fvOptions/fvOptions/fvOption.H b/src/fvOptions/fvOptions/fvOption.H
index 0f8896827b..ec19296565 100644
--- a/src/fvOptions/fvOptions/fvOption.H
+++ b/src/fvOptions/fvOptions/fvOption.H
@@ -427,6 +427,12 @@ public:
                 //- Make the given absolute flux relative
                 virtual void makeRelative(surfaceScalarField& phi) const;
 
+                //- Make the given absolute boundary flux relative
+                virtual void makeRelative
+                (
+                    FieldField<fvsPatchField, scalar>& phi
+                ) const;
+
                 //- Make the given absolute mass-flux relative
                 virtual void makeRelative
                 (
diff --git a/src/fvOptions/fvOptions/fvOptionList.C b/src/fvOptions/fvOptions/fvOptionList.C
index 5329aa1243..b3ccc28f87 100644
--- a/src/fvOptions/fvOptions/fvOptionList.C
+++ b/src/fvOptions/fvOptions/fvOptionList.C
@@ -26,6 +26,7 @@ License
 #include "fvOptionList.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
+#include "surfaceFields.H"
 #include "Time.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -129,6 +130,23 @@ void Foam::fv::optionList::makeRelative
 }
 
 
+Foam::tmp<Foam::FieldField<Foam::fvsPatchField, Foam::scalar> >
+Foam::fv::optionList::relative
+(
+    const tmp<FieldField<fvsPatchField, scalar> >& phi
+) const
+{
+    tmp<FieldField<fvsPatchField, scalar> > rphi(phi.ptr());
+
+    forAll(*this, i)
+    {
+        operator[](i).makeRelative(rphi());
+    }
+
+    return rphi;
+}
+
+
 void Foam::fv::optionList::makeAbsolute(surfaceScalarField& phi) const
 {
     forAll(*this, i)
diff --git a/src/fvOptions/fvOptions/fvOptionList.H b/src/fvOptions/fvOptions/fvOptionList.H
index 0d5b6abdb5..15e9733108 100644
--- a/src/fvOptions/fvOptions/fvOptionList.H
+++ b/src/fvOptions/fvOptions/fvOptionList.H
@@ -166,6 +166,12 @@ public:
                 surfaceScalarField& phi
             ) const;
 
+            //- Return the given absolute boundary flux relative
+            tmp<FieldField<fvsPatchField, scalar> > relative
+            (
+                const tmp<FieldField<fvsPatchField, scalar> >& tphi
+            ) const;
+
             //- Make the given relative flux absolute
             void makeAbsolute(surfaceScalarField& phi) const;
 
diff --git a/src/fvOptions/sources/derived/MRFSource/MRFSource.C b/src/fvOptions/sources/derived/MRFSource/MRFSource.C
index d61e10f362..9f9fcca442 100644
--- a/src/fvOptions/sources/derived/MRFSource/MRFSource.C
+++ b/src/fvOptions/sources/derived/MRFSource/MRFSource.C
@@ -126,6 +126,15 @@ void Foam::fv::MRFSource::makeRelative(surfaceScalarField& phi) const
 }
 
 
+void Foam::fv::MRFSource::makeRelative
+(
+    FieldField<fvsPatchField, scalar>& phi
+) const
+{
+    mrfPtr_->makeRelative(phi);
+}
+
+
 void Foam::fv::MRFSource::makeRelative
 (
     const surfaceScalarField& rho,
diff --git a/src/fvOptions/sources/derived/MRFSource/MRFSource.H b/src/fvOptions/sources/derived/MRFSource/MRFSource.H
index 0b3da28cc4..8f99788bc9 100644
--- a/src/fvOptions/sources/derived/MRFSource/MRFSource.H
+++ b/src/fvOptions/sources/derived/MRFSource/MRFSource.H
@@ -140,6 +140,12 @@ public:
             //- Make the given absolute flux relative
             virtual void makeRelative(surfaceScalarField& phi) const;
 
+            //- Make the given absolute boundary flux relative
+            virtual void makeRelative
+            (
+                FieldField<fvsPatchField, scalar>& phi
+            ) const;
+
             //- Make the given absolute mass-flux relative
             virtual void makeRelative
             (
diff --git a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/0/p_rgh.org b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/0/p_rgh.org
index 396d380dad..8bd9793d9a 100644
--- a/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/0/p_rgh.org
+++ b/tutorials/multiphase/compressibleInterFoam/laminar/depthCharge2D/0/p_rgh.org
@@ -23,7 +23,6 @@ boundaryField
     walls
     {
         type            fixedFluxPressure;
-        value           uniform 1e5;
     }
 
     defaultFaces
diff --git a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/p_rgh b/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/p_rgh
index 85bfb4aef1..1b9491f31e 100644
--- a/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/p_rgh
+++ b/tutorials/multiphase/interDyMFoam/ras/damBreakWithObstacle/constant/p_rgh
@@ -30,14 +30,3 @@ boundaryField
         p0              uniform 0;
         value           uniform 0;
     }
-    walls
-    {
-        type            fixedFluxPressure;
-        gradient        uniform 0;
-        phi             phiAbs;
-        value           uniform 0;
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/0/p_rgh b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/0/p_rgh
index dafcfc526a..66122dcff0 100644
--- a/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/0/p_rgh
+++ b/tutorials/multiphase/interDyMFoam/ras/sloshingTank2D3DoF/0/p_rgh
@@ -31,8 +31,6 @@ boundaryField
     walls
     {
         type            fixedFluxPressure;
-        phi             phiAbs;
-        value           uniform 0;
     }
 }