diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
index 9c8c2d4f18c2bbdf4267a19af53d7d2ecb690f37..6dc67ce353c5ca4b2080cac6cf28ef4663084492 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H
@@ -159,11 +159,10 @@
     forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
     {
         phaseModel& phase = iter();
-        Dp += alphafs[phasei]*rAlphaAUfs[phasei]/phase.rho();
+        Dp += mag(alphafs[phasei]*rAlphaAUfs[phasei])/phase.rho();
 
         phasei++;
     }
-    Dp = mag(Dp);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -173,6 +172,8 @@
           - fvm::laplacian(Dp, p)
         );
 
+        pEqnIncomp.setReference(pRefCell, pRefValue);
+
         solve
         (
             // (
@@ -196,7 +197,10 @@
                 phase.phi() =
                     phiHbyAs[phasei]
                   + rAlphaAUfs[phasei]*mSfGradp/phase.rho();
-                phi += alphafs[phasei]*phase.phi();
+                phi +=
+                    alphafs[phasei]*phiHbyAs[phasei]
+                  + mag(alphafs[phasei]*rAlphaAUfs[phasei])
+                   *mSfGradp/phase.rho();
 
                 phasei++;
             }