diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index fc700f81fa..9e67a47c4f 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -2,8 +2,8 @@
     rho1 = eos1->rho(p, T);
     rho2 = eos2->rho(p, T);
 
-    volScalarField rAU = 1.0/UEqn.A();
-    surfaceScalarField rAUf = fvc::interpolate(rAU);
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
     tmp<fvScalarMatrix> p_rghEqnComp1;
     tmp<fvScalarMatrix> p_rghEqnComp2;
@@ -27,20 +27,25 @@
           - fvc::Sp(fvc::div(phid2), p_rgh);
     }
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
     (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    phi = phiU +
+    surfaceScalarField phig
+    (
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf();
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
 
     // Thermodynamic density needs to be updated by psi*d(p) after the
     // pressure solution - done in 2 parts. Part 1:
@@ -50,7 +55,7 @@
     {
         fvScalarMatrix p_rghEqnIncomp
         (
-            fvc::div(phi)
+            fvc::div(phiHbyA)
           - fvm::laplacian(rAUf, p_rgh)
         );
 
@@ -74,13 +79,15 @@
                 pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
               - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
             );
-            phi += p_rghEqnIncomp.flux();
+
+            phi = phiHbyA + p_rghEqnIncomp.flux();
+
+            U = HbyA
+              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
     p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
 
     rho1 = eos1->rho(p, T);
diff --git a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
index 13d6b60662..9c56c4d26f 100644
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
-    volScalarField rAU("Dp", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dpf", fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -11,7 +11,6 @@
         (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-
     mrfZones.relativeFlux(phiHbyA);
 
     adjustPhi(phiHbyA, U, p_rgh);
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
index c18d5aba41..f6681086d0 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
-    volScalarField rAU("Dp", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dpf", fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H
index 935ed02c15..d1bc7bf1ed 100644
--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@@ -1,6 +1,6 @@
 {
-    volScalarField rAU("Dp", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dpf", fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
index c9b65eb9d5..e3ef0dc9ef 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@@ -1,24 +1,28 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
     (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);
 
-    phi =
-        phiU
-      + (
+    surfaceScalarField phig
+    (
+        (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf();
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
 
     Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
     const volScalarField& vDotcP = vDotP[0]();
@@ -28,7 +32,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvc::div(phi) - fvm::laplacian(rAUf, p_rgh)
+            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
           - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
         );
 
@@ -38,13 +42,13 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += p_rghEqn.flux();
+            phi = phiHbyA + p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
     #include "continuityErrs.H"
 
     p == p_rgh + rho*gh;
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
index efb181f3f2..7bc1eeb12c 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
@@ -1,31 +1,32 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
     (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
-      //+ fvc::ddtPhiCorr(rAU, rho, U, phi)
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
-    mrfZones.relativeFlux(phiU);
+    mrfZones.relativeFlux(phiHbyA);
 
-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);
 
-    phi =
-        phiU
-      + (
-            mixture.surfaceTensionForce()
-          - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf();
+    surfaceScalarField phig
+    (
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -34,13 +35,13 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
     #include "continuityErrs.H"
 
     p == p_rgh + rho*gh;
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
index e7e5bc14cf..c1e346ce43 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
@@ -1,29 +1,34 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
     (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);
 
-    phi = phiU +
+    surfaceScalarField phig
     (
-        mixture.surfaceTensionForce()
-      - ghf*fvc::snGrad(rho)
-    )*rAUf*mesh.magSf();
+        (
+            mixture.surfaceTensionForce()
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -32,13 +37,13 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
     #include "continuityErrs.H"
 
     p == p_rgh + rho*gh;
diff --git a/applications/solvers/multiphase/settlingFoam/pEqn.H b/applications/solvers/multiphase/settlingFoam/pEqn.H
index d89fd3ef4e..bf632adb4d 100644
--- a/applications/solvers/multiphase/settlingFoam/pEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/pEqn.H
@@ -1,54 +1,60 @@
-volScalarField rAU(1.0/UEqn.A());
-
-surfaceScalarField rAUf
-(
-    "(rho*(1|A(U)))",
-    fvc::interpolate(rho)*fvc::interpolate(rAU)
-);
-
-U = rAU*UEqn.H();
-phi =
-    fvc::interpolate(rho)
-   *(
-       (fvc::interpolate(U) & mesh.Sf())
-     + fvc::ddtPhiCorr(rAU, rho, U, phi)
-    );
-
-surfaceScalarField phiU("phiU", phi);
-phi -= ghf*fvc::snGrad(rho)*rAUf*mesh.magSf();
-
-while (pimple.correctNonOrthogonal())
 {
-    fvScalarMatrix p_rghEqn
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rho*rAU));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
     (
-        fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phi)
+        fvc::interpolate(rho)
+       *(
+           (fvc::interpolate(HbyA) & mesh.Sf())
+         + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        )
     );
 
-    p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+    surfaceScalarField phig
+    (
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+    );
 
-    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+    phiHbyA += phig;
 
-    if (pimple.finalNonOrthogonalIter())
+    while (pimple.correctNonOrthogonal())
     {
-        phi -= p_rghEqn.flux();
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
+        }
     }
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+
+    #include "rhoEqn.H"
+    #include "compressibleContinuityErrs.H"
 }
-
-p == p_rgh + rho*gh;
-
-if (p_rgh.needReference())
-{
-    p += dimensionedScalar
-    (
-        "p",
-        p.dimensions(),
-        pRefValue - getRefCellValue(p, pRefCell)
-    );
-    p_rgh = p - rho*gh;
-}
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
-
-U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-U.correctBoundaryConditions();
diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
index ac7fc35f68..ddfca4e3ea 100644
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
@@ -1,24 +1,31 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
-    surfaceScalarField phiU
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
     (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);
 
-    phi = phiU - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf();
+    surfaceScalarField phig
+    (
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -27,13 +34,13 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
         }
     }
 
-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
     #include "continuityErrs.H"
 
     p == p_rgh + rho*gh;
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
index aaa6ff8bde..ebfc0a8386 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
@@ -41,6 +41,7 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
     phiHbyAName_("phiHbyA"),
     phiName_("phi"),
     rhoName_("rho"),
+    DpName_("Dp"),
     adjoint_(false)
 {}
 
@@ -57,6 +58,7 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
     phiHbyAName_(ptf.phiHbyAName_),
     phiName_(ptf.phiName_),
     rhoName_(ptf.rhoName_),
+    DpName_(ptf.rhoName_),
     adjoint_(ptf.adjoint_)
 {}
 
@@ -72,6 +74,7 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
     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))
 {
     if (dict.found("gradient"))
@@ -97,6 +100,7 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
     phiHbyAName_(wbppsf.phiHbyAName_),
     phiName_(wbppsf.phiName_),
     rhoName_(wbppsf.rhoName_),
+    DpName_(wbppsf.DpName_),
     adjoint_(wbppsf.adjoint_)
 {}
 
@@ -111,6 +115,7 @@ Foam::fixedFluxPressureFvPatchScalarField::fixedFluxPressureFvPatchScalarField
     phiHbyAName_(wbppsf.phiHbyAName_),
     phiName_(wbppsf.phiName_),
     rhoName_(wbppsf.rhoName_),
+    DpName_(wbppsf.DpName_),
     adjoint_(wbppsf.adjoint_)
 {}
 
@@ -136,6 +141,7 @@ void Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs()
     fvsPatchField<scalar> phip =
         patch().patchField<surfaceScalarField, scalar>(phi);
 
+    /*
     if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
     {
         const fvPatchField<scalar>& rhop =
@@ -144,16 +150,38 @@ void Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs()
         phip /= rhop;
     }
 
-    const fvPatchField<scalar>& Dpp =
-        patch().lookupPatchField<volScalarField, scalar>("Dp");
+    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()/Dpp;
+        gradient() = (phip - phiHbyAp)/patch().magSf()/(*DppPtr);
     }
     else
     {
-        gradient() = (phiHbyAp - phip)/patch().magSf()/Dpp;
+        gradient() = (phiHbyAp - phip)/patch().magSf()/(*DppPtr);
     }
 
     fixedGradientFvPatchScalarField::updateCoeffs();
@@ -166,6 +194,7 @@ void Foam::fixedFluxPressureFvPatchScalarField::write(Ostream& os) const
     writeEntryIfDifferent<word>(os, "phiHbyA", "phiHbyA", phiHbyAName_);
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
+    writeEntryIfDifferent<word>(os, "Dp", "Dp", rhoName_);
     os.writeKeyword("adjoint") << adjoint_ << token::END_STATEMENT << nl;
     gradient().writeEntry("gradient", os);
 }
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.H
index 82a7f8a6b4..de8f837abd 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.H
@@ -71,6 +71,9 @@ class fixedFluxPressureFvPatchScalarField
         //  if neccessary
         word rhoName_;
 
+        //- Name of the pressure diffusivity field
+        word DpName_;
+
         //- Is the pressure adjoint, i.e. has the opposite sign
         Switch adjoint_;
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C
index aa84755155..644c8b161c 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C
@@ -136,6 +136,7 @@ void Foam::multiphaseFixedFluxPressureFvPatchScalarField::updateCoeffs()
     fvsPatchField<scalar> phip =
         patch().patchField<surfaceScalarField, scalar>(phi);
 
+    /*
     if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
     {
         const fvPatchField<scalar>& rhop =
@@ -144,6 +145,15 @@ void Foam::multiphaseFixedFluxPressureFvPatchScalarField::updateCoeffs()
         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>("Dp");
 
diff --git a/tutorials/multiphase/settlingFoam/ras/dahl/system/fvSchemes b/tutorials/multiphase/settlingFoam/ras/dahl/system/fvSchemes
index bfdc69ada0..026c549f6c 100644
--- a/tutorials/multiphase/settlingFoam/ras/dahl/system/fvSchemes
+++ b/tutorials/multiphase/settlingFoam/ras/dahl/system/fvSchemes
@@ -39,7 +39,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p_rgh) Gauss linear corrected;
+    laplacian(Dp,p_rgh) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian(mut,Alpha) Gauss linear corrected;
diff --git a/tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSchemes b/tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSchemes
index 65b5f239ef..ffdb15771e 100644
--- a/tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSchemes
+++ b/tutorials/multiphase/settlingFoam/ras/tank3D/system/fvSchemes
@@ -39,7 +39,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p_rgh) Gauss linear corrected;
+    laplacian(Dp,p_rgh) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian(mut,Alpha) Gauss linear corrected;