Rewrite of ddtPhiCorr - ddtCorr and density-weight HbyA on compressible solvers.

For DyM solvers phiAbs is replaced by Uf but this conversion is currently not complete

applications/solvers/DNS/dnsFoam/dnsFoam.C

@@ -85,7 +85,8 @@ int main(int argc, char *argv[])
 
         for (int corr=1; corr<=1; corr++)
         {
-            volScalarField rAU("Dp", 1.0/UEqn.A());
+            volScalarField rAU(1.0/UEqn.A());
+            surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
             volVectorField HbyA("HbyA", U);
             HbyA = rAU*UEqn.H();
 
@@ -93,12 +94,12 @@ int main(int argc, char *argv[])
             (
                 "phiHbyA",
                 (fvc::interpolate(HbyA) & mesh.Sf())
-              + fvc::ddtPhiCorr(rAU, U, phi)
+              + Dp*fvc::ddtCorr(U, phi)
             );
 
             fvScalarMatrix pEqn
             (
-                fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
+                fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
             );
 
             pEqn.solve();

applications/solvers/combustion/PDRFoam/pEqn.H

@@ -1,6 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
+
 volVectorField HbyA("HbyA", U);
 HbyA = invA & UEqn.H();
 
@@ -11,9 +12,9 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        )
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+        )/fvc::interpolate(rho)
     );
 
     while (pimple.correctNonOrthogonal())
@@ -38,10 +39,9 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)*
         (
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
         )
     );
 

applications/solvers/combustion/XiFoam/pEqn.H

@@ -1,6 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 
@@ -11,9 +13,9 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        )
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
+        )/fvc::interpolate(rho)
     );
 
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -24,7 +26,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(Dp, p)
          ==
             fvOptions(psi, p, rho.name())
         );
@@ -44,10 +46,9 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        (
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
         )
     );
 
@@ -59,7 +60,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(Dp, p)
          ==
             fvOptions(psi, p, rho.name())
         );

applications/solvers/combustion/engineFoam/pEqn.H

@@ -1,6 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 
@@ -10,7 +12,13 @@ if (pimple.transonic())
     (
         "phid",
         fvc::interpolate(psi)
-       *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
+       *(
+            (
+                (fvc::interpolate(rho*HbyA) & mesh.Sf())
+              //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+            )/fvc::interpolate(rho)
+          - fvc::meshPhi(rho, U)
+        )
     );
 
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -41,8 +49,11 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
+        (
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+        )
+      - fvc::interpolate(rho)*fvc::meshPhi(rho, U)
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);

applications/solvers/combustion/fireFoam/pEqn.H

@@ -1,22 +1,22 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 phi.boundaryField() =
     fvc::interpolate(rho.boundaryField()*U.boundaryField())
   & mesh.Sf().boundaryField();
 
-surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
 
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
-    fvc::interpolate(rho)
-   *(
-        (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+    (
+        (fvc::interpolate(rho*HbyA) & mesh.Sf())
+      + Dp*fvc::ddtCorr(rho, U, phi)
     )
   + phig
 );
@@ -30,7 +30,7 @@ while (pimple.correctNonOrthogonal())
         fvc::ddt(psi, rho)*gh
       + fvc::div(phiHbyA)
       + fvm::ddt(psi, p_rgh)
-      - fvm::laplacian(rhorAUf, p_rgh)
+      - fvm::laplacian(Dp, p_rgh)
      ==
         parcels.Srho()
       + surfaceFilm.Srho()
@@ -44,7 +44,7 @@ while (pimple.correctNonOrthogonal())
     if (pimple.finalNonOrthogonalIter())
     {
         phi = phiHbyA + p_rghEqn.flux();
-        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/Dp);
         U.correctBoundaryConditions();
         fvOptions.correct(U);
     }

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -1,6 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 
@@ -11,9 +13,9 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        )
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
+        )/fvc::interpolate(rho)
     );
 
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -44,10 +46,9 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        (
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
         )
     );
 

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H

@@ -6,20 +6,19 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+    surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        (
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
         )
       + phig
     );
@@ -39,7 +38,7 @@
         fvScalarMatrix p_rghEqn
         (
             p_rghDDtEqn
-          - fvm::laplacian(rhorAUf, p_rgh)
+          - fvm::laplacian(Dp, p_rgh)
         );
 
         fvOptions.constrain(p_rghEqn);
@@ -56,7 +55,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
             K = 0.5*magSqr(U);

applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H

@@ -6,6 +6,8 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
+    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
@@ -14,8 +16,10 @@
         surfaceScalarField phiHbyA
         (
             "phiHbyA",
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+            (
+                (fvc::interpolate(rho*HbyA) & mesh.Sf())
+              + Dp*fvc::ddtCorr(rho, U, phi)
+            )/fvc::interpolate(rho)
         );
 
         fvOptions.makeRelative(phiHbyA);
@@ -35,7 +39,7 @@
             fvScalarMatrix pEqn
             (
                 pDDtEqn
-              - fvm::laplacian(rho*rAU, p)
+              - fvm::laplacian(Dp, p)
              ==
                 fvOptions(psi, p, rho.name())
             );
@@ -55,10 +59,9 @@
         surfaceScalarField phiHbyA
         (
             "phiHbyA",
-            fvc::interpolate(rho)
-           *(
-                (fvc::interpolate(HbyA) & mesh.Sf())
-              + fvc::ddtPhiCorr(rAU, rho, U, phi)
+            (
+                (fvc::interpolate(rho*HbyA) & mesh.Sf())
+              + Dp*fvc::ddtCorr(rho, U, phi)
             )
         );
 
@@ -77,7 +80,7 @@
             fvScalarMatrix pEqn
             (
                 pDDtEqn
-              - fvm::laplacian(rho*rAU, p)
+              - fvm::laplacian(Dp, p)
             );
 
             fvOptions.constrain(pEqn);

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -4,6 +4,8 @@ rho = min(rho, rhoMax);
 rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
 
@@ -19,15 +21,13 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        )
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
+        )/fvc::interpolate(rho)
     );
 
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
 
-    volScalarField Dp("Dp", rho*rAU);
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -54,10 +54,9 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        (
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
         )
     );
 

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

@@ -4,6 +4,8 @@ rho = min(rho, rhoMax);
 rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
 
@@ -19,15 +21,13 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
-        )
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phiAbs)
+        )/fvc::interpolate(rho)
     );
 
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
 
-    volScalarField Dp("Dp", rho*rAU);
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -54,18 +54,12 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          - fvc::meshPhi(rho, U)
-          + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
-        )
+        (fvc::interpolate(rho*HbyA) & mesh.Sf())
+      + Dp*fvc::ddtCorr(rho, U, phiAbs)
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
-    volScalarField Dp("Dp", rho*rAU);
-
     while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C

@@ -59,7 +59,7 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
-    // Create old-time absolute flux for ddtPhiCorr
+    // Create old-time absolute flux for ddtCorr
     surfaceScalarField phiAbs("phiAbs", phi);
 
 
@@ -75,7 +75,7 @@ int main(int argc, char *argv[])
         // Make the fluxes absolute before mesh-motion
         fvc::makeAbsolute(phi, rho, U);
 
-        // Update absolute flux for ddtPhiCorr
+        // Update absolute flux for ddtCorr
         phiAbs = phi;
 
         #include "setDeltaT.H"

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

@@ -20,9 +20,9 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        )
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+        )/fvc::interpolate(rho)
     );
 
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -64,10 +64,9 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        (
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
         )
     );
 

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -12,7 +12,9 @@
         surfaceScalarField phid
         (
             "phid",
-            fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::interpolate(psi)
+           *(fvc::interpolate(rho*HbyA) & mesh.Sf())
+           /fvc::interpolate(rho)
         );
 
         fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -47,7 +49,7 @@
         surfaceScalarField phiHbyA
         (
             "phiHbyA",
-            fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+            fvc::interpolate(rho*HbyA) & mesh.Sf()
         );
 
         fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H

@@ -13,7 +13,9 @@ if (simple.transonic())
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::interpolate(psi)
+       *(fvc::interpolate(rho*HbyA) & mesh.Sf())
+       /fvc::interpolate(rho)
     );
 
     surfaceScalarField phic
@@ -57,7 +59,7 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+        fvc::interpolate(rho*HbyA) & mesh.Sf()
     );
 
     closedVolume = adjustPhi(phiHbyA, U, p);

applications/solvers/compressible/sonicFoam/pEqn.H

@@ -1,22 +1,21 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 
 surfaceScalarField phid
 (
     "phid",
-    fvc::interpolate(psi)
-   *(
-        (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
-    )
+    fvc::interpolate(psi)*
+    (
+        (mesh.Sf() & fvc::interpolate(rho*HbyA))
+      + Dp*fvc::ddtCorr(rho, U, phi)
+    )/fvc::interpolate(rho)
 );
 
-
-volScalarField Dp("Dp", rho*rAU);
-
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())
 {

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

@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    createUf
+
+Description
+    Creates and initialises the velocity velocity field Uf.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef createUf_H
+#define createUf_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Info<< "Reading/calculating face velocity rhoUf\n" << endl;
+
+surfaceVectorField rhoUf
+(
+    IOobject
+    (
+        "rhoUf",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(rho*U)
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -1,6 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 
@@ -9,13 +11,14 @@ surfaceScalarField phid
     "phid",
     fvc::interpolate(psi)
    *(
-        (fvc::interpolate(HbyA) & mesh.Sf())
+        (
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+        )/fvc::interpolate(rho)
       - fvc::meshPhi(rho, U)
     )
 );
 
-volScalarField Dp("Dp", rho*rAU);
-
 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
     fvScalarMatrix pEqn

applications/solvers/compressible/sonicFoam/sonicDyMFoam/sonicDyMFoam.C

@@ -44,6 +44,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
+    #include "createRhoUf.H"
     #include "initContinuityErrs.H"
 
     pimpleControl pimple(mesh);
@@ -63,6 +64,12 @@ 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

applications/solvers/compressible/sonicFoam/sonicLiquidFoam/sonicLiquidFoam.C

@@ -76,6 +76,8 @@ int main(int argc, char *argv[])
             while (pimple.correct())
             {
                 volScalarField rAU(1.0/UEqn.A());
+                surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
                 U = rAU*UEqn.H();
 
                 surfaceScalarField phid
@@ -83,13 +85,12 @@ int main(int argc, char *argv[])
                     "phid",
                     psi
                    *(
-                       (fvc::interpolate(U) & mesh.Sf())
-                     + fvc::ddtPhiCorr(rAU, rho, U, phi)
-                   )
+                       (fvc::interpolate(rho*U) & mesh.Sf())
+                     + Dp*fvc::ddtCorr(rho, U, phi)
+                    )/fvc::interpolate(rho)
                 );
 
                 phi = (rhoO/psi)*phid;
-                volScalarField Dp("Dp", rho*rAU);
 
                 fvScalarMatrix pEqn
                 (

applications/solvers/electromagnetics/mhdFoam/mhdFoam.C

@@ -93,6 +93,8 @@ int main(int argc, char *argv[])
             for (int corr=0; corr<nCorr; corr++)
             {
                 volScalarField rAU(1.0/UEqn.A());
+                surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+
                 volVectorField HbyA("HbyA", U);
                 HbyA = rAU*UEqn.H();
 
@@ -100,14 +102,14 @@ int main(int argc, char *argv[])
                 (
                     "phiHbyA",
                     (fvc::interpolate(HbyA) & mesh.Sf())
-                  + fvc::ddtPhiCorr(rAU, U, phi)
+                  + Dp*fvc::ddtCorr(U, phi)
                 );
 
                 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
                 {
                     fvScalarMatrix pEqn
                     (
-                        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
+                        fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
                     );
 
                     pEqn.setReference(pRefCell, pRefValue);
@@ -140,14 +142,15 @@ int main(int argc, char *argv[])
 
             BEqn.solve();
 
-            volScalarField rBA(1.0/BEqn.A());
+            volScalarField rAB(1.0/BEqn.A());
+            surfaceScalarField DpB("DpB", fvc::interpolate(rAB));
 
             phiB = (fvc::interpolate(B) & mesh.Sf())
-                + fvc::ddtPhiCorr(rBA, B, phiB);
+                + DpB*fvc::ddtCorr(B, phiB);
 
             fvScalarMatrix pBEqn
             (
-                fvm::laplacian(rBA, pB) == fvc::div(phiB)
+                fvm::laplacian(DpB, pB) == fvc::div(phiB)
             );
             pBEqn.solve();
 

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H

@@ -1,17 +1,17 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
+    surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rhok)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, U, phi)
+      + Dp*fvc::ddtCorr(U, phi)
       + phig
     );
 
@@ -19,7 +19,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -36,7 +36,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -6,20 +6,19 @@
     thermo.rho() -= psi*p_rgh;
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+    surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        (
+            (fvc::interpolate(rho*U) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
         )
       + phig
     );
@@ -39,7 +38,7 @@
         fvScalarMatrix p_rghEqn
         (
             p_rghDDtEqn
-          - fvm::laplacian(rhorAUf, p_rgh)
+          - fvm::laplacian(Dp, p_rgh)
         );
 
         fvOptions.constrain(p_rghEqn);
@@ -56,7 +55,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
             K = 0.5*magSqr(U);

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -14,7 +14,7 @@
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+        (fvc::interpolate(rho*HbyA) & mesh.Sf())
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H

@@ -16,7 +16,7 @@
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+        (fvc::interpolate(rho*HbyA) & mesh.Sf())
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H

@@ -6,20 +6,19 @@
     rho = thermo.rho();
 
     volScalarField rAU(1.0/UEqn().A());
-    surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn().H();
 
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+    surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        (
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
         )
       + phig
     );
@@ -42,7 +41,7 @@
             fvScalarMatrix p_rghEqn
             (
                 p_rghDDtEqn
-              - fvm::laplacian(rhorAUf, p_rgh)
+              - fvm::laplacian(Dp, p_rgh)
             );
 
             p_rghEqn.solve
@@ -64,7 +63,7 @@
             {
                 phi = phiHbyA + p_rghEqn.flux();
                 U = HbyA
-                  + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+                  + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
                 U.correctBoundaryConditions();
                 fvOptions.correct(U);
                 K = 0.5*magSqr(U);

applications/solvers/incompressible/icoFoam/icoFoam.C

@@ -74,7 +74,7 @@ int main(int argc, char *argv[])
             (
                 "phiHbyA",
                 (fvc::interpolate(HbyA) & mesh.Sf())
-              + fvc::ddtPhiCorr(rAU, U, phi)
+              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
             );
 
             adjustPhi(phiHbyA, U, p);

applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C

@@ -78,7 +78,7 @@ int main(int argc, char *argv[])
             (
                 "phiHbyA",
                 (fvc::interpolate(HbyA) & mesh.Sf())
-              + fvc::ddtPhiCorr(rAU, U, phi)
+              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
             );
 
             adjustPhi(phiHbyA, U, p);

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

@@ -11,7 +11,7 @@ surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (fvc::interpolate(HbyA) & mesh.Sf())
-  + fvc::ddtPhiCorr(rAUrel, Urel, phi)
+  + fvc::interpolate(rAUrel)*fvc::ddtCorr(Urel, phi)
 );
 
 adjustPhi(phiHbyA, Urel, p);

applications/solvers/incompressible/pimpleFoam/pEqn.H

@@ -1,3 +1,5 @@
+surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
 
@@ -10,7 +12,7 @@ surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (fvc::interpolate(HbyA) & mesh.Sf())
-  + fvc::ddtPhiCorr(rAU, U, phi)
+  + Dp*fvc::ddtCorr(U, phi)
 );
 
 fvOptions.makeRelative(phiHbyA);
@@ -23,7 +25,7 @@ while (pimple.correctNonOrthogonal())
     // Pressure corrector
     fvScalarMatrix pEqn
     (
-        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
+        fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
     );
 
     pEqn.setReference(pRefCell, pRefValue);

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

@@ -1,6 +1,5 @@
+if (mesh.changing())
 {
-    if (mesh.changing())
-    {
     forAll(U.boundaryField(), patchI)
     {
         if (U.boundaryField()[patchI].fixesValue())
@@ -20,8 +19,9 @@
               & mesh.Sf().boundaryField()[patchI];
         }
     }
-    }
+}
 
+{
     volScalarField pcorr
     (
         IOobject

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

@@ -1,42 +0,0 @@
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    #include "createPhi.H"
-
-
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
-
-    singlePhaseTransportModel laminarTransport(U, phi);
-
-    autoPtr<incompressible::turbulenceModel> turbulence
-    (
-        incompressible::turbulenceModel::New(U, phi, laminarTransport)
-    );

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

@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    createUf
+
+Description
+    Creates and initialises the velocity velocity field Uf.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef createUf_H
+#define createUf_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Info<< "Reading/calculating face velocity Uf\n" << endl;
+
+surfaceVectorField Uf
+(
+    IOobject
+    (
+        "Uf",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(U)
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -1,3 +1,5 @@
+surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
 
@@ -10,7 +12,7 @@ surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (fvc::interpolate(HbyA) & mesh.Sf())
-  + fvc::ddtPhiCorr(rAU, U, phiAbs)
+  + Dp*fvc::ddtCorr(U, Uf)
 );
 
 if (p.needReference())
@@ -24,7 +26,7 @@ while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix pEqn
     (
-        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
+        fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
     );
 
     pEqn.setReference(pRefCell, pRefValue);
@@ -42,9 +44,15 @@ while (pimple.correctNonOrthogonal())
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-// Make the fluxes relative to the mesh motion
-fvc::makeRelative(phi, U);
-
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
+
+{
+    Uf = fvc::interpolate(U);
+    surfaceVectorField n(mesh.Sf()/mesh.magSf());
+    Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+}
+
+// Make the fluxes relative to the mesh motion
+fvc::makeRelative(phi, U);

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C

@@ -51,15 +51,13 @@ int main(int argc, char *argv[])
     pimpleControl pimple(mesh);
 
     #include "createFields.H"
+    #include "createUf.H"
     #include "createFvOptions.H"
     #include "readTimeControls.H"
     #include "createPcorrTypes.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
-    // Create old-time absolute flux for ddtPhiCorr
-    surfaceScalarField phiAbs("phiAbs", phi);
-
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -69,12 +67,6 @@ int main(int argc, char *argv[])
         #include "readControls.H"
         #include "CourantNo.H"
 
-        // Make the fluxes absolute
-        fvc::makeAbsolute(phi, U);
-
-        // Update absolute flux for ddtPhiCorr
-        phiAbs = phi;
-
         #include "setDeltaT.H"
 
         runTime++;
@@ -83,12 +75,15 @@ int main(int argc, char *argv[])
 
         mesh.update();
 
+        // Calculate absolute flux from the mapped surface velocity
+        phi = mesh.Sf() & Uf;
+
         if (mesh.changing() && correctPhi)
         {
             #include "correctPhi.H"
         }
 
-        // Make the fluxes relative to the mesh motion
+        // Make the flux relative to the mesh motion
         fvc::makeRelative(phi, U);
 
         if (mesh.changing() && checkMeshCourantNo)

applications/solvers/incompressible/pisoFoam/pisoFoam.C

@@ -87,7 +87,7 @@ int main(int argc, char *argv[])
                 (
                     "phiHbyA",
                     (fvc::interpolate(HbyA) & mesh.Sf())
-                  + fvc::ddtPhiCorr(rAU, U, phi)
+                  + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
                 );
 
                 adjustPhi(phiHbyA, U, p);

applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C

@@ -111,7 +111,7 @@ int main(int argc, char *argv[])
                 (
                     "phiHbyA",
                     (fvc::interpolate(HbyA) & mesh.Sf())
-                  + fvc::ddtPhiCorr(rAU, h, hU, phi)
+                  + fvc::interpolate(rAU)*fvc::ddtCorr(h, hU, phi)
                   - phih0
                 );
 

applications/solvers/lagrangian/coalChemistryFoam/pEqn.H

@@ -1,6 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 
@@ -11,9 +13,9 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        )
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
+        )/fvc::interpolate(rho)
     );
 
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -24,7 +26,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(Dp, p)
          ==
             coalParcels.Srho()
           + fvOptions(psi, p, rho.name())
@@ -45,10 +47,9 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        (
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
         )
     );
 
@@ -60,7 +61,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(Dp, p)
          ==
             coalParcels.Srho()
           + fvOptions(psi, p, rho.name())

applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H

@@ -1,19 +1,19 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 
-surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
 
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
-    fvc::interpolate(rho)
-   *(
-        (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+    (
+        (fvc::interpolate(rho*HbyA) & mesh.Sf())
+      + Dp*fvc::ddtCorr(rho, U, phi)
     )
   + phig
 );
@@ -27,7 +27,7 @@ while (pimple.correctNonOrthogonal())
         fvc::ddt(psi, rho)*gh
       + fvc::div(phiHbyA)
       + fvm::ddt(psi, p_rgh)
-      - fvm::laplacian(rhorAUf, p_rgh)
+      - fvm::laplacian(Dp, p_rgh)
      ==
         parcels.Srho()
       + surfaceFilm.Srho()
@@ -41,7 +41,7 @@ while (pimple.correctNonOrthogonal())
     if (pimple.finalNonOrthogonalIter())
     {
         phi = phiHbyA + p_rghEqn.flux();
-        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/Dp);
         U.correctBoundaryConditions();
         fvOptions.correct(U);
     }

applications/solvers/lagrangian/reactingParcelFoam/LTSReactingParcelFoam/LTSReactingParcelFoam.C

@@ -30,9 +30,6 @@ Description
     parcels and porous media, including run-time selectable finitite volume
     options, e.g. sources, constraints
 
-    Note: ddtPhiCorr not used here when porous zones are active
-    - not well defined for porous calculations
-
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"

applications/solvers/lagrangian/reactingParcelFoam/pEqn.H

@@ -6,16 +6,17 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
+    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        (
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
         )
     );
 
@@ -35,7 +36,7 @@
         fvScalarMatrix pEqn
         (
             pDDtEqn
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(Dp, p)
         );
 
         fvOptions.constrain(pEqn);

applications/solvers/lagrangian/sprayFoam/pEqn.H

@@ -1,6 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 
@@ -11,9 +13,9 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        )
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
+        )/fvc::interpolate(rho)
     );
 
     fvOptions.makeRelative(fvc::interpolate(psi), phid);
@@ -24,7 +26,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(Dp, p)
          ==
             parcels.Srho()
           + fvOptions(psi, p, rho.name())
@@ -45,10 +47,9 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-            (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        (
+            (fvc::interpolate(rho*HbyA) & mesh.Sf())
+          + Dp*fvc::ddtCorr(rho, U, phi)
         )
     );
 
@@ -60,7 +61,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(Dp, p)
          ==
             parcels.Srho()
           + fvOptions(psi, p, rho.name())

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

@@ -1,6 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
 
@@ -11,8 +13,11 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            ((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+            (
+                (fvc::interpolate(rho*HbyA) & mesh.Sf())
+              //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+            )/fvc::interpolate(rho)
+          - fvc::meshPhi(rho, U)
         )
     );
 
@@ -24,7 +29,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(Dp, p)
          ==
             parcels.Srho()
           + fvOptions(psi, p, rho.name())
@@ -45,11 +50,11 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-            ((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        (
+            (fvc::interpolate(HbyA) & mesh.Sf())
+          //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
         )
+      - fvc::interpolate(rho)*fvc::meshPhi(rho, U)
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
@@ -60,7 +65,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(rho*rAU, p)
+          - fvm::laplacian(Dp, p)
          ==
             parcels.Srho()
           + fvOptions(psi, p, rho.name())

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

@@ -12,16 +12,16 @@
     surfaceScalarField rhof("rhof", fvc::interpolate(rho));
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf("Dp", rhof*fvc::interpolate(rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
     phiv = (fvc::interpolate(HbyA) & mesh.Sf())
-         + fvc::ddtPhiCorr(rAU, rho, U, phivAbs);
+         + Dp*fvc::ddtCorr(U, phivAbs);
     fvc::makeRelative(phiv, U);
 
-    surfaceScalarField phiGradp(rAUf*mesh.magSf()*fvc::snGrad(p));
+    surfaceScalarField phiGradp(Dp*mesh.magSf()*fvc::snGrad(p));
 
     phiv -= phiGradp/rhof;
 
@@ -35,7 +35,7 @@
           + psi*correction(fvm::ddt(p))
           + fvc::div(phiv, rho)
           + fvc::div(phiGradp)
-          - fvm::laplacian(rAUf, p)
+          - fvm::laplacian(Dp, p)
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

applications/solvers/multiphase/cavitatingFoam/pEqn.H

@@ -12,15 +12,15 @@
     surfaceScalarField rhof("rhof", fvc::interpolate(rho));
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf("Dp", rhof*fvc::interpolate(rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
     phiv = (fvc::interpolate(HbyA) & mesh.Sf())
-         + fvc::ddtPhiCorr(rAU, rho, U, phiv);
+         + Dp*fvc::ddtCorr(U, phiv);
 
-    surfaceScalarField phiGradp(rAUf*mesh.magSf()*fvc::snGrad(p));
+    surfaceScalarField phiGradp(Dp*mesh.magSf()*fvc::snGrad(p));
 
     phiv -= phiGradp/rhof;
 
@@ -32,7 +32,7 @@
           - (rhol0 + (psil - psiv)*pSat)*fvc::ddt(alphav) - pSat*fvc::ddt(psi)
           + fvc::div(phiv, rho)
           + fvc::div(phiGradp)
-          - fvm::laplacian(rAUf, p)
+          - fvm::laplacian(Dp, p)
         );
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
     );
     phi = phiHbyA;
 
@@ -18,7 +18,7 @@
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf()
+        )*Dp*mesh.magSf()
     );
 
     phiHbyA += phig;
@@ -70,7 +70,7 @@
         fvScalarMatrix p_rghEqnIncomp
         (
             fvc::div(phiHbyA)
-          - fvm::laplacian(rAUf, p_rgh)
+          - fvm::laplacian(Dp, p_rgh)
         );
 
         solve
@@ -97,7 +97,7 @@
             phi = phiHbyA + p_rghEqnIncomp.flux();
 
             U = HbyA
-              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/Dp);
             U.correctBoundaryConditions();
         }
     }

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

@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
     );
     adjustPhi(phiHbyA, U, p_rgh);
     mrfZones.makeRelative(phiHbyA);
@@ -20,7 +20,7 @@
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf()
+        )*Dp*mesh.magSf()
     );
 
     phiHbyA += phig;
@@ -29,7 +29,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -40,7 +40,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }

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

@@ -1,20 +1,27 @@
+if (mesh.changing())
 {
-    #include "continuityErrs.H"
-
-    wordList pcorrTypes
-    (
-        p_rgh.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    forAll (p_rgh.boundaryField(), i)
+    forAll(U.boundaryField(), patchI)
     {
-        if (p_rgh.boundaryField()[i].fixesValue())
+        if (U.boundaryField()[patchI].fixesValue())
         {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+            U.boundaryField()[patchI].initEvaluate();
         }
     }
 
+    forAll(U.boundaryField(), patchI)
+    {
+        if (U.boundaryField()[patchI].fixesValue())
+        {
+            U.boundaryField()[patchI].evaluate();
+
+            phi.boundaryField()[patchI] =
+                U.boundaryField()[patchI]
+              & mesh.Sf().boundaryField()[patchI];
+        }
+    }
+}
+
+{
     volScalarField pcorr
     (
         IOobject
@@ -30,17 +37,13 @@
         pcorrTypes
     );
 
-    dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
-
-    adjustPhi(phi, U, pcorr);
-
-    fvc::makeAbsolute(phi, U);
+    dimensionedScalar Dp("Dp", dimTime/rho.dimensions(), 1.0);
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
+            fvm::laplacian(Dp, pcorr) == fvc::div(phi)
         );
 
         pcorrEqn.setReference(pRefCell, pRefValue);
@@ -49,9 +52,6 @@
         if (pimple.finalNonOrthogonalIter())
         {
             phi -= pcorrEqn.flux();
-            phiAbs = phi;
-            phiAbs.oldTime() = phi;
-            fvc::makeRelative(phi, U);
         }
     }
 }

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

@@ -0,0 +1,13 @@
+    wordList pcorrTypes
+    (
+        p_rgh.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+    for (label i=0; i<p_rgh.boundaryField().size(); i++)
+    {
+        if (p_rgh.boundaryField()[i].fixesValue())
+        {
+            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+        }
+    }

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

@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    createUf
+
+Description
+    Creates and initialises the velocity velocity field Uf.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef createUf_H
+#define createUf_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Info<< "Reading/calculating face velocity Uf\n" << endl;
+
+surfaceVectorField Uf
+(
+    IOobject
+    (
+        "Uf",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(U)
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -50,15 +50,13 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
     #include "initContinuityErrs.H"
-    #include "createFields.H"
 
     pimpleControl pimple(mesh);
 
+    #include "createFields.H"
+    #include "createUf.H"
     #include "readTimeControls.H"
-
-    surfaceScalarField phiAbs("phiAbs", phi);
-    fvc::makeAbsolute(phiAbs, U);
-
+    #include "createPcorrTypes.H"
     #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
@@ -80,21 +78,7 @@ int main(int argc, char *argv[])
 
         scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
 
-        {
-            // Ensure old-time U exists for mapping
-            U.oldTime();
-
-            // Calculate the relative velocity used to map the relative flux phi
-            volVectorField Urel("Urel", U);
-
-            if (mesh.moving())
-            {
-                Urel -= fvc::reconstruct(fvc::meshPhi(U));
-            }
-
-            // Do any mesh changes
         mesh.update();
-        }
 
         if (mesh.changing())
         {
@@ -108,7 +92,13 @@ int main(int argc, char *argv[])
 
         if (mesh.changing() && correctPhi)
         {
+            // Calculate absolute flux from the mapped surface velocity
+            phi = mesh.Sf() & Uf;
+
             #include "correctPhi.H"
+
+            // Make the flux relative to the mesh motion
+            fvc::makeRelative(phi, U);
         }
 
         if (mesh.changing() && checkMeshCourantNo)

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

@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
     );
 
     if (p_rgh.needReference())
@@ -19,14 +19,14 @@
         fvc::makeAbsolute(phiHbyA, U);
     }
 
-    phiAbs = phiHbyA;
+    surfaceScalarField phiAbs("phiAbs", phiHbyA);
 
     surfaceScalarField phig
     (
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf()
+        )*Dp*mesh.magSf()
     );
 
     phiHbyA += phig;
@@ -35,7 +35,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -46,7 +46,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
@@ -54,7 +54,11 @@
 
     #include "continuityErrs.H"
 
-    phiAbs = phi;
+    {
+        Uf = fvc::interpolate(U);
+        surfaceVectorField n(mesh.Sf()/mesh.magSf());
+        Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+    }
 
     // Make the fluxes relative to the mesh motion
     fvc::makeRelative(phi, U);

applications/solvers/multiphase/interFoam/pEqn.H

@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
     );
 
     adjustPhi(phiHbyA, U, p_rgh);
@@ -20,7 +20,7 @@
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf()
+        )*Dp*mesh.magSf()
     );
 
     phiHbyA += phig;
@@ -29,7 +29,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -40,7 +40,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }

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

@@ -60,14 +60,13 @@ int main(int argc, char *argv[])
     #include "createDynamicFvMesh.H"
     #include "readGravitationalAcceleration.H"
     #include "initContinuityErrs.H"
-    #include "createFields.H"
-    #include "readTimeControls.H"
 
     pimpleControl pimple(mesh);
 
-    surfaceScalarField phiAbs("phiAbs", phi);
-    fvc::makeAbsolute(phiAbs, U);
-
+    #include "createFields.H"
+    #include "../interFoam/interDyMFoam/createUf.H"
+    #include "readTimeControls.H"
+    #include "../interFoam/interDyMFoam/createPcorrTypes.H"
     #include "../interFoam/interDyMFoam/correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
@@ -88,21 +87,7 @@ int main(int argc, char *argv[])
 
         scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
 
-        {
-            // Ensure old-time U exists for mapping
-            U.oldTime();
-
-            // Calculate the relative velocity used to map the relative flux phi
-            volVectorField Urel("Urel", U);
-
-            if (mesh.moving())
-            {
-                Urel -= fvc::reconstruct(fvc::meshPhi(U));
-            }
-
-            // Do any mesh changes
         mesh.update();
-        }
 
         if (mesh.changing())
         {
@@ -116,7 +101,13 @@ int main(int argc, char *argv[])
 
         if (mesh.changing() && correctPhi)
         {
+            // Calculate absolute flux from the mapped surface velocity
+            phi = mesh.Sf() & Uf;
+
             #include "../interFoam/interDyMFoam/correctPhi.H"
+
+            // Make the flux relative to the mesh motion
+            fvc::makeRelative(phi, U);
         }
 
         if (mesh.changing() && checkMeshCourantNo)

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

@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, Uf)
     );
 
     if (p_rgh.needReference())
@@ -19,14 +19,14 @@
         fvc::makeAbsolute(phiHbyA, U);
     }
 
-    phiAbs = phiHbyA;
+    surfaceScalarField phiAbs("phiAbs", phiHbyA);
 
     surfaceScalarField phig
     (
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf()
+        )*Dp*mesh.magSf()
     );
 
     phiHbyA += phig;
@@ -39,7 +39,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
+            fvc::div(phiHbyA) - fvm::laplacian(Dp, p_rgh)
           - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
         );
 
@@ -51,13 +51,17 @@
         {
             phi = phiHbyA + p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
     }
 
-    phiAbs = phi;
+    {
+        Uf = fvc::interpolate(U);
+        surfaceVectorField n(mesh.Sf()/mesh.magSf());
+        Uf += mesh.Sf()*(phi - (mesh.Sf() & Uf))/sqr(mesh.magSf());
+    }
 
     // Make the fluxes relative to the mesh motion
     fvc::makeRelative(phi, U);

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
     );
     adjustPhi(phiHbyA, U, p_rgh);
     phi = phiHbyA;
@@ -19,7 +19,7 @@
         (
             fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf()
+        )*Dp*mesh.magSf()
     );
 
     phiHbyA += phig;
@@ -32,7 +32,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
+            fvc::div(phiHbyA) - fvm::laplacian(Dp, p_rgh)
           - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
         );
 
@@ -44,7 +44,7 @@
         {
             phi = phiHbyA + p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }

applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H

@@ -94,7 +94,7 @@
         phiHbyAs[phasei] =
         (
             (fvc::interpolate(HbyAs[phasei]) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAUs[phasei], alpha, phase.U(), phase.phi())
+          + rAlphaAUfs[phasei]*fvc::ddtCorr(phase.U(), phase.phi())
         );
         mrfZones.makeRelative(phiHbyAs[phasei]);
 

applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H

@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
     );
     adjustPhi(phiHbyA, U, p_rgh);
     mrfZones.makeRelative(phiHbyA);
@@ -17,7 +17,7 @@
 
     surfaceScalarField phig
     (
-        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+        - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
     );
 
     phiHbyA += phig;
@@ -26,7 +26,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -37,7 +37,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
             U.correctBoundaryConditions();
         }
     }

applications/solvers/multiphase/multiphaseInterFoam/pEqn.H

@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -9,7 +9,7 @@
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
     );
     adjustPhi(phiHbyA, U, p_rgh);
     phi = phiHbyA;
@@ -19,7 +19,7 @@
         (
             mixture.surfaceTensionForce()
           - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf()
+        )*Dp*mesh.magSf()
     );
 
     phiHbyA += phig;
@@ -28,7 +28,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -39,7 +39,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
             U.correctBoundaryConditions();
         }
     }

applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H

@@ -13,7 +13,7 @@ surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (fvc::interpolate(HbyA) & mesh.Sf())
-  + fvc::ddtPhiCorr(rAU, U, phi)
+  + rAUf*fvc::ddtCorr(U, phi)
 );
 
 adjustPhi(phiHbyA, U, p_gh);

applications/solvers/multiphase/settlingFoam/pEqn.H

@@ -1,7 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-
-    surfaceScalarField rAUf("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -9,17 +8,16 @@
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::interpolate(rho)
-       *(
-           (fvc::interpolate(HbyA) & mesh.Sf())
-         + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        (
+           (fvc::interpolate(rho*HbyA) & mesh.Sf())
+         + Dp*fvc::ddtCorr(rho, U, phi)
         )
     );
     phi = phiHbyA;
 
     surfaceScalarField phig
     (
-        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+        - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
     );
 
     phiHbyA += phig;
@@ -28,7 +26,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
+            fvm::laplacian(Dp, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -39,7 +37,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
             U.correctBoundaryConditions();
         }
     }

applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H

@@ -1,6 +1,6 @@
 {
     volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
+    surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -9,14 +9,14 @@
     (
         "phiHbyA",
         (fvc::interpolate(HbyA) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
     );
     adjustPhi(phiHbyA, U, p_rgh);
     phi = phiHbyA;
 
     surfaceScalarField phig
     (
-        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+        - ghf*fvc::snGrad(rho)*Dp*mesh.magSf()
     );
 
     phiHbyA += phig;
@@ -25,7 +25,7 @@
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+            fvm::laplacian(Dp, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -36,7 +36,7 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
-            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/Dp);
             U.correctBoundaryConditions();
         }
     }

applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H

@@ -47,14 +47,14 @@
     (
         IOobject::groupName("phiHbyA", phase1.name()),
         (fvc::interpolate(HbyA1) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
+      + rAlphaAU1f*fvc::ddtCorr(U1, phi1)
     );
 
     surfaceScalarField phiHbyA2
     (
         IOobject::groupName("phiHbyA", phase2.name()),
         (fvc::interpolate(HbyA2) & mesh.Sf())
-      + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2)
+      + rAlphaAU2f*fvc::ddtCorr(U2, phi2)
     );
 
     phiHbyA1 +=

applications/test/LduMatrix/LduMatrixTest3.C

@@ -105,7 +105,7 @@ int main(int argc, char *argv[])
 
             U = rAU*UEqn.H();
             phi = (fvc::interpolate(U) & mesh.Sf())
-                + fvc::ddtPhiCorr(rAU, U, phi);
+                + fvc::ddtCorr(rAU, U, phi);
 
             adjustPhi(phi, U, p);
 

applications/test/PisoFoam/PisoFoam.C

@@ -87,7 +87,7 @@ int main(int argc, char *argv[])
                 (
                     "phiHbyA",
                     (fvc::interpolate(HbyA) & mesh.Sf())
-                  + fvc::ddtPhiCorr(rAU, U, phi)
+                  + fvc::ddtCorr(rAU, U, phi)
                 );
 
                 adjustPhi(phiHbyA, U, p);

applications/test/RhoPimpleFoam/pEqn.H

@@ -20,7 +20,7 @@ if (pimple.transonic())
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+          + fvc::ddtCorr(rAU, rho, U, phi)
         )
     );
 
@@ -57,7 +57,7 @@ else
         fvc::interpolate(rho)
        *(
             (fvc::interpolate(HbyA) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+          + fvc::ddtCorr(rAU, rho, U, phi)
         )
     );
 

src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C

@@ -480,122 +480,26 @@ CoEulerDdtScheme<Type>::fvmDdt
 
 template<class Type>
 tmp<typename CoEulerDdtScheme<Type>::fluxFieldType>
-CoEulerDdtScheme<Type>::fvcDdtPhiCorr
+CoEulerDdtScheme<Type>::fvcDdtUfCorr
 (
-    const volScalarField& rA,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
 )
 {
     IOobject ddtIOobject
     (
-        "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
+        "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
         mesh().time().timeName(),
         mesh()
     );
 
-    if (mesh().moving())
-    {
-        return tmp<fluxFieldType>
-        (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rA.dimensions()*phi.dimensions()/dimTime,
-                    pTraits<typename flux<Type>::type>::zero
-                )
-            )
-        );
-    }
-    else
-    {
-        const volScalarField rDeltaT(CorDeltaT());
+    const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
 
-        return tmp<fluxFieldType>
+    fluxFieldType phiCorr
     (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
-                (
-                    fvc::interpolate(rDeltaT*rA)*phi.oldTime()
-                  - (fvc::interpolate(rDeltaT*rA*U.oldTime()) & mesh().Sf())
-                )
-            )
-        );
-    }
-}
-
-
-template<class Type>
-tmp<typename CoEulerDdtScheme<Type>::fluxFieldType>
-CoEulerDdtScheme<Type>::fvcDdtPhiCorr
-(
-    const volScalarField& rA,
-    const volScalarField& rho,
-    const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
-)
-{
-    IOobject ddtIOobject
-    (
-        "ddtPhiCorr("
-      + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
-        mesh().time().timeName(),
-        mesh()
+        mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
     );
 
-    if (mesh().moving())
-    {
-        return tmp<fluxFieldType>
-        (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
-                    pTraits<typename flux<Type>::type>::zero
-                )
-            )
-        );
-    }
-    else
-    {
-        const volScalarField rDeltaT(CorDeltaT());
-
-        if
-        (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == dimVelocity*dimArea
-        )
-        {
-            return tmp<fluxFieldType>
-            (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
-                   *(
-                        fvc::interpolate(rDeltaT*rA*rho.oldTime())*phi.oldTime()
-                      - (fvc::interpolate(rDeltaT*rA*rho.oldTime()*U.oldTime())
-                      & mesh().Sf())
-                    )
-                )
-            );
-        }
-        else if
-        (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == dimDensity*dimVelocity*dimArea
-        )
-        {
     return tmp<fluxFieldType>
     (
         new fluxFieldType
@@ -604,40 +508,199 @@ CoEulerDdtScheme<Type>::fvcDdtPhiCorr
             this->fvcDdtPhiCoeff
             (
                 U.oldTime(),
-                        phi.oldTime()/fvc::interpolate(rho.oldTime())
-                    )
-                   *(
-                        fvc::interpolate(rDeltaT*rA*rho.oldTime())
-                       *phi.oldTime()/fvc::interpolate(rho.oldTime())
-                      - (
-                            fvc::interpolate
-                            (
-                                rDeltaT*rA*rho.oldTime()*U.oldTime()
-                            ) & mesh().Sf()
-                        )
+                (mesh().Sf() & Uf.oldTime()),
+                phiCorr
             )
+           *rDeltaT*phiCorr
         )
     );
-        }
-        else if
+}
+
+
+template<class Type>
+tmp<typename CoEulerDdtScheme<Type>::fluxFieldType>
+CoEulerDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    IOobject ddtIOobject
     (
-            U.dimensions() == dimDensity*dimVelocity
-         && phi.dimensions() == dimDensity*dimVelocity*dimArea
+        "ddtCorr(" + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
+
+    fluxFieldType phiCorr
+    (
+        phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+    );
+
+    return tmp<fluxFieldType>
+    (
+        new fluxFieldType
+        (
+            ddtIOobject,
+            this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+           *rDeltaT*phiCorr
+        )
+    );
+}
+
+
+template<class Type>
+tmp<typename CoEulerDdtScheme<Type>::fluxFieldType>
+CoEulerDdtScheme<Type>::fvcDdtUfCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+)
+{
+    IOobject ddtIOobject
+    (
+        "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    const surfaceScalarField rDeltaT(fvc::interpolate(CorDeltaT()));
+
+    if
+    (
+        U.dimensions() == dimVelocity
+     && Uf.dimensions() == dimDensity*dimVelocity
     )
     {
+        GeometricField<Type, fvPatchField, volMesh> rhoU0
+        (
+            rho.oldTime()*U.oldTime()
+        );
+
+        fluxFieldType phiCorr
+        (
+            mesh().Sf() & (Uf.oldTime() - fvc::interpolate(rhoU0))
+        );
+
         return tmp<fluxFieldType>
         (
             new fluxFieldType
             (
                 ddtIOobject,
                 this->fvcDdtPhiCoeff
-                    (rho.oldTime(), U.oldTime(), phi.oldTime())
-                  * (
-                        fvc::interpolate(rDeltaT*rA)*phi.oldTime()
-                      - (
-                            fvc::interpolate(rDeltaT*rA*U.oldTime())&mesh().Sf()
+                (
+                    rhoU0,
+                    mesh().Sf() & Uf.oldTime(),
+                    phiCorr
                 )
+               *rDeltaT*phiCorr
             )
+        );
+    }
+    else if
+    (
+        U.dimensions() == dimDensity*dimVelocity
+     && Uf.dimensions() == dimDensity*dimVelocity
+    )
+    {
+        fluxFieldType phiCorr
+        (
+             mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+        );
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff
+                (
+                    U.oldTime(),
+                    mesh().Sf() & Uf.oldTime(),
+                    phiCorr
+                )
+               *rDeltaT*phiCorr
+            )
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "CoEulerDdtScheme<Type>::fvcDdtPhiCorr"
+        )   << "dimensions of Uf are not correct"
+            << abort(FatalError);
+
+        return fluxFieldType::null();
+    }
+}
+
+
+template<class Type>
+tmp<typename CoEulerDdtScheme<Type>::fluxFieldType>
+CoEulerDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    if
+    (
+        U.dimensions() == dimVelocity
+     && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        GeometricField<Type, fvPatchField, volMesh> rhoU0
+        (
+            rho.oldTime()*U.oldTime()
+        );
+
+        fluxFieldType phiCorr
+        (
+            phi.oldTime() - (mesh().Sf() & fvc::interpolate(rhoU0))
+        );
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr)
+               *rDeltaT*phiCorr
+            )
+        );
+    }
+    else if
+    (
+        U.dimensions() == rho.dimensions()*dimVelocity
+     && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        fluxFieldType phiCorr
+        (
+            phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+        );
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+               *rDeltaT*phiCorr
             )
         );
     }
@@ -651,7 +714,6 @@ CoEulerDdtScheme<Type>::fvcDdtPhiCorr
 
         return fluxFieldType::null();
     }
-    }
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.H

@@ -158,16 +158,27 @@ public:
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
 
-        tmp<fluxFieldType> fvcDdtPhiCorr
+        tmp<fluxFieldType> fvcDdtUfCorr
         (
-            const volScalarField& rA,
             const GeometricField<Type, fvPatchField, volMesh>& U,
-            const fluxFieldType& phi
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+        );
+
+        tmp<fluxFieldType> fvcDdtPhiCorr
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const fluxFieldType& phi
+        );
+
+        tmp<fluxFieldType> fvcDdtUfCorr
+        (
+            const volScalarField& rho,
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
         );
 
         tmp<fluxFieldType> fvcDdtPhiCorr
         (
-            const volScalarField& rA,
             const volScalarField& rho,
             const GeometricField<Type, fvPatchField, volMesh>& U,
             const fluxFieldType& phi
@@ -180,10 +191,17 @@ public:
 };
 
 
+template<>
+tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtUfCorr
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& U,
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
+);
+
+
 template<>
 tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& U,
     const surfaceScalarField& phi
 );
@@ -192,7 +210,6 @@ tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtPhiCorr
 template<>
 tmp<surfaceScalarField> CoEulerDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& rho,
     const volScalarField& U,
     const surfaceScalarField& phi

src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C

@@ -646,7 +646,6 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
     fvMatrix<Type>& fvm = tfvm();
 
     scalar rDtCoef = rDtCoef_(ddt0).value();
-
     fvm.diag() = rDtCoef*mesh().V();
 
     vf.oldTime().oldTime();
@@ -875,215 +874,34 @@ CrankNicolsonDdtScheme<Type>::fvmDdt
 }
 
 
-
 template<class Type>
 tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
-CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
+CrankNicolsonDdtScheme<Type>::fvcDdtUfCorr
 (
-    const volScalarField& rA,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
 )
 {
-    DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& dUdt0 =
+    DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
         ddt0_<GeometricField<Type, fvPatchField, volMesh> >
         (
             "ddt0(" + U.name() + ')',
             U.dimensions()
         );
 
-    DDt0Field<fluxFieldType>& dphidt0 =
-        ddt0_<fluxFieldType>
-        (
-            "ddt0(" + phi.name() + ')',
-            phi.dimensions()
-        );
+    dimensionedScalar rDeltaT = rDtCoef_(ddt0);
 
     IOobject ddtIOobject
     (
-        "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
+        "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
         mesh().time().timeName(),
         mesh()
     );
 
-    dimensionedScalar rDtCoef = rDtCoef_(dUdt0);
-
-    if (mesh().moving())
-    {
-        return tmp<fluxFieldType>
+    fluxFieldType phiCorr
     (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rDtCoef.dimensions()*rA.dimensions()*phi.dimensions(),
-                    pTraits<typename flux<Type>::type>::zero
-                )
-            )
+        mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
     );
-    }
-    else
-    {
-        if (evaluate(dUdt0))
-        {
-            dUdt0 =
-                rDtCoef0_(dUdt0)*(U.oldTime() - U.oldTime().oldTime())
-              - offCentre_(dUdt0());
-        }
-
-        if (evaluate(dphidt0))
-        {
-            dphidt0 =
-                rDtCoef0_(dphidt0)*(phi.oldTime() - phi.oldTime().oldTime())
-              - offCentre_(dphidt0());
-        }
-
-        return tmp<fluxFieldType>
-        (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
-               *fvc::interpolate(rA)
-               *(
-                    (rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
-                  - (
-                        fvc::interpolate
-                        (
-                            (rDtCoef*U.oldTime() + offCentre_(dUdt0()))
-                        ) & mesh().Sf()
-                    )
-                )
-            )
-        );
-    }
-}
-
-
-template<class Type>
-tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
-CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
-(
-    const volScalarField& rA,
-    const volScalarField& rho,
-    const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
-)
-{
-    DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& dUdt0 =
-        ddt0_<GeometricField<Type, fvPatchField, volMesh> >
-        (
-            "ddt0(" + U.name() + ')',
-            U.dimensions()
-        );
-
-    DDt0Field<fluxFieldType>& dphidt0 =
-        ddt0_<fluxFieldType>
-        (
-            "ddt0(" + phi.name() + ')',
-            U.dimensions()*dimArea
-        );
-
-    IOobject ddtIOobject
-    (
-        "ddtPhiCorr("
-      + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
-        mesh().time().timeName(),
-        mesh()
-    );
-
-    dimensionedScalar rDtCoef = rDtCoef_(dUdt0);
-
-    if (mesh().moving())
-    {
-        return tmp<fluxFieldType>
-        (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
-                    pTraits<typename flux<Type>::type>::zero
-                )
-            )
-        );
-    }
-    else
-    {
-        if
-        (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == dimVelocity*dimArea
-        )
-        {
-            if (evaluate(dUdt0))
-            {
-                dUdt0 = rDtCoef0_(dUdt0)*
-                (
-                    U.oldTime() - U.oldTime().oldTime()
-                ) - offCentre_(dUdt0());
-            }
-
-            if (evaluate(dphidt0))
-            {
-                dphidt0 = rDtCoef0_(dphidt0)*
-                (
-                    phi.oldTime()
-                  - fvc::interpolate(rho.oldTime().oldTime()/rho.oldTime())
-                   *phi.oldTime().oldTime()
-                ) - offCentre_(dphidt0());
-            }
-
-            return tmp<fluxFieldType>
-            (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
-                    (
-                        fvc::interpolate(rA*rho.oldTime())
-                       *(rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
-                      - (
-                            fvc::interpolate
-                            (
-                                rA*rho.oldTime()
-                               *(rDtCoef*U.oldTime() + offCentre_(dUdt0()))
-                            ) & mesh().Sf()
-                        )
-                    )
-                )
-            );
-        }
-        else if
-        (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
-        )
-        {
-            if (evaluate(dUdt0))
-            {
-                dUdt0 = rDtCoef0_(dUdt0)*
-                (
-                    U.oldTime() - U.oldTime().oldTime()
-                ) - offCentre_(dUdt0());
-            }
-
-            if (evaluate(dphidt0))
-            {
-                dphidt0 = rDtCoef0_(dphidt0)*
-                (
-                    phi.oldTime()
-                   /fvc::interpolate(rho.oldTime())
-                  - phi.oldTime().oldTime()
-                   /fvc::interpolate(rho.oldTime().oldTime())
-                ) - offCentre_(dphidt0());
-            }
 
     return tmp<fluxFieldType>
     (
@@ -1093,24 +911,216 @@ CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
             this->fvcDdtPhiCoeff
             (
                 U.oldTime(),
-                        phi.oldTime()/fvc::interpolate(rho.oldTime())
-                    )
-                   *(
-                        fvc::interpolate(rA*rho.oldTime())
-                       *(
-                           rDtCoef*phi.oldTime()/fvc::interpolate(rho.oldTime())
-                         + offCentre_(dphidt0())
-                        )
-                      - (
-                            fvc::interpolate
-                            (
-                                rA*rho.oldTime()
-                               *(rDtCoef*U.oldTime() + offCentre_(dUdt0()))
-                            ) & mesh().Sf()
-                        )
+                mesh().Sf() & Uf.oldTime(),
+                phiCorr
             )
+           *rDeltaT*phiCorr
         )
     );
+}
+
+
+template<class Type>
+tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
+CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
+        ddt0_<GeometricField<Type, fvPatchField, volMesh> >
+        (
+            "ddt0(" + U.name() + ')',
+            U.dimensions()
+        );
+
+    dimensionedScalar rDeltaT = rDtCoef_(ddt0);
+
+    IOobject ddtIOobject
+    (
+        "ddtCorr(" + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    fluxFieldType phiCorr
+    (
+        phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+    );
+
+    return tmp<fluxFieldType>
+    (
+        new fluxFieldType
+        (
+            ddtIOobject,
+            this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+           *rDeltaT*phiCorr
+        )
+    );
+}
+
+
+template<class Type>
+tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
+CrankNicolsonDdtScheme<Type>::fvcDdtUfCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+)
+{
+    IOobject ddtIOobject
+    (
+        "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    if
+    (
+        U.dimensions() == dimVelocity
+     && Uf.dimensions() == rho.dimensions()*dimVelocity
+    )
+    {
+        DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
+            ddt0_<GeometricField<Type, fvPatchField, volMesh> >
+            (
+                "ddt0(" + rho.name() + ',' + U.name() + ')',
+                U.dimensions()
+            );
+
+        dimensionedScalar rDeltaT = rDtCoef_(ddt0);
+
+        GeometricField<Type, fvPatchField, volMesh> rhoU0
+        (
+            rho.oldTime()*U.oldTime()
+        );
+
+        fluxFieldType phiCorr
+        (
+            mesh().Sf() & (Uf.oldTime() - fvc::interpolate(rhoU0))
+        );
+
+        tmp<fluxFieldType> ddtCorr
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff
+                (
+                    rhoU0,
+                    mesh().Sf() & Uf.oldTime(),
+                    phiCorr
+                )
+               *rDeltaT*phiCorr
+            )
+        );
+
+        return ddtCorr;
+    }
+    else if
+    (
+        U.dimensions() == rho.dimensions()*dimVelocity
+     && Uf.dimensions() == rho.dimensions()*dimVelocity
+    )
+    {
+        DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
+            ddt0_<GeometricField<Type, fvPatchField, volMesh> >
+            (
+                "ddt0(" + U.name() + ')',
+                U.dimensions()
+            );
+
+        dimensionedScalar rDeltaT = rDtCoef_(ddt0);
+
+        fluxFieldType phiCorr
+        (
+             mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+        );
+
+        tmp<fluxFieldType> ddtCorr
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff
+                (
+                    U.oldTime(),
+                    mesh().Sf() & Uf.oldTime(),
+                    phiCorr
+                )
+               *rDeltaT*phiCorr
+            )
+        );
+
+        return ddtCorr;
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr"
+        )   << "dimensions of Uf are not correct"
+            << abort(FatalError);
+
+        return fluxFieldType::null();
+    }
+}
+
+
+template<class Type>
+tmp<typename CrankNicolsonDdtScheme<Type>::fluxFieldType>
+CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    IOobject ddtIOobject
+    (
+        "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    if
+    (
+        U.dimensions() == dimVelocity
+     && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
+            ddt0_<GeometricField<Type, fvPatchField, volMesh> >
+            (
+                "ddt0(" + rho.name() + ',' + U.name() + ')',
+                U.dimensions()
+            );
+
+        dimensionedScalar rDeltaT = rDtCoef_(ddt0);
+
+        GeometricField<Type, fvPatchField, volMesh> rhoU0
+        (
+            rho.oldTime()*U.oldTime()
+        );
+
+        fluxFieldType phiCorr
+        (
+            phi.oldTime() - (mesh().Sf() & fvc::interpolate(rhoU0))
+        );
+
+        tmp<fluxFieldType> ddtCorr
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr)
+               *rDeltaT*phiCorr
+            )
+        );
+
+        return ddtCorr;
     }
     else if
     (
@@ -1118,41 +1128,31 @@ CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
      && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
     )
     {
-            if (evaluate(dUdt0))
-            {
-                dUdt0 = rDtCoef0_(dUdt0)*
+        DDt0Field<GeometricField<Type, fvPatchField, volMesh> >& ddt0 =
+            ddt0_<GeometricField<Type, fvPatchField, volMesh> >
             (
-                    U.oldTime() - U.oldTime().oldTime()
-                ) - offCentre_(dUdt0());
-            }
+                "ddt0(" + U.name() + ')',
+                U.dimensions()
+            );
 
-            if (evaluate(dphidt0))
-            {
-                dphidt0 = rDtCoef0_(dphidt0)*
+        dimensionedScalar rDeltaT = rDtCoef_(ddt0);
+
+        fluxFieldType phiCorr
         (
-                    phi.oldTime() - phi.oldTime().oldTime()
-                ) - offCentre_(dphidt0());
-            }
+            phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+        );
 
-            return tmp<fluxFieldType>
+        tmp<fluxFieldType> ddtCorr
         (
             new fluxFieldType
             (
                 ddtIOobject,
-                    this->fvcDdtPhiCoeff
-                    (rho.oldTime(), U.oldTime(), phi.oldTime())
-                  * (
-                        fvc::interpolate(rA)
-                       *(rDtCoef*phi.oldTime() + offCentre_(dphidt0()))
-                      - (
-                            fvc::interpolate
-                            (
-                                rA*(rDtCoef*U.oldTime() + offCentre_(dUdt0()))
-                            ) & mesh().Sf()
-                        )
-                    )
+                this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+               *rDeltaT*phiCorr
             )
         );
+
+        return ddtCorr;
     }
     else
     {
@@ -1164,7 +1164,6 @@ CrankNicolsonDdtScheme<Type>::fvcDdtPhiCorr
 
         return fluxFieldType::null();
     }
-    }
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.H

@@ -229,16 +229,27 @@ public:
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
 
-        tmp<fluxFieldType> fvcDdtPhiCorr
+        tmp<fluxFieldType> fvcDdtUfCorr
         (
-            const volScalarField& rA,
             const GeometricField<Type, fvPatchField, volMesh>& U,
-            const fluxFieldType& phi
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+        );
+
+        tmp<fluxFieldType> fvcDdtPhiCorr
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const fluxFieldType& phi
+        );
+
+        tmp<fluxFieldType> fvcDdtUfCorr
+        (
+            const volScalarField& rho,
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
         );
 
         tmp<fluxFieldType> fvcDdtPhiCorr
         (
-            const volScalarField& rA,
             const volScalarField& rho,
             const GeometricField<Type, fvPatchField, volMesh>& U,
             const fluxFieldType& phi
@@ -252,10 +263,17 @@ public:
 };
 
 
+template<>
+tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtUfCorr
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& U,
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
+);
+
+
 template<>
 tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& U,
     const surfaceScalarField& phi
 );
@@ -264,7 +282,6 @@ tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtPhiCorr
 template<>
 tmp<surfaceScalarField> CrankNicolsonDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& rho,
     const volScalarField& U,
     const surfaceScalarField& phi

src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C

@@ -370,34 +370,38 @@ EulerDdtScheme<Type>::fvmDdt
 
 template<class Type>
 tmp<typename EulerDdtScheme<Type>::fluxFieldType>
-EulerDdtScheme<Type>::fvcDdtPhiCorr
+EulerDdtScheme<Type>::fvcDdtUfCorr
 (
-    const volScalarField& rA,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phiAbs
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
 )
 {
     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
 
     IOobject ddtIOobject
     (
-        "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phiAbs.name() + ')',
+        "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
         mesh().time().timeName(),
         mesh()
     );
 
-    tmp<fluxFieldType> phiCorr =
-        phiAbs.oldTime() - (fvc::interpolate(U.oldTime()) & mesh().Sf());
-
-    phiCorr().boundaryField() = pTraits<typename flux<Type>::type>::zero;
+    fluxFieldType phiCorr
+    (
+        mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+    );
 
     return tmp<fluxFieldType>
     (
         new fluxFieldType
         (
             ddtIOobject,
-            this->fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime(), phiCorr())
-          * fvc::interpolate(rDeltaT*rA)*phiCorr
+            this->fvcDdtPhiCoeff
+            (
+                U.oldTime(),
+                mesh().Sf() & Uf.oldTime(),
+                phiCorr
+            )
+           *rDeltaT*phiCorr
         )
     );
 }
