DyM solvers: correct Uf using phi after construction

applications/solvers/combustion/engineFoam/pEqn.H

@@ -91,10 +91,7 @@ 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());
+    rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf));
 }
 
 if (thermo.dpdt())

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H

@@ -107,10 +107,7 @@ 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());
+    rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf));
 }
 
 if (thermo.dpdt())

applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H

@@ -42,8 +42,5 @@ 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());
+    rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf));
 }

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H

@@ -51,7 +51,7 @@ fvOptions.correct(U);
 {
     Uf = fvc::interpolate(U);
     surfaceVectorField n(mesh.Sf()/mesh.magSf());
-    Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+    Uf += n*(phi/mesh.magSf() - (n & Uf));
 }
 
 // Make the fluxes relative to the mesh motion

applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H

@@ -93,10 +93,7 @@ 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());
+    rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf));
 }
 
 if (thermo.dpdt())

applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H

@@ -85,6 +85,6 @@
     {
         Uf = fvc::interpolate(U);
         surfaceVectorField n(mesh.Sf()/mesh.magSf());
-        Uf += mesh.Sf()*(phiv - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+        Uf += n*(phiv/mesh.magSf() - (n & Uf));
     }
 }

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H

@@ -114,7 +114,7 @@
     {
         Uf = fvc::interpolate(U);
         surfaceVectorField n(mesh.Sf()/mesh.magSf());
-        Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+        Uf += n*(phi/mesh.magSf() - (n & Uf));
     }
 
     // Make the fluxes relative to the mesh motion

applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C

@@ -55,10 +55,10 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
-    #include "createUf.H"
     #include "readTimeControls.H"
     #include "createPrghCorrTypes.H"
     #include "correctPhi.H"
+    #include "createUf.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 

applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H

@@ -67,7 +67,7 @@
     {
         Uf = fvc::interpolate(U);
         surfaceVectorField n(mesh.Sf()/mesh.magSf());
-        Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+        Uf += n*(phi/mesh.magSf() - (n & Uf));
     }
 
     // Make the fluxes relative to the mesh motion

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C

@@ -65,10 +65,10 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
-    #include "createUf.H"
     #include "readTimeControls.H"
     #include "createPcorrTypes.H"
     #include "../interFoam/interDyMFoam/correctPhi.H"
+    #include "createUf.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H

@@ -68,7 +68,7 @@
     {
         Uf = fvc::interpolate(U);
         surfaceVectorField n(mesh.Sf()/mesh.magSf());
-        Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+        Uf += n*(phi/mesh.magSf() - (n & Uf));
     }
 
     // Make the fluxes relative to the mesh motion

src/finiteVolume/cfdTools/compressible/createRhoUf.H

@@ -46,9 +46,14 @@ surfaceVectorField rhoUf
         IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
-    linearInterpolate(rho*U)
+    fvc::interpolate(rho*U)
 );
 
+{
+    surfaceVectorField n(mesh.Sf()/mesh.magSf());
+    rhoUf += n*(phi/mesh.magSf() - (n & rhoUf));
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/finiteVolume/cfdTools/incompressible/createUf.H

@@ -46,9 +46,15 @@ surfaceVectorField Uf
         IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
-    linearInterpolate(U)
+    fvc::interpolate(U)
 );
 
+{
+    surfaceVectorField n(mesh.Sf()/mesh.magSf());
+    Uf += n*(phi/mesh.magSf() - (n & Uf));
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif