diff --git a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
index bd5e2c3b99d3779fdd0ffce94cf68ccecfe8d657..43eb616e17376390e190ed965c566a3ec9a193db 100644
--- a/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
+++ b/applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C
@@ -48,6 +48,7 @@ int main(int argc, char *argv[])
     #include "createEngineMesh.H"
     #include "createFields.H"
     #include "createFvOptions.H"
+    #include "createRhoUf.H"
     #include "initContinuityErrs.H"
     #include "readEngineTimeControls.H"
     #include "compressibleCourantNo.H"
diff --git a/applications/solvers/combustion/engineFoam/engineFoam.C b/applications/solvers/combustion/engineFoam/engineFoam.C
index 592abbd3e67ba23a79d05c4432ed1452af46b1d6..a35e22dc58d9cb2595766221d30ed1b778e1cf16 100644
--- a/applications/solvers/combustion/engineFoam/engineFoam.C
+++ b/applications/solvers/combustion/engineFoam/engineFoam.C
@@ -72,6 +72,7 @@ int main(int argc, char *argv[])
     #include "readCombustionProperties.H"
     #include "createFields.H"
     #include "createFvOptions.H"
+    #include "createRhoUf.H"
     #include "initContinuityErrs.H"
     #include "readEngineTimeControls.H"
     #include "compressibleCourantNo.H"
diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index 384ecbee598a7f03806862c4c74645f1ba142c68..2f6333cc9761b89834245073ac3edfad744a318f 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/applications/solvers/combustion/engineFoam/pEqn.H
@@ -15,12 +15,12 @@ if (pimple.transonic())
        *(
             (
                 (fvc::interpolate(rho*HbyA) & mesh.Sf())
-              //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
+              + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
             )/fvc::interpolate(rho)
-          - fvc::meshPhi(rho, U)
         )
     );
 
+    fvc::makeRelative(phid, psi, U);
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
@@ -51,11 +51,11 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
         )
-      - fvc::interpolate(rho)*fvc::meshPhi(rho, U)
     );
 
+    fvc::makeRelative(phiHbyA, rho, U);
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
@@ -88,6 +88,15 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+{
+    rhoUf = fvc::interpolate(rho*U);
+    surfaceVectorField n(mesh.Sf()/mesh.magSf());
+    rhoUf +=
+        mesh.Sf()
+       *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf))
+       /sqr(mesh.magSf());
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p) - fvc::div(fvc::meshPhi(rho, U), p);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index 40cf109457916f765b7b56e4ebb1414e900862bf..14ba97d1c2dac1c2a74a6a94de1e1d6b172c5197 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -26,6 +26,7 @@ if (pimple.transonic())
         )/fvc::interpolate(rho)
     );
 
+    fvc::makeRelative(phid, psi, U);
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
@@ -58,6 +59,7 @@ else
       + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
     );
 
+    fvc::makeRelative(phiHbyA, rho, U);
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
@@ -102,6 +104,15 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+{
+    rhoUf = fvc::interpolate(rho*U);
+    surfaceVectorField n(mesh.Sf()/mesh.magSf());
+    rhoUf +=
+        mesh.Sf()
+       *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf))
+       /sqr(mesh.magSf());
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
index aae32851863f83ce2f6cb9128fa6b7946a1377a9..e633ebee602bd230152bd201135a052a55c02e90 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
@@ -76,6 +76,9 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         {
+            // Store momentum to set rhoUf for introduced faces.
+            volVectorField rhoU("rhoU", rho*U);
+
             // Store divrhoU from the previous time-step/mesh for the correctPhi
             volScalarField divrhoU(fvc::div(fvc::absolute(phi, rho, U)));
 
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
index f8f99dafa545fa0758ce83ece4778e51083b0a39..65afd1aeb88642cae2b880ac21ba6951bab763f0 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
@@ -13,12 +13,13 @@ surfaceScalarField phid
    *(
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
         )/fvc::interpolate(rho)
-      - fvc::meshPhi(rho, U)
     )
 );
 
+fvc::makeRelative(phid, psi, U);
+
 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
     fvScalarMatrix pEqn
@@ -37,3 +38,12 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
+
+{
+    rhoUf = fvc::interpolate(rho*U);
+    surfaceVectorField n(mesh.Sf()/mesh.magSf());
+    rhoUf +=
+        mesh.Sf()
+       *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf))
+       /sqr(mesh.magSf());
+}
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
index 11ec1de915e09c0d7c93fce3e04e9b7b652fe28b..62d4716fb5a70691651f5304e4fb4adcfabf165a 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C
@@ -64,12 +64,6 @@ int main(int argc, char *argv[])
 
         mesh.movePoints(motionPtr->newPoints());
 
-        // Calculate absolute flux from the mapped surface velocity
-        phi = mesh.Sf() & rhoUf;
-
-        // Make the flux relative to the mesh motion
-        fvc::makeRelative(phi, rho, U);
-
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
index 62c3e0376e363ef3f9a8d9a459d33bf10773dc33..b10c24aa0b7c86eb281cea29e612929b31e1ad50 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
@@ -15,12 +15,12 @@ if (pimple.transonic())
        *(
             (
                 (fvc::interpolate(rho*HbyA) & mesh.Sf())
-              //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
+              + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
             )/fvc::interpolate(rho)
-          - fvc::meshPhi(rho, U)
         )
     );
 
+    fvc::makeRelative(phid, psi, U);
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
 
     while (pimple.correctNonOrthogonal())
@@ -52,11 +52,11 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(HbyA) & mesh.Sf())
-          //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
         )
-      - fvc::interpolate(rho)*fvc::meshPhi(rho, U)
     );
 
+    fvc::makeRelative(phiHbyA, rho, U);
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
@@ -90,6 +90,15 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+{
+    rhoUf = fvc::interpolate(rho*U);
+    surfaceVectorField n(mesh.Sf()/mesh.magSf());
+    rhoUf +=
+        mesh.Sf()
+       *(fvc::absolute(phi, rho, U) - (mesh.Sf() & rhoUf))
+       /sqr(mesh.magSf());
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p) - fvc::div(fvc::meshPhi(rho, U), p);
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C
index baafc8d5e3ff9625d987d0ae1c7d6d1a3e37abb7..b91a3ad79de5e5c8341e571564c65b8453994ea9 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C
+++ b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C
@@ -52,6 +52,7 @@ int main(int argc, char *argv[])
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createFvOptions.H"
+    #include "createRhoUf.H"
     #include "createClouds.H"
     #include "createRadiationModel.H"
     #include "initContinuityErrs.H"