@@ -407,21 +411,50 @@ template<class Type>
 tmp<typename EulerDdtScheme<Type>::fluxFieldType>
 EulerDdtScheme<Type>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
-    const volScalarField& rho,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phiAbs
+    const fluxFieldType& phi
 )
 {
     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
 
     IOobject ddtIOobject
     (
-        "ddtPhiCorr("
-      + rA.name() + ','
-      + rho.name() + ','
-      + U.name() + ','
-      + phiAbs.name() + ')',
+        "ddtCorr(" + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    fluxFieldType phiCorr
+    (
+        phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+    );
+
+    return tmp<fluxFieldType>
+    (
+        new fluxFieldType
+        (
+            ddtIOobject,
+            this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+           *rDeltaT*phiCorr
+        )
+    );
+}
+
+
+template<class Type>
+tmp<typename EulerDdtScheme<Type>::fluxFieldType>
+EulerDdtScheme<Type>::fvcDdtUfCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+)
+{
+    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
         mesh().time().timeName(),
         mesh()
     );
@@ -429,95 +462,144 @@ EulerDdtScheme<Type>::fvcDdtPhiCorr
     if
     (
         U.dimensions() == dimVelocity
-     && phiAbs.dimensions() == dimVelocity*dimArea
+     && Uf.dimensions() == rho.dimensions()*dimVelocity
     )
     {
-        tmp<fluxFieldType> ddtPhiCorr
+        GeometricField<Type, fvPatchField, volMesh> rhoU0
+        (
+            rho.oldTime()*U.oldTime()
+        );
+
+        fluxFieldType phiCorr
+        (
+            mesh().Sf() & (Uf.oldTime() - fvc::interpolate(rhoU0))
+        );
+
+        return tmp<fluxFieldType>
         (
             new fluxFieldType
             (
                 ddtIOobject,
-                rDeltaT
-              * this->fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime())
-              * (
-                    fvc::interpolate(rA*rho.oldTime())*phiAbs.oldTime()
-                  - (
-                        fvc::interpolate(rA*rho.oldTime()*U.oldTime())
-                      & mesh().Sf()
-                    )
+                this->fvcDdtPhiCoeff
+                (
+                    rhoU0,
+                    mesh().Sf() & Uf.oldTime(),
+                    phiCorr
                 )
+               *rDeltaT*phiCorr
             )
         );
-
-        ddtPhiCorr().boundaryField() = pTraits<typename flux<Type>::type>::zero;
-
-        return ddtPhiCorr;
-    }
-    else if
-    (
-        U.dimensions() == dimVelocity
-     && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
-    )
-    {
-        tmp<fluxFieldType> ddtPhiCorr
-        (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                rDeltaT
-              * this->fvcDdtPhiCoeff
-                (
-                    U.oldTime(),
-                    phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
-                )
-              * (
-                    fvc::interpolate(rA*rho.oldTime())
-                  * phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
-                 - (
-                       fvc::interpolate
-                       (
-                           rA*rho.oldTime()*U.oldTime()
-                       ) & mesh().Sf()
-                   )
-                )
-            )
-        );
-
-        ddtPhiCorr().boundaryField() = pTraits<typename flux<Type>::type>::zero;
-
-        return ddtPhiCorr;
     }
     else if
     (
         U.dimensions() == rho.dimensions()*dimVelocity
-     && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+     && Uf.dimensions() == rho.dimensions()*dimVelocity
     )
     {
-        tmp<fluxFieldType> ddtPhiCorr
+        fluxFieldType phiCorr
+        (
+             mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+        );
+
+        return tmp<fluxFieldType>
         (
             new fluxFieldType
             (
                 ddtIOobject,
-                rDeltaT
-              * this->fvcDdtPhiCoeff
-                (rho.oldTime(), U.oldTime(), phiAbs.oldTime())
-              * (
-                    fvc::interpolate(rA)*phiAbs.oldTime()
-                  - (fvc::interpolate(rA*U.oldTime()) & mesh().Sf())
+                this->fvcDdtPhiCoeff
+                (
+                    U.oldTime(),
+                    mesh().Sf() & Uf.oldTime(),
+                    phiCorr
                 )
+               *rDeltaT*phiCorr
             )
         );
-
-        ddtPhiCorr().boundaryField() = pTraits<typename flux<Type>::type>::zero;
-
-        return ddtPhiCorr;
     }
     else
     {
         FatalErrorIn
         (
             "EulerDdtScheme<Type>::fvcDdtPhiCorr"
-        )   << "dimensions of phiAbs are not correct"
+        )   << "dimensions of Uf are not correct"
+            << abort(FatalError);
+
+        return fluxFieldType::null();
+    }
+}
+
+
+template<class Type>
+tmp<typename EulerDdtScheme<Type>::fluxFieldType>
+EulerDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    if
+    (
+        U.dimensions() == dimVelocity
+     && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        GeometricField<Type, fvPatchField, volMesh> rhoU0
+        (
+            rho.oldTime()*U.oldTime()
+        );
+
+        fluxFieldType phiCorr
+        (
+            phi.oldTime() - (mesh().Sf() & fvc::interpolate(rhoU0))
+        );
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr)
+               *rDeltaT*phiCorr
+            )
+        );
+    }
+    else if
+    (
+        U.dimensions() == rho.dimensions()*dimVelocity
+     && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        fluxFieldType phiCorr
+        (
+            phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+        );
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+               *rDeltaT*phiCorr
+            )
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "EulerDdtScheme<Type>::fvcDdtPhiCorr"
+        )   << "dimensions of phi are not correct"
             << abort(FatalError);
 
         return fluxFieldType::null();

src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.H

@@ -136,16 +136,27 @@ public:
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
 
-        tmp<fluxFieldType> fvcDdtPhiCorr
+        tmp<fluxFieldType> fvcDdtUfCorr
         (
-            const volScalarField& rA,
             const GeometricField<Type, fvPatchField, volMesh>& U,
-            const fluxFieldType& phi
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+        );
+
+        tmp<fluxFieldType> fvcDdtPhiCorr
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const fluxFieldType& phi
+        );
+
+        tmp<fluxFieldType> fvcDdtUfCorr
+        (
+            const volScalarField& rho,
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
         );
 
         tmp<fluxFieldType> fvcDdtPhiCorr
         (
-            const volScalarField& rA,
             const volScalarField& rho,
             const GeometricField<Type, fvPatchField, volMesh>& U,
             const fluxFieldType& phi
@@ -158,10 +169,17 @@ public:
 };
 
 
+template<>
+tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtUfCorr
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& U,
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
+);
+
+
 template<>
 tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& U,
     const surfaceScalarField& phi
 );
@@ -170,7 +188,6 @@ tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtPhiCorr
 template<>
 tmp<surfaceScalarField> EulerDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& rho,
     const volScalarField& U,
     const surfaceScalarField& phi

src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C

@@ -483,122 +483,26 @@ SLTSDdtScheme<Type>::fvmDdt
 
 template<class Type>
 tmp<typename SLTSDdtScheme<Type>::fluxFieldType>
-SLTSDdtScheme<Type>::fvcDdtPhiCorr
+SLTSDdtScheme<Type>::fvcDdtUfCorr
 (
-    const volScalarField& rA,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
 )
 {
     IOobject ddtIOobject
     (
-        "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
+        "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
         mesh().time().timeName(),
         mesh()
     );
 
-    if (mesh().moving())
-    {
-        return tmp<fluxFieldType>
-        (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rA.dimensions()*phi.dimensions()/dimTime,
-                    pTraits<typename flux<Type>::type>::zero
-                )
-            )
-        );
-    }
-    else
-    {
-        const volScalarField rDeltaT(SLrDeltaT());
+    const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
 
-        return tmp<fluxFieldType>
+    fluxFieldType phiCorr
     (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
-                (
-                    fvc::interpolate(rDeltaT*rA)*phi.oldTime()
-                  - (fvc::interpolate(rDeltaT*rA*U.oldTime()) & mesh().Sf())
-                )
-            )
-        );
-    }
-}
-
-
-template<class Type>
-tmp<typename SLTSDdtScheme<Type>::fluxFieldType>
-SLTSDdtScheme<Type>::fvcDdtPhiCorr
-(
-    const volScalarField& rA,
-    const volScalarField& rho,
-    const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
-)
-{
-    IOobject ddtIOobject
-    (
-        "ddtPhiCorr("
-      + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
-        mesh().time().timeName(),
-        mesh()
+        mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
     );
 
-    if (mesh().moving())
-    {
-        return tmp<fluxFieldType>
-        (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
-                    pTraits<typename flux<Type>::type>::zero
-                )
-            )
-        );
-    }
-    else
-    {
-        const volScalarField rDeltaT(SLrDeltaT());
-
-        if
-        (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == dimVelocity*dimArea
-        )
-        {
-            return tmp<fluxFieldType>
-            (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
-                   *(
-                        fvc::interpolate(rDeltaT*rA*rho.oldTime())*phi.oldTime()
-                      - (fvc::interpolate(rDeltaT*rA*rho.oldTime()*U.oldTime())
-                      & mesh().Sf())
-                    )
-                )
-            );
-        }
-        else if
-        (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == dimDensity*dimVelocity*dimArea
-        )
-        {
     return tmp<fluxFieldType>
     (
         new fluxFieldType
@@ -607,40 +511,199 @@ SLTSDdtScheme<Type>::fvcDdtPhiCorr
             this->fvcDdtPhiCoeff
             (
                 U.oldTime(),
-                        phi.oldTime()/fvc::interpolate(rho.oldTime())
-                    )
-                   *(
-                        fvc::interpolate(rDeltaT*rA*rho.oldTime())
-                       *phi.oldTime()/fvc::interpolate(rho.oldTime())
-                      - (
-                            fvc::interpolate
-                            (
-                                rDeltaT*rA*rho.oldTime()*U.oldTime()
-                            ) & mesh().Sf()
-                        )
+                (mesh().Sf() & Uf.oldTime()),
+                phiCorr
             )
+           *rDeltaT*phiCorr
         )
     );
-        }
-        else if
+}
+
+
+template<class Type>
+tmp<typename SLTSDdtScheme<Type>::fluxFieldType>
+SLTSDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    IOobject ddtIOobject
     (
-            U.dimensions() == dimDensity*dimVelocity
-         && phi.dimensions() == dimDensity*dimVelocity*dimArea
+        "ddtCorr(" + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
+
+    fluxFieldType phiCorr
+    (
+        phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+    );
+
+    return tmp<fluxFieldType>
+    (
+        new fluxFieldType
+        (
+            ddtIOobject,
+            this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+           *rDeltaT*phiCorr
+        )
+    );
+}
+
+
+template<class Type>
+tmp<typename SLTSDdtScheme<Type>::fluxFieldType>
+SLTSDdtScheme<Type>::fvcDdtUfCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+)
+{
+    IOobject ddtIOobject
+    (
+        "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    const surfaceScalarField rDeltaT(fvc::interpolate(SLrDeltaT()));
+
+    if
+    (
+        U.dimensions() == dimVelocity
+     && Uf.dimensions() == dimDensity*dimVelocity
     )
     {
+        GeometricField<Type, fvPatchField, volMesh> rhoU0
+        (
+            rho.oldTime()*U.oldTime()
+        );
+
+        fluxFieldType phiCorr
+        (
+            mesh().Sf() & (Uf.oldTime() - fvc::interpolate(rhoU0))
+        );
+
         return tmp<fluxFieldType>
         (
             new fluxFieldType
             (
                 ddtIOobject,
                 this->fvcDdtPhiCoeff
-                    (rho.oldTime(), U.oldTime(), phi.oldTime())
-                  * (
-                        fvc::interpolate(rDeltaT*rA)*phi.oldTime()
-                      - (
-                            fvc::interpolate(rDeltaT*rA*U.oldTime())&mesh().Sf()
+                (
+                    rhoU0,
+                    mesh().Sf() & Uf.oldTime(),
+                    phiCorr
                 )
+               *rDeltaT*phiCorr
             )
+        );
+    }
+    else if
+    (
+        U.dimensions() == dimDensity*dimVelocity
+     && Uf.dimensions() == dimDensity*dimVelocity
+    )
+    {
+        fluxFieldType phiCorr
+        (
+             mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+        );
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff
+                (
+                    U.oldTime(),
+                    mesh().Sf() & Uf.oldTime(),
+                    phiCorr
+                )
+               *rDeltaT*phiCorr
+            )
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "SLTSDdtScheme<Type>::fvcDdtPhiCorr"
+        )   << "dimensions of Uf are not correct"
+            << abort(FatalError);
+
+        return fluxFieldType::null();
+    }
+}
+
+
+template<class Type>
+tmp<typename SLTSDdtScheme<Type>::fluxFieldType>
+SLTSDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    if
+    (
+        U.dimensions() == dimVelocity
+     && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        GeometricField<Type, fvPatchField, volMesh> rhoU0
+        (
+            rho.oldTime()*U.oldTime()
+        );
+
+        fluxFieldType phiCorr
+        (
+            phi.oldTime() - (mesh().Sf() & fvc::interpolate(rhoU0))
+        );
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr)
+               *rDeltaT*phiCorr
+            )
+        );
+    }
+    else if
+    (
+        U.dimensions() == rho.dimensions()*dimVelocity
+     && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        fluxFieldType phiCorr
+        (
+            phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+        );
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+               *rDeltaT*phiCorr
             )
         );
     }
@@ -654,7 +717,6 @@ SLTSDdtScheme<Type>::fvcDdtPhiCorr
 
         return fluxFieldType::null();
     }
-    }
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.H

@@ -159,16 +159,27 @@ public:
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
 
-        tmp<fluxFieldType> fvcDdtPhiCorr
+        tmp<fluxFieldType> fvcDdtUfCorr
         (
-            const volScalarField& rA,
             const GeometricField<Type, fvPatchField, volMesh>& U,
-            const fluxFieldType& phi
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+        );
+
+        tmp<fluxFieldType> fvcDdtPhiCorr
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const fluxFieldType& phi
+        );
+
+        tmp<fluxFieldType> fvcDdtUfCorr
+        (
+            const volScalarField& rho,
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
         );
 
         tmp<fluxFieldType> fvcDdtPhiCorr
         (
-            const volScalarField& rA,
             const volScalarField& rho,
             const GeometricField<Type, fvPatchField, volMesh>& U,
             const fluxFieldType& phi
@@ -181,10 +192,17 @@ public:
 };
 
 
+template<>
+tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtUfCorr
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& U,
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
+);
+
+
 template<>
 tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& U,
     const surfaceScalarField& phi
 );
@@ -193,7 +211,6 @@ tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtPhiCorr
 template<>
 tmp<surfaceScalarField> SLTSDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& rho,
     const volScalarField& U,
     const surfaceScalarField& phi

src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C

@@ -514,18 +514,17 @@ backwardDdtScheme<Type>::fvmDdt
 
 template<class Type>
 tmp<typename backwardDdtScheme<Type>::fluxFieldType>
