From a8c917fd4b69e0cc8e0295ce921cbacfa630365d Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Wed, 11 Dec 2013 17:26:34 +0000
Subject: [PATCH] DyM solvers: correct Uf using phi after construction

---
 applications/solvers/combustion/engineFoam/pEqn.H         | 5 +----
 .../compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H    | 5 +----
 .../solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H    | 5 +----
 .../incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H        | 2 +-
 .../solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H   | 5 +----
 .../multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H    | 2 +-
 .../compressibleInterFoam/compressibleInterDyMFoam/pEqn.H | 2 +-
 .../multiphase/interFoam/interDyMFoam/interDyMFoam.C      | 2 +-
 .../solvers/multiphase/interFoam/interDyMFoam/pEqn.H      | 2 +-
 .../interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C     | 2 +-
 .../interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H   | 2 +-
 src/finiteVolume/cfdTools/compressible/createRhoUf.H      | 7 ++++++-
 src/finiteVolume/cfdTools/incompressible/createUf.H       | 8 +++++++-
 13 files changed, 24 insertions(+), 25 deletions(-)

diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
index 2f6333cc976..bd12a82495b 100644
--- a/applications/solvers/combustion/engineFoam/pEqn.H
+++ b/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())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index fd892bb9327..1c58a02c339 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/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())
diff --git a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
index 65afd1aeb88..4799f645672 100644
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/pEqn.H
+++ b/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));
 }
diff --git a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H b/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
index 7ecfe1f6623..e2938bdb949 100644
--- a/applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H
+++ b/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
diff --git a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H b/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
index b10c24aa0b7..d23ec0f9505 100644
--- a/applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/pEqn.H
+++ b/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())
diff --git a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
index 0f0d5768414..a0e0e6c74dc 100644
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/pEqn.H
+++ b/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));
     }
 }
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index c86921371f5..a014a5ab4ed 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/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
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
index 0a3d7ed11e9..ba7b4efa9ae 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C
+++ b/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"
 
diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
index c30b361e72f..76441207919 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/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
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
index 0bfd8ad56ea..810edeff308 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/interPhaseChangeDyMFoam.C
+++ b/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"
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
index 7863f5d7485..6fac1df95c2 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/pEqn.H
+++ b/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
diff --git a/src/finiteVolume/cfdTools/compressible/createRhoUf.H b/src/finiteVolume/cfdTools/compressible/createRhoUf.H
index e91115ae75b..312aec60dca 100644
--- a/src/finiteVolume/cfdTools/compressible/createRhoUf.H
+++ b/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
diff --git a/src/finiteVolume/cfdTools/incompressible/createUf.H b/src/finiteVolume/cfdTools/incompressible/createUf.H
index aab92ce0e2a..0b65f2daed2 100644
--- a/src/finiteVolume/cfdTools/incompressible/createUf.H
+++ b/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
-- 
GitLab