-backwardDdtScheme<Type>::fvcDdtPhiCorr
+backwardDdtScheme<Type>::fvcDdtUfCorr
 (
-    const volScalarField& rA,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
 )
 {
     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
 
     IOobject ddtIOobject
     (
-        "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
+        "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
         mesh().time().timeName(),
         mesh()
     );
@@ -542,22 +541,16 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
         new fluxFieldType
         (
             ddtIOobject,
-            rDeltaT*this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
+            this->fvcDdtPhiCoeff(U.oldTime(), (mesh().Sf() & Uf.oldTime()))
+           *rDeltaT
            *(
-                fvc::interpolate(rA)
-               *(
-                   coefft0*phi.oldTime()
-                 - coefft00*phi.oldTime().oldTime()
-                )
-              - (
-                    fvc::interpolate
+                mesh().Sf()
+              & (
+                    (coefft0*Uf.oldTime() - coefft00*Uf.oldTime().oldTime())
+                  - fvc::interpolate
                     (
-                        rA*
-                        (
-                            coefft0*U.oldTime()
-                          - coefft00*U.oldTime().oldTime()
+                        coefft0*U.oldTime() - coefft00*U.oldTime().oldTime()
                     )
-                    ) & mesh().Sf()
                 )
             )
         )
@@ -569,21 +562,62 @@ template<class Type>
 tmp<typename backwardDdtScheme<Type>::fluxFieldType>
 backwardDdtScheme<Type>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
-    const volScalarField& rho,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phiAbs
+    const fluxFieldType& phi
 )
 {
     dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
 
     IOobject ddtIOobject
     (
-        "ddtPhiCorr("
-      + rA.name() + ','
-      + rho.name() + ','
-      + U.name() + ','
-      + phiAbs.name() + ')',
+        "ddtCorr(" + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    scalar deltaT = deltaT_();
+    scalar deltaT0 = deltaT0_(U);
+
+    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    scalar coefft0  = coefft + coefft00;
+
+    return tmp<fluxFieldType>
+    (
+        new fluxFieldType
+        (
+            ddtIOobject,
+            this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
+           *rDeltaT
+           *(
+                (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime())
+              - (
+                    mesh().Sf()
+                  & fvc::interpolate
+                    (
+                        coefft0*U.oldTime() - coefft00*U.oldTime().oldTime()
+                    )
+                )
+            )
+        )
+    );
+}
+
+
+template<class Type>
+tmp<typename backwardDdtScheme<Type>::fluxFieldType>
+backwardDdtScheme<Type>::fvcDdtUfCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+)
+{
+    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
         mesh().time().timeName(),
         mesh()
     );
@@ -598,68 +632,31 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
     if
     (
         U.dimensions() == dimVelocity
-     && phiAbs.dimensions() == dimVelocity*dimArea
+     && Uf.dimensions() == rho.dimensions()*dimVelocity
     )
     {
-        return tmp<fluxFieldType>
+        GeometricField<Type, fvPatchField, volMesh> rhoU0
         (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                rDeltaT*this->fvcDdtPhiCoeff(U.oldTime(), phiAbs.oldTime())
-               *(
-                    coefft0*fvc::interpolate(rA*rho.oldTime())
-                   *phiAbs.oldTime()
-                  - coefft00*fvc::interpolate(rA*rho.oldTime().oldTime())
-                   *phiAbs.oldTime().oldTime()
-                  - (
-                        fvc::interpolate
-                        (
-                            rA*
-                            (
-                                coefft0*rho.oldTime()*U.oldTime()
-                              - coefft00*rho.oldTime().oldTime()
-                               *U.oldTime().oldTime()
-                            )
-                        ) & mesh().Sf()
-                    )
-                )
-            )
+            rho.oldTime()*U.oldTime()
         );
-    }
-    else if
+
+        GeometricField<Type, fvPatchField, volMesh> rhoU00
         (
-        U.dimensions() == dimVelocity
-     && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
-    )
-    {
+            rho.oldTime().oldTime()*U.oldTime().oldTime()
+        );
+
         return tmp<fluxFieldType>
         (
             new fluxFieldType
             (
                 ddtIOobject,
-                rDeltaT
-               *this->fvcDdtPhiCoeff
-                (
-                    U.oldTime(),
-                    phiAbs.oldTime()/fvc::interpolate(rho.oldTime())
-                )
+                this->fvcDdtPhiCoeff(rhoU0, mesh().Sf() & Uf.oldTime())
+               *rDeltaT
                *(
-                    fvc::interpolate(rA)
-                   *(
-                       coefft0*phiAbs.oldTime()
-                     - coefft00*phiAbs.oldTime().oldTime()
-                    )
-                  - (
-                        fvc::interpolate
-                        (
-                            rA*
-                            (
-                                coefft0*rho.oldTime()*U.oldTime()
-                              - coefft00*rho.oldTime().oldTime()
-                               *U.oldTime().oldTime()
-                            )
-                        ) & mesh().Sf()
+                    mesh().Sf()
+                  & (
+                        (coefft0*Uf.oldTime() - coefft00*Uf.oldTime().oldTime())
+                      - fvc::interpolate(coefft0*rhoU0 - coefft00*rhoU00)
                     )
                 )
             )
@@ -668,32 +665,25 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
     else if
     (
         U.dimensions() == rho.dimensions()*dimVelocity
-     && phiAbs.dimensions() == rho.dimensions()*dimVelocity*dimArea
+     && Uf.dimensions() == rho.dimensions()*dimVelocity
     )
     {
+
         return tmp<fluxFieldType>
         (
             new fluxFieldType
             (
                 ddtIOobject,
-                rDeltaT
-              * this->fvcDdtPhiCoeff
-                (rho.oldTime(), U.oldTime(), phiAbs.oldTime())
-              * (
-                    fvc::interpolate(rA)
-                  * (
-                       coefft0*phiAbs.oldTime()
-                     - coefft00*phiAbs.oldTime().oldTime()
-                    )
-                  - (
-                        fvc::interpolate
+                this->fvcDdtPhiCoeff(U.oldTime(), mesh().Sf() & Uf.oldTime())
+               *rDeltaT
+               *(
+                    mesh().Sf()
+                  & (
+                        (coefft0*Uf.oldTime() - coefft00*Uf.oldTime().oldTime())
+                      - fvc::interpolate
                         (
-                            rA*
-                            (
-                                coefft0*U.oldTime()
-                              - coefft00*U.oldTime().oldTime()
+                            coefft0*U.oldTime() - coefft00*U.oldTime().oldTime()
                         )
-                        ) & mesh().Sf()
                     )
                 )
             )
@@ -704,7 +694,105 @@ backwardDdtScheme<Type>::fvcDdtPhiCorr
         FatalErrorIn
         (
             "backwardDdtScheme<Type>::fvcDdtPhiCorr"
-        )   << "dimensions of phiAbs are not correct"
+        )   << "dimensions of phi are not correct"
+            << abort(FatalError);
+
+        return fluxFieldType::null();
+    }
+}
+
+
+template<class Type>
+tmp<typename backwardDdtScheme<Type>::fluxFieldType>
+backwardDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    scalar deltaT = deltaT_();
+    scalar deltaT0 = deltaT0_(U);
+
+    scalar coefft   = 1 + deltaT/(deltaT + deltaT0);
+    scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
+    scalar coefft0  = coefft + coefft00;
+
+    if
+    (
+        U.dimensions() == dimVelocity
+     && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        GeometricField<Type, fvPatchField, volMesh> rhoU0
+        (
+            rho.oldTime()*U.oldTime()
+        );
+
+        GeometricField<Type, fvPatchField, volMesh> rhoU00
+        (
+            rho.oldTime().oldTime()*U.oldTime().oldTime()
+        );
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff(rhoU0, phi.oldTime())
+               *rDeltaT
+               *(
+                    (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime())
+                  - (
+                        mesh().Sf()
+                      & fvc::interpolate(coefft0*rhoU0 - coefft00*rhoU00)
+                    )
+                )
+            )
+        );
+    }
+    else if
+    (
+        U.dimensions() == rho.dimensions()*dimVelocity
+     && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
+               *rDeltaT
+               *(
+                    (coefft0*phi.oldTime() - coefft00*phi.oldTime().oldTime())
+                  - (
+                        mesh().Sf()
+                      & fvc::interpolate
+                        (
+                            coefft0*U.oldTime() - coefft00*U.oldTime().oldTime()
+                        )
+                    )
+                )
+            )
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "backwardDdtScheme<Type>::fvcDdtPhiCorr"
+        )   << "dimensions of phi are not correct"
             << abort(FatalError);
 
         return fluxFieldType::null();

src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.H

@@ -147,22 +147,32 @@ public:
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
 
-        tmp<fluxFieldType> fvcDdtPhiCorr
+        tmp<fluxFieldType> fvcDdtUfCorr
         (
-            const volScalarField& rA,
             const GeometricField<Type, fvPatchField, volMesh>& U,
-            const fluxFieldType& phi
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+        );
+
+        tmp<fluxFieldType> fvcDdtPhiCorr
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const fluxFieldType& phi
+        );
+
+        tmp<fluxFieldType> fvcDdtUfCorr
+        (
+            const volScalarField& rho,
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
         );
 
         tmp<fluxFieldType> fvcDdtPhiCorr
         (
-            const volScalarField& rA,
             const volScalarField& rho,
             const GeometricField<Type, fvPatchField, volMesh>& U,
             const fluxFieldType& phi
         );
 
-
         tmp<surfaceScalarField> meshPhi
         (
             const GeometricField<Type, fvPatchField, volMesh>&
@@ -170,10 +180,17 @@ public:
 };
 
 
+template<>
+tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtUfCorr
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& U,
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
+);
+
+
 template<>
 tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& U,
     const surfaceScalarField& phi
 );
@@ -182,7 +199,6 @@ tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtPhiCorr
 template<>
 tmp<surfaceScalarField> backwardDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& rho,
     const volScalarField& U,
     const surfaceScalarField& phi

src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.C

@@ -124,14 +124,38 @@ boundedDdtScheme<Type>::fvmDdt
 
 template<class Type>
 tmp<typename boundedDdtScheme<Type>::fluxFieldType>
-boundedDdtScheme<Type>::fvcDdtPhiCorr
+boundedDdtScheme<Type>::fvcDdtUfCorr
 (
-    const volScalarField& rA,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
 )
 {
-    return scheme_().fvcDdtPhiCorr(rA, U, phi);
+    return scheme_().fvcDdtUfCorr(U, Uf);
+}
+
+
+template<class Type>
+tmp<typename boundedDdtScheme<Type>::fluxFieldType>
+boundedDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    return scheme_().fvcDdtPhiCorr(U, phi);
+}
+
+
+template<class Type>
+tmp<typename boundedDdtScheme<Type>::fluxFieldType>
+boundedDdtScheme<Type>::fvcDdtUfCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+)
+{
+    return scheme_().fvcDdtUfCorr(rho, U, Uf);
 }
 
 
@@ -139,13 +163,12 @@ template<class Type>
 tmp<typename boundedDdtScheme<Type>::fluxFieldType>
 boundedDdtScheme<Type>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& rho,
     const GeometricField<Type, fvPatchField, volMesh>& U,
     const fluxFieldType& phi
 )
 {
-    return scheme_().fvcDdtPhiCorr(rA, rho, U, phi);
+    return scheme_().fvcDdtPhiCorr(rho, U, phi);
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/boundedDdtScheme/boundedDdtScheme.H

@@ -145,16 +145,27 @@ public:
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
 
-        tmp<fluxFieldType> fvcDdtPhiCorr
+        tmp<fluxFieldType> fvcDdtUfCorr
         (
-            const volScalarField& rA,
             const GeometricField<Type, fvPatchField, volMesh>& U,
-            const fluxFieldType& phi
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+        );
+
+        tmp<fluxFieldType> fvcDdtPhiCorr
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const fluxFieldType& phi
+        );
+
+        tmp<fluxFieldType> fvcDdtUfCorr
+        (
+            const volScalarField& rho,
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
         );
 
         tmp<fluxFieldType> fvcDdtPhiCorr
         (
-            const volScalarField& rA,
             const volScalarField& rho,
             const GeometricField<Type, fvPatchField, volMesh>& U,
             const fluxFieldType& phi
@@ -167,10 +178,17 @@ public:
 };
 
 
+template<>
+tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtUfCorr
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& U,
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
+);
+
+
 template<>
 tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& U,
     const surfaceScalarField& phi
 );
@@ -179,7 +197,6 @@ tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtPhiCorr
 template<>
 tmp<surfaceScalarField> boundedDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& rho,
     const volScalarField& U,
     const surfaceScalarField& phi

src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C

@@ -141,82 +141,7 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
     const fluxFieldType& phi
 )
 {
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
-
-    tmp<surfaceScalarField> tddtCouplingCoeff = scalar(1)
-      - min
-        (
-            mag(phi - (mesh().Sf() & fvc::interpolate(U)))
-           /(mag(phi) + dimensionedScalar("small", phi.dimensions(), VSMALL)),
-           //(rDeltaT*mesh().magSf()/mesh().deltaCoeffs()),
-            scalar(1)
-        );
-
-    surfaceScalarField& ddtCouplingCoeff = tddtCouplingCoeff();
-
-    forAll(U.boundaryField(), patchi)
-    {
-        if (U.boundaryField()[patchi].fixesValue())
-        {
-            ddtCouplingCoeff.boundaryField()[patchi] = 0.0;
-        }
-    }
-
-    if (debug > 1)
-    {
-        Info<< "ddtCouplingCoeff mean max min = "
-            << gAverage(ddtCouplingCoeff.internalField())
-            << " " << gMax(ddtCouplingCoeff.internalField())
-            << " " << gMin(ddtCouplingCoeff.internalField())
-            << endl;
-    }
-
-    return tddtCouplingCoeff;
-}
-
-
-template<class Type>
-tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
-(
-    const volScalarField& rho,
-    const GeometricField<Type, fvPatchField, volMesh>& rhoU,
-    const fluxFieldType& phi
-)
-{
-    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
-
-    tmp<surfaceScalarField> tddtCouplingCoeff = scalar(1)
-      - min
-        (
-            mag(phi - (mesh().Sf() & fvc::interpolate(rhoU)))
-           /(
-                mag(phi) + dimensionedScalar("small", phi.dimensions(), VSMALL)
-                //fvc::interpolate(rho)*rDeltaT
-                //*mesh().magSf()/mesh().deltaCoeffs()
-            ),
-            scalar(1)
-        );
-
-    surfaceScalarField& ddtCouplingCoeff = tddtCouplingCoeff();
-
-    forAll(rhoU.boundaryField(), patchi)
-    {
-        if (rhoU.boundaryField()[patchi].fixesValue())
-        {
-            ddtCouplingCoeff.boundaryField()[patchi] = 0.0;
-        }
-    }
-
-    if (debug > 1)
-    {
-        Info<< "ddtCouplingCoeff mean max min = "
-            << gAverage(ddtCouplingCoeff.internalField())
-            << " " << gMax(ddtCouplingCoeff.internalField())
-            << " " << gMin(ddtCouplingCoeff.internalField())
-            << endl;
-    }
-
-    return tddtCouplingCoeff;
+    return fvcDdtPhiCoeff(U, phi, phi - (mesh().Sf() & fvc::interpolate(U)));
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H

@@ -198,29 +198,32 @@ public:
             const fluxFieldType& phi
         );
 
+        virtual tmp<fluxFieldType> fvcDdtUfCorr
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+        ) = 0;
+
         virtual tmp<fluxFieldType> fvcDdtPhiCorr
         (
-            const volScalarField& rA,
             const GeometricField<Type, fvPatchField, volMesh>& U,
             const fluxFieldType& phi
         ) = 0;
 
-        tmp<surfaceScalarField> fvcDdtPhiCoeff
+        virtual tmp<fluxFieldType> fvcDdtUfCorr
         (
             const volScalarField& rho,
-            const GeometricField<Type, fvPatchField, volMesh>& rhoU,
-            const fluxFieldType& phi
-        );
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+        ) = 0;
 
         virtual tmp<fluxFieldType> fvcDdtPhiCorr
         (
-            const volScalarField& rA,
             const volScalarField& rho,
             const GeometricField<Type, fvPatchField, volMesh>& U,
             const fluxFieldType& phi
         ) = 0;
 
-
         virtual tmp<surfaceScalarField> meshPhi
         (
             const GeometricField<Type, fvPatchField, volMesh>&
@@ -257,9 +260,19 @@ makeFvDdtTypeScheme(SS, symmTensor)                                            \
 makeFvDdtTypeScheme(SS, tensor)                                                \
                                                                                \
 template<>                                                                     \
+tmp<surfaceScalarField> SS<scalar>::fvcDdtUfCorr                               \
+(                                                                              \
+    const volScalarField& U,                                                   \
+    const surfaceScalarField& Uf                                               \
+)                                                                              \
+{                                                                              \
+    notImplemented(#SS"<scalar>::fvcDdtUfCorr");                               \
+    return surfaceScalarField::null();                                         \
+}                                                                              \
+                                                                               \
+template<>                                                                     \
 tmp<surfaceScalarField> SS<scalar>::fvcDdtPhiCorr                              \
 (                                                                              \
-    const volScalarField& rA,                                                  \
     const volScalarField& U,                                                   \
     const surfaceScalarField& phi                                              \
 )                                                                              \
@@ -269,9 +282,20 @@ tmp<surfaceScalarField> SS<scalar>::fvcDdtPhiCorr                              \
 }                                                                              \
                                                                                \
 template<>                                                                     \
+tmp<surfaceScalarField> SS<scalar>::fvcDdtUfCorr                               \
+(                                                                              \
+    const volScalarField& rho,                                                 \
+    const volScalarField& U,                                                   \
+    const surfaceScalarField& Uf                                               \
+)                                                                              \
+{                                                                              \
+    notImplemented(#SS"<scalar>::fvcDdtUfCorr");                               \
+    return surfaceScalarField::null();                                         \
+}                                                                              \
+                                                                               \
+template<>                                                                     \
 tmp<surfaceScalarField> SS<scalar>::fvcDdtPhiCorr                              \
 (                                                                              \
-    const volScalarField& rA,                                                  \
     const volScalarField& rho,                                                 \
     const volScalarField& U,                                                   \
     const surfaceScalarField& phi                                              \

src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C

@@ -377,122 +377,26 @@ localEulerDdtScheme<Type>::fvmDdt
 
 template<class Type>
 tmp<typename localEulerDdtScheme<Type>::fluxFieldType>
-localEulerDdtScheme<Type>::fvcDdtPhiCorr
+localEulerDdtScheme<Type>::fvcDdtUfCorr
 (
-    const volScalarField& rA,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
 )
 {
     IOobject ddtIOobject
     (
-        "ddtPhiCorr(" + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
+        "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
         mesh().time().timeName(),
         mesh()
     );
 
-    if (mesh().moving())
-    {
-        return tmp<fluxFieldType>
-        (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rA.dimensions()*phi.dimensions()/dimTime,
-                    pTraits<typename flux<Type>::type>::zero
-                )
-            )
-        );
-    }
-    else
-    {
-        const volScalarField& rDeltaT = localRDeltaT();
+    const surfaceScalarField rDeltaT(fvc::interpolate(localRDeltaT()));
 
-        return tmp<fluxFieldType>
+    fluxFieldType phiCorr
     (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())*
-                (
-                    fvc::interpolate(rDeltaT*rA)*phi.oldTime()
-                  - (fvc::interpolate(rDeltaT*rA*U.oldTime()) & mesh().Sf())
-                )
-            )
-        );
-    }
-}
-
-
-template<class Type>
-tmp<typename localEulerDdtScheme<Type>::fluxFieldType>
-localEulerDdtScheme<Type>::fvcDdtPhiCorr
-(
-    const volScalarField& rA,
-    const volScalarField& rho,
-    const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
-)
-{
-    IOobject ddtIOobject
-    (
-        "ddtPhiCorr("
-      + rA.name() + ',' + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
-        mesh().time().timeName(),
-        mesh()
+        mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
     );
 
-    if (mesh().moving())
-    {
-        return tmp<fluxFieldType>
-        (
-            new fluxFieldType
-            (
-                ddtIOobject,
-                mesh(),
-                dimensioned<typename flux<Type>::type>
-                (
-                    "0",
-                    rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
-                    pTraits<typename flux<Type>::type>::zero
-                )
-            )
-        );
-    }
-    else
-    {
-        const volScalarField& rDeltaT = localRDeltaT();
-
-        if
-        (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == dimVelocity*dimArea
-        )
-        {
-            return tmp<fluxFieldType>
-            (
-                new fluxFieldType
-                (
-                    ddtIOobject,
-                    this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime())
-                   *(
-                        fvc::interpolate(rDeltaT*rA*rho.oldTime())*phi.oldTime()
-                      - (fvc::interpolate(rDeltaT*rA*rho.oldTime()*U.oldTime())
-                      & mesh().Sf())
-                    )
-                )
-            );
-        }
-        else if
-        (
-            U.dimensions() == dimVelocity
-         && phi.dimensions() == dimDensity*dimVelocity*dimArea
-        )
-        {
     return tmp<fluxFieldType>
     (
         new fluxFieldType
@@ -501,40 +405,199 @@ localEulerDdtScheme<Type>::fvcDdtPhiCorr
             this->fvcDdtPhiCoeff
             (
                 U.oldTime(),
-                        phi.oldTime()/fvc::interpolate(rho.oldTime())
-                    )
-                   *(
-                        fvc::interpolate(rDeltaT*rA*rho.oldTime())
-                       *phi.oldTime()/fvc::interpolate(rho.oldTime())
-                      - (
-                            fvc::interpolate
-                            (
-                                rDeltaT*rA*rho.oldTime()*U.oldTime()
-                            ) & mesh().Sf()
-                        )
+                (mesh().Sf() & Uf.oldTime()),
+                phiCorr
             )
+           *rDeltaT*phiCorr
         )
     );
-        }
-        else if
+}
+
+
+template<class Type>
+tmp<typename localEulerDdtScheme<Type>::fluxFieldType>
+localEulerDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    IOobject ddtIOobject
     (
-            U.dimensions() == dimDensity*dimVelocity
-         && phi.dimensions() == dimDensity*dimVelocity*dimArea
+        "ddtCorr(" + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    const surfaceScalarField rDeltaT(fvc::interpolate(localRDeltaT()));
+
+    fluxFieldType phiCorr
+    (
+        phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+    );
+
+    return tmp<fluxFieldType>
+    (
+        new fluxFieldType
+        (
+            ddtIOobject,
+            this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+           *rDeltaT*phiCorr
+        )
+    );
+}
+
+
+template<class Type>
+tmp<typename localEulerDdtScheme<Type>::fluxFieldType>
+localEulerDdtScheme<Type>::fvcDdtUfCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+)
+{
+    IOobject ddtIOobject
+    (
+        "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    const surfaceScalarField rDeltaT(fvc::interpolate(localRDeltaT()));
+
+    if
+    (
+        U.dimensions() == dimVelocity
+     && Uf.dimensions() == dimDensity*dimVelocity
     )
     {
+        GeometricField<Type, fvPatchField, volMesh> rhoU0
+        (
+            rho.oldTime()*U.oldTime()
+        );
+
+        fluxFieldType phiCorr
+        (
+            mesh().Sf() & (Uf.oldTime() - fvc::interpolate(rhoU0))
+        );
+
         return tmp<fluxFieldType>
         (
             new fluxFieldType
             (
                 ddtIOobject,
                 this->fvcDdtPhiCoeff
-                    (rho.oldTime(), U.oldTime(), phi.oldTime())
-                  * (
-                        fvc::interpolate(rDeltaT*rA)*phi.oldTime()
-                      - (
-                            fvc::interpolate(rDeltaT*rA*U.oldTime())&mesh().Sf()
+                (
+                    rhoU0,
+                    mesh().Sf() & Uf.oldTime(),
+                    phiCorr
                 )
+               *rDeltaT*phiCorr
             )
+        );
+    }
+    else if
+    (
+        U.dimensions() == dimDensity*dimVelocity
+     && Uf.dimensions() == dimDensity*dimVelocity
+    )
+    {
+        fluxFieldType phiCorr
+        (
+             mesh().Sf() & (Uf.oldTime() - fvc::interpolate(U.oldTime()))
+        );
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff
+                (
+                    U.oldTime(),
+                    mesh().Sf() & Uf.oldTime(),
+                    phiCorr
+                )
+               *rDeltaT*phiCorr
+            )
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "localEulerDdtScheme<Type>::fvcDdtPhiCorr"
+        )   << "dimensions of Uf are not correct"
+            << abort(FatalError);
+
+        return fluxFieldType::null();
+    }
+}
+
+
+template<class Type>
+tmp<typename localEulerDdtScheme<Type>::fluxFieldType>
+localEulerDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
+
+    IOobject ddtIOobject
+    (
+        "ddtCorr(" + rho.name() + ',' + U.name() + ',' + phi.name() + ')',
+        mesh().time().timeName(),
+        mesh()
+    );
+
+    if
+    (
+        U.dimensions() == dimVelocity
+     && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        GeometricField<Type, fvPatchField, volMesh> rhoU0
+        (
+            rho.oldTime()*U.oldTime()
+        );
+
+        fluxFieldType phiCorr
+        (
+            phi.oldTime() - (mesh().Sf() & fvc::interpolate(rhoU0))
+        );
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff(rhoU0, phi.oldTime(), phiCorr)
+               *rDeltaT*phiCorr
+            )
+        );
+    }
+    else if
+    (
+        U.dimensions() == rho.dimensions()*dimVelocity
+     && phi.dimensions() == rho.dimensions()*dimVelocity*dimArea
+    )
+    {
+        fluxFieldType phiCorr
+        (
+            phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
+        );
+
+        return tmp<fluxFieldType>
+        (
+            new fluxFieldType
+            (
+                ddtIOobject,
+                this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
+               *rDeltaT*phiCorr
             )
         );
     }
@@ -548,7 +611,6 @@ localEulerDdtScheme<Type>::fvcDdtPhiCorr
 
         return fluxFieldType::null();
     }
-    }
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.H

@@ -148,16 +148,27 @@ public:
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
 
-        tmp<fluxFieldType> fvcDdtPhiCorr
+        tmp<fluxFieldType> fvcDdtUfCorr
         (
-            const volScalarField& rA,
             const GeometricField<Type, fvPatchField, volMesh>& U,
-            const fluxFieldType& phi
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+        );
+
+        tmp<fluxFieldType> fvcDdtPhiCorr
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const fluxFieldType& phi
+        );
+
+        tmp<fluxFieldType> fvcDdtUfCorr
+        (
+            const volScalarField& rho,
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
         );
 
         tmp<fluxFieldType> fvcDdtPhiCorr
         (
-            const volScalarField& rA,
             const volScalarField& rho,
             const GeometricField<Type, fvPatchField, volMesh>& U,
             const fluxFieldType& phi
@@ -170,10 +181,17 @@ public:
 };
 
 
+template<>
+tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtUfCorr
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& U,
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
+);
+
+
 template<>
 tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& U,
     const surfaceScalarField& phi
 );
@@ -182,7 +200,6 @@ tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtPhiCorr
 template<>
 tmp<surfaceScalarField> localEulerDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& rho,
     const volScalarField& U,
     const surfaceScalarField& phi

src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.C

@@ -221,11 +221,10 @@ steadyStateDdtScheme<Type>::fvmDdt
 
 template<class Type>
 tmp<typename steadyStateDdtScheme<Type>::fluxFieldType>
-steadyStateDdtScheme<Type>::fvcDdtPhiCorr
+steadyStateDdtScheme<Type>::fvcDdtUfCorr
 (
-    const volScalarField& rA,
     const GeometricField<Type, fvPatchField, volMesh>& U,
-    const fluxFieldType& phi
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
 )
 {
     return tmp<fluxFieldType>
@@ -234,8 +233,7 @@ steadyStateDdtScheme<Type>::fvcDdtPhiCorr
         (
             IOobject
             (
-                "ddtPhiCorr("
-              + rA.name() + ',' + U.name() + ',' + phi.name() + ')',
+                "ddtCorr(" + U.name() + ',' + Uf.name() + ')',
                 mesh().time().timeName(),
                 mesh()
             ),
@@ -243,7 +241,70 @@ steadyStateDdtScheme<Type>::fvcDdtPhiCorr
             dimensioned<typename flux<Type>::type>
             (
                 "0",
-                rA.dimensions()*phi.dimensions()/dimTime,
+                mesh().Sf().dimensions()*Uf.dimensions()*dimArea/dimTime,
+                pTraits<typename flux<Type>::type>::zero
+            )
+        )
+    );
+}
+
+
+template<class Type>
+tmp<typename steadyStateDdtScheme<Type>::fluxFieldType>
+steadyStateDdtScheme<Type>::fvcDdtPhiCorr
+(
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const fluxFieldType& phi
+)
+{
+    return tmp<fluxFieldType>
+    (
+        new fluxFieldType
+        (
+            IOobject
+            (
+                "ddtCorr(" + U.name() + ',' + phi.name() + ')',
+                mesh().time().timeName(),
+                mesh()
+            ),
+            mesh(),
+            dimensioned<typename flux<Type>::type>
+            (
+                "0",
+                phi.dimensions()/dimTime,
+                pTraits<typename flux<Type>::type>::zero
+            )
+        )
+    );
+}
+
+
+template<class Type>
+tmp<typename steadyStateDdtScheme<Type>::fluxFieldType>
+steadyStateDdtScheme<Type>::fvcDdtUfCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+)
+{
+    return tmp<fluxFieldType>
+    (
+        new fluxFieldType
+        (
+            IOobject
+            (
+                "ddtCorr("
+              + rho.name()
+              + ',' + U.name() + ',' + Uf.name() + ')',
+                mesh().time().timeName(),
+                mesh()
+            ),
+            mesh(),
+            dimensioned<typename flux<Type>::type>
+            (
+                "0",
+                rho.dimensions()*Uf.dimensions()*dimArea/dimTime,
                 pTraits<typename flux<Type>::type>::zero
             )
         )
@@ -255,7 +316,6 @@ template<class Type>
 tmp<typename steadyStateDdtScheme<Type>::fluxFieldType>
 steadyStateDdtScheme<Type>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& rho,
     const GeometricField<Type, fvPatchField, volMesh>& U,
     const fluxFieldType& phi
@@ -267,8 +327,8 @@ steadyStateDdtScheme<Type>::fvcDdtPhiCorr
         (
             IOobject
             (
-                "ddtPhiCorr("
-              + rA.name() + ',' + rho.name()
+                "ddtCorr("
+              + rho.name()
               + ',' + U.name() + ',' + phi.name() + ')',
                 mesh().time().timeName(),
                 mesh()
@@ -277,7 +337,7 @@ steadyStateDdtScheme<Type>::fvcDdtPhiCorr
             dimensioned<typename flux<Type>::type>
             (
                 "0",
-                rA.dimensions()*rho.dimensions()*phi.dimensions()/dimTime,
+                rho.dimensions()*phi.dimensions()/dimTime,
                 pTraits<typename flux<Type>::type>::zero
             )
         )

src/finiteVolume/finiteVolume/ddtSchemes/steadyStateDdtScheme/steadyStateDdtScheme.H

@@ -135,16 +135,27 @@ public:
 
         typedef typename ddtScheme<Type>::fluxFieldType fluxFieldType;
 
-        tmp<fluxFieldType> fvcDdtPhiCorr
+        tmp<fluxFieldType> fvcDdtUfCorr
         (
-            const volScalarField& rA,
             const GeometricField<Type, fvPatchField, volMesh>& U,
-            const fluxFieldType& phi
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+        );
+
+        tmp<fluxFieldType> fvcDdtPhiCorr
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const fluxFieldType& phi
+        );
+
+        tmp<fluxFieldType> fvcDdtUfCorr
+        (
+            const volScalarField& rho,
+            const GeometricField<Type, fvPatchField, volMesh>& U,
+            const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
         );
 
         tmp<fluxFieldType> fvcDdtPhiCorr
         (
-            const volScalarField& rA,
             const volScalarField& rho,
             const GeometricField<Type, fvPatchField, volMesh>& U,
             const fluxFieldType& phi
@@ -157,10 +168,17 @@ public:
 };
 
 
+template<>
+tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtUfCorr
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& U,
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& Uf
+);
+
+
 template<>
 tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& U,
     const surfaceScalarField& phi
 );
@@ -169,7 +187,6 @@ tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtPhiCorr
 template<>
 tmp<surfaceScalarField> steadyStateDdtScheme<scalar>::fvcDdtPhiCorr
 (
-    const volScalarField& rA,
     const volScalarField& rho,
     const volScalarField& U,
     const surfaceScalarField& phi

src/finiteVolume/finiteVolume/fvc/fvcDdt.C

@@ -128,9 +128,24 @@ ddt
 
 template<class Type>
 tmp<GeometricField<typename flux<Type>::type, fvsPatchField, surfaceMesh> >
-ddtPhiCorr
+ddtCorr
+(
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+)
+{
+    return fv::ddtScheme<Type>::New
+    (
+        U.mesh(),
+        U.mesh().ddtScheme("ddt(" + U.name() + ')')
+    )().fvcDdtUfCorr(U, Uf);
+}
+
+
+template<class Type>
+tmp<GeometricField<typename flux<Type>::type, fvsPatchField, surfaceMesh> >
+ddtCorr
 (
-    const volScalarField& rA,
     const GeometricField<Type, fvPatchField, volMesh>& U,
     const GeometricField
     <
@@ -144,15 +159,32 @@ ddtPhiCorr
     (
         U.mesh(),
         U.mesh().ddtScheme("ddt(" + U.name() + ')')
-    )().fvcDdtPhiCorr(rA, U, phi);
+    )().fvcDdtPhiCorr(U, phi);
 }
 
 
 template<class Type>
 tmp<GeometricField<typename flux<Type>::type, fvsPatchField, surfaceMesh> >
-ddtPhiCorr
+ddtCorr
+(
+    const volScalarField& rho,
+    const GeometricField<Type, fvPatchField, volMesh>& U,
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+)
+{
+    return fv::ddtScheme<Type>::New
+    (
+        U.mesh(),
+        U.mesh().ddtScheme("ddt(" + U.name() + ')')
+    )().fvcDdtPhiCorr(rho, U, U.mesh().Sf() & Uf);
+    //***HGW fvcDdtUfCorr(rho, U, Uf);
+}
+
+
+template<class Type>
+tmp<GeometricField<typename flux<Type>::type, fvsPatchField, surfaceMesh> >
+ddtCorr
 (
-    const volScalarField& rA,
     const volScalarField& rho,
     const GeometricField<Type, fvPatchField, volMesh>& U,
     const GeometricField
@@ -167,7 +199,7 @@ ddtPhiCorr
     (
         U.mesh(),
         U.mesh().ddtScheme("ddt(" + rho.name() + ',' + U.name() + ')')
-    )().fvcDdtPhiCorr(rA, rho, U, phi);
+    )().fvcDdtPhiCorr(rho, U, phi);
 }
 
 

src/finiteVolume/finiteVolume/fvc/fvcDdt.H

@@ -113,9 +113,24 @@ namespace fvc
             surfaceMesh
         >
     >
-    ddtPhiCorr
+    ddtCorr
+    (
+        const GeometricField<Type, fvPatchField, volMesh>& U,
+        const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+    );
+
+    template<class Type>
+    tmp
+    <
+        GeometricField
+        <
+            typename Foam::flux<Type>::type,
+            fvsPatchField,
+            surfaceMesh
+        >
+    >
+    ddtCorr
     (
-        const volScalarField& rA,
         const GeometricField<Type, fvPatchField, volMesh>& U,
         const GeometricField
         <
@@ -135,9 +150,25 @@ namespace fvc
             surfaceMesh
         >
     >
-    ddtPhiCorr
+    ddtCorr
+    (
+        const volScalarField& rho,
+        const GeometricField<Type, fvPatchField, volMesh>& U,
+        const GeometricField<Type, fvsPatchField, surfaceMesh>& Uf
+    );
+
+    template<class Type>
+    tmp
+    <
+        GeometricField
+        <
+            typename Foam::flux<Type>::type,
+            fvsPatchField,
+            surfaceMesh
+        >
+    >
+    ddtCorr
     (
-        const volScalarField& rA,
         const volScalarField& rho,
         const GeometricField<Type, fvPatchField, volMesh>& U,
         const GeometricField

tutorials/DNS/dnsFoam/boxTurb16/system/fvSchemes

@@ -23,8 +23,6 @@ ddtSchemes
 gradSchemes
 {
     default         Gauss linear;
-    grad(p)         Gauss linear;
-    grad(U)         Gauss linear;
 }
 
 divSchemes
@@ -35,15 +33,12 @@ divSchemes
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(nu,U) Gauss linear corrected;
-    laplacian(Dp,p) Gauss linear corrected;
+    default         Gauss linear corrected;
 }
 
 interpolationSchemes
 {
     default         linear;
-    interpolate(U)  linear;
 }
 
 snGradSchemes

tutorials/basic/potentialFoam/cylinder/system/fvSchemes

@@ -32,8 +32,7 @@ divSchemes
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(1,p)  Gauss linear corrected;
+    default         Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/basic/potentialFoam/pitzDaily/system/fvSchemes

@@ -32,8 +32,7 @@ divSchemes
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(1,p)  Gauss linear corrected;
+    default         Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/fvSchemes

@@ -56,15 +56,7 @@ divSchemes
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(muEff,U) Gauss linear limited corrected 0.333;
-    laplacian(DkEff,k) Gauss linear limited corrected 0.333;
-    laplacian(DepsilonEff,epsilon) Gauss linear limited corrected 0.333;
-    laplacian((rho*inv((((1)*A(U))+tDragDcu))),p) Gauss linear limited corrected 0.333;
-    laplacian(Db,b) Gauss linear limited corrected 0.333;
-    laplacian(Db,ft) Gauss linear limited corrected 0.333;
-    laplacian(Db,ha) Gauss linear limited corrected 0.333;
-    laplacian(Db,hau) Gauss linear limited corrected 0.333;
+    default         Gauss linear limited corrected 0.333;
 }
 
 interpolationSchemes

tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/system/fvSchemes

@@ -54,16 +54,7 @@ divSchemes
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(muEff,U) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
-    laplacian(DREff,R) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
-    laplacian(alphaEff,b) Gauss linear corrected;
-    laplacian(alphaEff,ft) Gauss linear corrected;
-    laplacian(alphaEff,ha) Gauss linear corrected;
-    laplacian(alphaEff,hau) Gauss linear corrected;
+    default         Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/combustion/fireFoam/les/oppositeBurningPanels/system/fvSchemes

@@ -46,14 +46,7 @@ divSchemes
 
 laplacianSchemes
 {
-    //default         none;
     default         Gauss linear corrected;
-    laplacian(muEff,U) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(thermo:alpha,h) Gauss linear corrected;
-    laplacian((((rho*(1|A(U)))*rho)*gh)) Gauss linear corrected;
-    laplacian(interpolate((rho*(1|A(U)))),p) Gauss linear corrected;
-    laplacian(gammaRad,G) Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/system/fvSchemes

@@ -23,7 +23,6 @@ ddtSchemes
 gradSchemes
 {
     default         Gauss linear;
-    grad(p)         Gauss linear;
 }
 
 divSchemes
@@ -45,15 +44,7 @@ divSchemes
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(muEff,U) Gauss linear corrected;
-    laplacian(mut,U) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
-    laplacian(DREff,R) Gauss linear corrected;
-    laplacian(DomegaEff,omega) Gauss linear corrected;
-    laplacian(Dp,p) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
+    default         Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/compressible/rhoPimpleDyMFoam/annularThermalMixer/constant/polyMesh/boundary

@@ -22,77 +22,77 @@ FoamFile
         type            patch;
         inGroups        1(inlet);
         nFaces          544;
-        startFace       244948;
+        startFace       245836;
     }
     outerInlet
     {
         type            patch;
         inGroups        1(inlet);
-        nFaces          1396;
-        startFace       245492;
+        nFaces          1404;
+        startFace       246380;
     }
     innerOutlet
     {
         type            patch;
         inGroups        1(outlet);
         nFaces          544;
-        startFace       246888;
+        startFace       247784;
     }
     outerOutlet
     {
         type            patch;
         inGroups        1(outlet);
-        nFaces          1396;
-        startFace       247432;
+        nFaces          1404;
+        startFace       248328;
     }
     rotorBlades
     {
         type            wall;
         inGroups        1(movingWalls);
         nFaces          540;
-        startFace       248828;
+        startFace       249732;
     }
     rotorBlades_slave
     {
         type            wall;
         inGroups        1(movingWalls);
         nFaces          540;
-        startFace       249368;
+        startFace       250272;
     }
     shaft
     {
         type            wall;
         inGroups        1(movingWalls);
-        nFaces          1052;
-        startFace       249908;
+        nFaces          1044;
+        startFace       250812;
     }
     statorBlades
     {
         type            wall;
         inGroups        1(staticWalls);
         nFaces          2128;
-        startFace       250960;
+        startFace       251856;
     }
     statorBlades_slave
     {
         type            wall;
         inGroups        1(staticWalls);
         nFaces          2128;
-        startFace       253088;
+        startFace       253984;
     }
     wall
     {
         type            wall;
         inGroups        1(staticWalls);
-        nFaces          5808;
-        startFace       255216;
+        nFaces          6165;
+        startFace       256112;
     }
     AMI1
     {
         type            cyclicAMI;
         inGroups        1(cyclicAMI);
         nFaces          10944;
-        startFace       261024;
+        startFace       262277;
         matchTolerance  0.0001;
         transform       noOrdering;
         neighbourPatch  AMI2;
@@ -102,7 +102,7 @@ FoamFile
         type            cyclicAMI;
         inGroups        1(cyclicAMI);
         nFaces          10944;
-        startFace       271968;
+        startFace       273221;
         matchTolerance  0.0001;
         transform       noOrdering;
         neighbourPatch  AMI1;

tutorials/compressible/rhoPimpleFoam/les/pitzDaily/constant/polyMesh/boundary

@@ -44,6 +44,7 @@ FoamFile
     frontAndBack
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          24450;
         startFace       24730;
     }

tutorials/compressible/rhoPimpleFoam/les/pitzDaily/system/fvSchemes

@@ -40,13 +40,7 @@ divSchemes
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(muEff,U) Gauss linear corrected;
-    laplacian(Dp,p) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(DBEff,B) Gauss linear corrected;
-    laplacian(DmuTildaEff,muTilda) Gauss linear corrected;
+    default         Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/compressible/rhoPimpleFoam/ras/angledDuct/system/fvSchemes

@@ -23,7 +23,6 @@ ddtSchemes
 gradSchemes
 {
     default         Gauss linear;
-    grad(p)         Gauss linear;
 }
 
 divSchemes
@@ -45,15 +44,7 @@ divSchemes
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(muEff,U) Gauss linear corrected;
-    laplacian(mut,U) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
-    laplacian(DREff,R) Gauss linear corrected;
-    laplacian(DomegaEff,omega) Gauss linear corrected;
-    laplacian(Dp,p) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
+    default         Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/compressible/rhoPimpleFoam/ras/cavity/system/fvSchemes

@@ -23,7 +23,6 @@ ddtSchemes
 gradSchemes
 {
     default         Gauss linear;
-    grad(p)         Gauss linear;
 }
 
 divSchemes
@@ -45,15 +44,7 @@ divSchemes
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(muEff,U) Gauss linear orthogonal;
-    laplacian(mut,U) Gauss linear orthogonal;
-    laplacian(DkEff,k) Gauss linear orthogonal;
-    laplacian(DepsilonEff,epsilon) Gauss linear orthogonal;
-    laplacian(DREff,R) Gauss linear orthogonal;
-    laplacian(DomegaEff,omega) Gauss linear orthogonal;
-    laplacian(Dp,p) Gauss linear orthogonal;
-    laplacian(alphaEff,h) Gauss linear orthogonal;
+    default         Gauss linear orthogonal;
 }
 
 interpolationSchemes

tutorials/compressible/rhoPimpleFoam/ras/mixerVessel2D/system/fvSchemes

@@ -38,18 +38,12 @@ divSchemes
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(muEff,U) Gauss linear corrected;
-    laplacian(Dp,p) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+    default         Gauss linear corrected;
 }
 
 interpolationSchemes
 {
     default         linear;
-    interpolate(U)  linear;
 }
 
 snGradSchemes

tutorials/compressible/rhoPimplecFoam/angledDuct/system/fvSchemes

@@ -23,7 +23,6 @@ ddtSchemes
 gradSchemes
 {
     default         Gauss linear;
-    grad(p)         Gauss linear;
 }
 
 divSchemes
@@ -45,15 +44,7 @@ divSchemes
 
 laplacianSchemes
 {
-    default         none;
-    laplacian(muEff,U) Gauss linear corrected;
-    laplacian(mut,U) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
-    laplacian(DREff,R) Gauss linear corrected;
-    laplacian(DomegaEff,omega) Gauss linear corrected;
-    laplacian(Dp,p) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
+    default         Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSchemes

@@ -37,12 +37,7 @@ divSchemes
 
 laplacianSchemes
 {
-    laplacian(muEff,U) Gauss linear corrected;
-    laplacian(alphaEff,e) Gauss linear corrected;
-    laplacian((rho*rAU),p) Gauss linear corrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(1,p)  Gauss linear corrected;
+    default         Gauss linear corrected;
 }
 
 interpolationSchemes