Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/DNS/dnsFoam/dnsFoam.C

@@ -86,23 +86,28 @@ int main(int argc, char *argv[])
         for (int corr=1; corr<=1; corr++)
         {
             volScalarField rAU(1.0/UEqn.A());
+            volVectorField HbyA("HbyA", U);
+            HbyA = rAU*UEqn.H();
 
-            U = rAU*UEqn.H();
+            surfaceScalarField phiHbyA
-            phi = (fvc::interpolate(U) & mesh.Sf())
+            (
-                + fvc::ddtPhiCorr(rAU, U, phi);
+                "phiHbyA",
+                (fvc::interpolate(HbyA) & mesh.Sf())
+              + fvc::ddtPhiCorr(rAU, U, phi)
+            );
 
             fvScalarMatrix pEqn
             (
-                fvm::laplacian(rAU, p) == fvc::div(phi)
+                fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
             );
 
             pEqn.solve();
 
-            phi -= pEqn.flux();
+            phi = phiHbyA - pEqn.flux();
 
             #include "continuityErrs.H"
 
-            U -= rAU*fvc::grad(p);
+            U = HbyA - rAU*fvc::grad(p);
             U.correctBoundaryConditions();
         }
 

applications/solvers/combustion/XiFoam/pEqn.H

@@ -1,7 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -10,7 +11,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -34,19 +35,22 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
         fvc::interpolate(rho)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
         );
 
@@ -54,7 +58,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -62,7 +66,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 

applications/solvers/combustion/engineFoam/pEqn.H

@@ -1,7 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -9,7 +10,7 @@ if (pimple.transonic())
     (
         "phid",
         fvc::interpolate(psi)
-       *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
+       *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 
     while (pimple.correctNonOrthogonal())
@@ -31,15 +32,19 @@ if (pimple.transonic())
 }
 else
 {
-    phi = fvc::interpolate(rho)
+    surfaceScalarField phiHbyA
-         *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)
+       *((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
         );
 
@@ -47,7 +52,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -55,7 +60,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 

applications/solvers/combustion/fireFoam/pEqn.H

@@ -2,25 +2,29 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU));
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
-surfaceScalarField phiU
+surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+surfaceScalarField phiHbyA
 (
+    "phiHbyA",
     fvc::interpolate(rho)
    *(
-        (fvc::interpolate(U) & mesh.Sf())
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     )
+  - phig
 );
 
-phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
 while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix p_rghEqn
     (
         fvc::ddt(psi, rho)*gh
-      + fvc::div(phi)
+      + fvc::div(phiHbyA)
       + fvm::ddt(psi, p_rgh)
       - fvm::laplacian(rhorAUf, p_rgh)
      ==
@@ -32,7 +36,9 @@ while (pimple.correctNonOrthogonal())
 
     if (pimple.finalNonOrthogonalIter())
     {
-        phi += p_rghEqn.flux();
+        phi = phiHbyA + p_rghEqn.flux();
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() - phig)/rhorAUf);
+        U.correctBoundaryConditions();
     }
 }
 
@@ -41,8 +47,6 @@ p = p_rgh + rho*gh;
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
-U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 
 dpdt = fvc::ddt(p);

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -1,7 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -10,7 +11,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -34,19 +35,22 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
         fvc::interpolate(rho)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
         );
 
@@ -54,7 +58,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -62,7 +66,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 

applications/solvers/combustion/rhoReactingFoam/pEqn.H

@@ -6,27 +6,25 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
     if (pimple.transonic())
     {
-        surfaceScalarField phiv
+        surfaceScalarField phiHbyA
         (
-            (fvc::interpolate(U) & mesh.Sf())
+            "phiHbyA",
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-        phi = fvc::interpolate(rho)*phiv;
+        surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);
 
-        surfaceScalarField phid
+        phiHbyA *= fvc::interpolate(rho);
-        (
-            "phid",
-            fvc::interpolate(thermo.psi())*phiv
-        );
 
         fvScalarMatrix pDDtEqn
         (
-            fvc::ddt(rho) + fvc::div(phi)
+            fvc::ddt(rho) + fvc::div(phiHbyA)
           + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
         );
 
@@ -42,23 +40,26 @@
 
             if (pimple.finalNonOrthogonalIter())
             {
-                phi += pEqn.flux();
+                phi = phiHbyA + pEqn.flux();
             }
         }
     }
     else
     {
-        phi =
+        surfaceScalarField phiHbyA
+        (
+            "phiHbyA",
             fvc::interpolate(rho)
            *(
-                (fvc::interpolate(U) & mesh.Sf())
+                (fvc::interpolate(HbyA) & mesh.Sf())
               + fvc::ddtPhiCorr(rAU, rho, U, phi)
-            );
+            )
+        );
 
         fvScalarMatrix pDDtEqn
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
         );
 
         while (pimple.correctNonOrthogonal())
@@ -73,7 +74,7 @@
 
             if (pimple.finalNonOrthogonalIter())
             {
-                phi += pEqn.flux();
+                phi = phiHbyA + pEqn.flux();
             }
         }
     }
@@ -84,7 +85,7 @@
     #include "rhoEqn.H"
     #include "compressibleContinuityErrs.H"
 
-    U -= rAU*fvc::grad(p);
+    U = HbyA - rAU*fvc::grad(p);
     U.correctBoundaryConditions();
     K = 0.5*magSqr(U);
 

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -3,7 +3,8 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-U = rAU*UEqn().H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -17,7 +18,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -41,12 +42,15 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
+    surfaceScalarField phiHbyA
-        fvc::interpolate(rho)*
+    (
-        (
+        "phiHbyA",
-            (fvc::interpolate(U) & mesh.Sf())
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
+    );
 
     while (pimple.correctNonOrthogonal())
     {
@@ -54,7 +58,7 @@ else
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
         );
 
@@ -62,7 +66,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -81,7 +85,7 @@ rho.relax();
 Info<< "rho max/min : " << max(rho).value()
     << " " << min(rho).value() << endl;
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 

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

@@ -4,7 +4,8 @@ rho = min(rho, rhoMax);
 rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
-U = rAU*UEqn().H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -18,7 +19,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -43,13 +44,17 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
+    surfaceScalarField phiHbyA
-        fvc::interpolate(rho)*
+    (
-        (
+        "phiHbyA",
-            (fvc::interpolate(U) & mesh.Sf())
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
-    mrfZones.relativeFlux(fvc::interpolate(rho), phi);
+    );
+
+    mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -57,7 +62,7 @@ else
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
         );
 
@@ -65,7 +70,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -83,7 +88,7 @@ rho.relax();
 Info<< "rho max/min : " << max(rho).value()
     << " " << min(rho).value() << endl;
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -4,7 +4,9 @@ rho = min(rho, rhoMax);
 rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
-U = rAU*UEqn().H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn().H();
+
 UEqn.clear();
 
 bool closedVolume = false;
@@ -14,7 +16,7 @@ if (simple.transonic())
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
+        fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
     );
 
     while (simple.correctNonOrthogonal())
@@ -40,14 +42,19 @@ if (simple.transonic())
 }
 else
 {
-    phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
+    surfaceScalarField phiHbyA
-    closedVolume = adjustPhi(phi, U, p);
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+    );
+
+    closedVolume = adjustPhi(phiHbyA, U, p);
 
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(rho*rAU, p) == fvc::div(phi)
+            fvm::laplacian(rho*rAU, p) == fvc::div(phiHbyA)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -56,7 +63,7 @@ else
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi -= pEqn.flux();
+            phi = phiHbyA - pEqn.flux();
         }
     }
 }
@@ -67,7 +74,7 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
 // For closed-volume cases adjust the pressure and density levels

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

@@ -1,10 +1,12 @@
+volVectorField HbyA("HbyA", U);
+
 if (pressureImplicitPorosity)
 {
-    U = trTU()&UEqn().H();
+    HbyA = trTU() & UEqn().H();
 }
 else
 {
-    U = trAU()*UEqn().H();
+    HbyA = trAU()*UEqn().H();
 }
 
 UEqn.clear();
@@ -16,7 +18,7 @@ if (simple.transonic())
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
+        fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
     );
     mrfZones.relativeFlux(fvc::interpolate(psi), phid);
 
@@ -45,10 +47,15 @@ if (simple.transonic())
 }
 else
 {
-    phi = fvc::interpolate(rho*U) & mesh.Sf();
+    surfaceScalarField phiHbyA
-    mrfZones.relativeFlux(fvc::interpolate(rho), phi);
+    (
+        "phiHbyA",
+        fvc::interpolate(rho*HbyA) & mesh.Sf()
+    );
 
-    closedVolume = adjustPhi(phi, U, p);
+    mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
+    closedVolume = adjustPhi(phiHbyA, U, p);
 
     while (simple.correctNonOrthogonal())
     {
@@ -56,11 +63,11 @@ else
 
         if (pressureImplicitPorosity)
         {
-            tpEqn = (fvm::laplacian(rho*trTU(), p) == fvc::div(phi));
+            tpEqn = (fvm::laplacian(rho*trTU(), p) == fvc::div(phiHbyA));
         }
         else
         {
-            tpEqn = (fvm::laplacian(rho*trAU(), p) == fvc::div(phi));
+            tpEqn = (fvm::laplacian(rho*trAU(), p) == fvc::div(phiHbyA));
         }
 
         tpEqn().setReference(pRefCell, pRefValue);
@@ -69,7 +76,7 @@ else
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi -= tpEqn().flux();
+            phi = phiHbyA - tpEqn().flux();
         }
     }
 }
@@ -81,11 +88,11 @@ p.relax();
 
 if (pressureImplicitPorosity)
 {
-    U -= trTU()&fvc::grad(p);
+    U = HbyA - (trTU() & fvc::grad(p));
 }
 else
 {
-    U -= trAU()*fvc::grad(p);
+    U = HbyA - trAU()*fvc::grad(p);
 }
 
 U.correctBoundaryConditions();

applications/solvers/compressible/rhoSimplecFoam/pEqn.H

@@ -7,7 +7,10 @@ volScalarField p0(p);
 
 volScalarField AU(UEqn().A());
 volScalarField AtU(AU - UEqn().H1());
-U = UEqn().H()/AU;
+
+volVectorField HbyA("HbyA", U);
+HbyA = UEqn().H()/AU;
+
 UEqn.clear();
 
 bool closedVolume = false;
@@ -19,7 +22,7 @@ if (simple.transonic())
         surfaceScalarField phid
         (
             "phid",
-            fvc::interpolate(psi*U) & mesh.Sf()
+            fvc::interpolate(psi*HbyA) & mesh.Sf()
         );
 
         surfaceScalarField phic
@@ -56,13 +59,18 @@ else
 {
     while (simple.correctNonOrthogonal())
     {
-        phi = fvc::interpolate(rho*U) & mesh.Sf();
+        surfaceScalarField phiHbyA
+        (
+            "phiHbyA",
+            fvc::interpolate(rho*HbyA) & mesh.Sf()
+        );
+
         closedVolume = adjustPhi(phi, U, p);
         phi += fvc::interpolate(rho/AtU - rho/AU)*fvc::snGrad(p)*mesh.magSf();
 
         fvScalarMatrix pEqn
         (
-            fvc::div(phi)
+            fvc::div(phiHbyA)
         //- fvm::laplacian(rho/AU, p)
           - fvm::laplacian(rho/AtU, p)
         );
@@ -73,7 +81,7 @@ else
 
         if (simple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -85,8 +93,8 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-U -= (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU);
+U = HbyA - (fvc::grad(p0)*(1.0/AU - 1.0/AtU) + fvc::grad(p)/AtU);
-//U -= fvc::grad(p)/AU;
+//U = HbyA - fvc::grad(p)/AU;
 
 U.correctBoundaryConditions();
 

applications/solvers/compressible/sonicFoam/pEqn.H

@@ -1,14 +1,15 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 surfaceScalarField phid
 (
     "phid",
     fvc::interpolate(psi)
    *(
-        (fvc::interpolate(U) & mesh.Sf())
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     )
 );
@@ -33,5 +34,5 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();

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

@@ -1,14 +1,15 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = UEqn.H()/UEqn.A();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 surfaceScalarField phid
 (
     "phid",
     fvc::interpolate(psi)
    *(
-        (fvc::interpolate(U) & mesh.Sf())
+        (fvc::interpolate(HbyA) & mesh.Sf())
       - fvc::meshPhi(rho, U)
     )
 );
@@ -29,5 +30,5 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();

applications/solvers/electromagnetics/mhdFoam/mhdFoam.C

@@ -93,17 +93,21 @@ int main(int argc, char *argv[])
             for (int corr=0; corr<nCorr; corr++)
             {
                 volScalarField rAU(1.0/UEqn.A());
+                volVectorField HbyA("HbyA", U);
+                HbyA = rAU*UEqn.H();
 
-                U = rAU*UEqn.H();
+                surfaceScalarField phiHbyA
-
+                (
-                phi = (fvc::interpolate(U) & mesh.Sf())
+                    "phiHbyA",
-                    + fvc::ddtPhiCorr(rAU, U, phi);
+                    (fvc::interpolate(HbyA) & mesh.Sf())
+                  + fvc::ddtPhiCorr(rAU, U, phi)
+                );
 
                 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
                 {
                     fvScalarMatrix pEqn
                     (
-                        fvm::laplacian(rAU, p) == fvc::div(phi)
+                        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                     );
 
                     pEqn.setReference(pRefCell, pRefValue);
@@ -111,13 +115,13 @@ int main(int argc, char *argv[])
 
                     if (nonOrth == nNonOrthCorr)
                     {
-                        phi -= pEqn.flux();
+                        phi = phiHbyA - pEqn.flux();
                     }
                 }
 
                 #include "continuityErrs.H"
 
-                U -= rAU*fvc::grad(p);
+                U = HbyA - rAU*fvc::grad(p);
                 U.correctBoundaryConditions();
             }
         }

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H

@@ -2,19 +2,24 @@
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    phi = (fvc::interpolate(U) & mesh.Sf())
+    surfaceScalarField phig(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
-        + fvc::ddtPhiCorr(rAU, U, phi);
 
-    surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
+    surfaceScalarField phiHbyA
-    phi -= buoyancyPhi;
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU, U, phi)
+      - phig
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -24,14 +29,14 @@
         if (pimple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
             p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf);
+            U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
         }
     }

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H

@@ -2,20 +2,27 @@
     volScalarField rAU("rAU", 1.0/UEqn().A());
     surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
 
-    U = rAU*UEqn().H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn().H();
     UEqn.clear();
 
-    phi = fvc::interpolate(U) & mesh.Sf();
+    surfaceScalarField phig(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
-    adjustPhi(phi, U, p_rgh);
 
-    surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
+    surfaceScalarField phiHbyA
-    phi -= buoyancyPhi;
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+    );
+
+    adjustPhi(phiHbyA, U, p_rgh);
+
+    phiHbyA -= phig;
 
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -25,14 +32,14 @@
         if (simple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
             p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rAUf);
+            U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
         }
     }

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -8,21 +8,26 @@
     volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
-    phi = fvc::interpolate(rho)*
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+    surfaceScalarField phiHbyA
     (
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
-      + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        )
+      + phig
     );
 
-    surfaceScalarField buoyancyPhi(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-    phi += buoyancyPhi;
-
     fvScalarMatrix p_rghDDtEqn
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-      + fvc::div(phi)
+      + fvc::div(phiHbyA)
     );
 
     while (pimple.correctNonOrthogonal())
@@ -38,14 +43,14 @@
         if (pimple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
-            phi += p_rghEqn.flux();
+            phi = phiHbyA + p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
             p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U += rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
             K = 0.5*magSqr(U);
         }

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -5,20 +5,27 @@
     volScalarField rAU(1.0/UEqn().A());
     surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
-    U = rAU*UEqn().H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn().H();
     UEqn.clear();
 
-    phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
+    surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-    bool closedVolume = adjustPhi(phi, U, p_rgh);
 
-    surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+    surfaceScalarField phiHbyA
-    phi -= buoyancyPhi;
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
+    );
+
+    bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);
+
+    phiHbyA -= phig
 
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@@ -27,14 +34,14 @@
         if (simple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
             p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
+            U = HbyA - rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
         }
     }

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

@@ -8,24 +8,27 @@
     volScalarField rAU(1.0/UEqn().A());
     surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
-    U = rAU*UEqn().H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn().H();
 
-    surfaceScalarField phiU
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+    surfaceScalarField phiHbyA
     (
+        "phiHbyA",
         fvc::interpolate(rho)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
+      + phig
     );
 
-    phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
-
     {
         fvScalarMatrix p_rghDDtEqn
         (
             fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
         );
 
         // Thermodynamic density needs to be updated by psi*d(p) after the
@@ -57,7 +60,11 @@
 
             if (nonOrth == nNonOrthCorr)
             {
-                phi += p_rghEqn.flux();
+                phi = phiHbyA + p_rghEqn.flux();
+                U = HbyA
+                  + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+                U.correctBoundaryConditions();
+                K = 0.5*magSqr(U);
             }
         }
 
@@ -65,9 +72,6 @@
         thermo.rho() += psi*p_rgh;
     }
 
-    // Correct velocity field
-    U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
-    U.correctBoundaryConditions();
     p = p_rgh + rho*gh;
 
     // Update pressure time derivative

applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H

@@ -6,27 +6,23 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    U = rAU*(UEqn == sources(rho, U))().H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*(UEqn == sources(rho, U))().H();
 
-    if (pZones.size() > 0)
+    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
+    if (pZones.size() == 0)
     {
-        // ddtPhiCorr not well defined for cases with porosity
+        // ddtPhiCorr only used without porosity
-        phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
+        phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi);
-    }
-    else
-    {
-        phi =
-            fvc::interpolate(rho)
-           *(
-                (fvc::interpolate(U) & mesh.Sf())
-              + fvc::ddtPhiCorr(rAU, rho, U, phi)
-            );
     }
 
+    phiHbyA *= fvc::interpolate(rho);
+
+
     fvScalarMatrix pDDtEqn
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-      + fvc::div(phi)
+      + fvc::div(phiHbyA)
      ==
         parcels.Srho()
       + sources(psi, p, rho.name())
@@ -46,7 +42,7 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 
@@ -58,7 +54,7 @@
     #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
     #include "compressibleContinuityErrs.H"
 
-    U -= rAU*fvc::grad(p);
+    U = HbyA - rAU*fvc::grad(p);
     U.correctBoundaryConditions();
     sources.correct(U);
 

applications/solvers/lagrangian/coalChemistryFoam/pEqn.H

@@ -1,7 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*(UEqn == sources(rho, U))().H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*(UEqn == sources(rho, U))().H();
 
 if (pimple.transonic())
 {
@@ -10,7 +11,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -39,19 +40,22 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
         fvc::interpolate(rho)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
          ==
             coalParcels.Srho()
@@ -64,7 +68,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -72,7 +76,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 sources.correct(U);
 

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H

@@ -6,27 +6,23 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    U = rAU*(UEqn == sources(rho, U))().H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*(UEqn == sources(rho, U))().H();
 
-    if (pZones.size() > 0)
+    surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf());
+    if (pZones.size() == 0)
     {
-        // ddtPhiCorr not well defined for cases with porosity
+        // ddtPhiCorr only used without porosity
-        phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
+        phiHbyA += fvc::ddtPhiCorr(rAU, rho, U, phi);
-    }
-    else
-    {
-        phi =
-            fvc::interpolate(rho)
-           *(
-                (fvc::interpolate(U) & mesh.Sf())
-              + fvc::ddtPhiCorr(rAU, rho, U, phi)
-            );
     }
 
+    phiHbyA *= fvc::interpolate(rho);
+
+
     fvScalarMatrix pDDtEqn
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-      + fvc::div(phi)
+      + fvc::div(phiHbyA)
      ==
         parcels.Srho()
       + sources(psi, p, rho.name())
@@ -47,7 +43,7 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 
@@ -57,7 +53,7 @@
     #include "rhoEqn.H"
     #include "compressibleContinuityErrs.H"
 
-    U -= rAU*fvc::grad(p);
+    U = HbyA - rAU*fvc::grad(p);
     U.correctBoundaryConditions();
     sources.correct(U);
 

applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H

@@ -2,25 +2,29 @@ rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf(rAU.name() + 'f', fvc::interpolate(rho*rAU));
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
-surfaceScalarField phiU
+surfaceScalarField phig(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+surfaceScalarField phiHbyA
 (
+    "phiHbyA",
     fvc::interpolate(rho)
    *(
-        (fvc::interpolate(U) & mesh.Sf())
+        (fvc::interpolate(HbyA) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     )
+  - phig
 );
 
-phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
 while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix p_rghEqn
     (
         fvc::ddt(psi, rho)*gh
-      + fvc::div(phi)
+      + fvc::div(phiHbyA)
       + fvm::ddt(psi, p_rgh)
       - fvm::laplacian(rhorAUf, p_rgh)
      ==
@@ -32,7 +36,9 @@ while (pimple.correctNonOrthogonal())
 
     if (pimple.finalNonOrthogonalIter())
     {
-        phi += p_rghEqn.flux();
+        phi = phiHbyA + p_rghEqn.flux();
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() - phig)/rhorAUf);
+        U.correctBoundaryConditions();
     }
 }
 
@@ -41,8 +47,6 @@ p = p_rgh + rho*gh;
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U += rAU*fvc::reconstruct((phi - phiU)/rhorAUf);
-U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 
 dpdt = fvc::ddt(p);

applications/solvers/lagrangian/reactingParcelFoam/pEqn.H

@@ -1,7 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -10,7 +11,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -36,19 +37,22 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
         fvc::interpolate(rho)
        *(
-            (fvc::interpolate(U) & mesh.Sf())
+            (fvc::interpolate(HbyA) & mesh.Sf())
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
          ==
             parcels.Srho()
@@ -58,7 +62,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -66,7 +70,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 

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

@@ -1,7 +1,8 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-U = rAU*UEqn.H();
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn.H();
 
 if (pimple.transonic())
 {
@@ -10,7 +11,7 @@ if (pimple.transonic())
         "phid",
         fvc::interpolate(psi)
        *(
-            ((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
+            ((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         )
     );
@@ -36,19 +37,22 @@ if (pimple.transonic())
 }
 else
 {
-    phi =
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
         fvc::interpolate(rho)
        *(
-            ((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
+            ((fvc::interpolate(HbyA) & mesh.Sf()) - fvc::meshPhi(rho, U))
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
-        );
+        )
+    );
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvc::div(phi)
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rho*rAU, p)
          ==
             parcels.Srho()
@@ -58,7 +62,7 @@ else
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi += pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
@@ -66,7 +70,7 @@ else
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
-U -= rAU*fvc::grad(p);
+U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
 

applications/solvers/multiphase/cavitatingFoam/pEqn.H

@@ -9,11 +9,13 @@
         )/psi;
     }
 
-    surfaceScalarField rhof(fvc::interpolate(rho, "rhof"));
+    surfaceScalarField rhof("rhof", fvc::interpolate(rho));
 
     volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rAUf("rAUf", rhof*fvc::interpolate(rAU));
-    volVectorField HbyA("HbyA", rAU*UEqn.H());
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
 
     phiv = (fvc::interpolate(HbyA) & mesh.Sf())
          + fvc::ddtPhiCorr(rAU, rho, U, phiv);
@@ -22,8 +24,6 @@
 
     phiv -= phiGradp/rhof;
 
-    #include "resetPhivPatches.H"
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn

applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H

@@ -1,15 +0,0 @@
-fvsPatchScalarFieldField& phiPatches = phi.boundaryField();
-const fvPatchScalarFieldField& rhoPatches = rho.boundaryField();
-const fvPatchVectorFieldField& Upatches = U.boundaryField();
-const fvsPatchVectorFieldField& SfPatches = mesh.Sf().boundaryField();
-
-forAll(phiPatches, patchI)
-{
-    if (phi.boundaryField().types()[patchI] == "calculated")
-    {
-        calculatedFvsPatchScalarField& phiPatch =
-            refCast<calculatedFvsPatchScalarField>(phiPatches[patchI]);
-
-        phiPatch == ((rhoPatches[patchI]*Upatches[patchI]) & SfPatches[patchI]);
-    }
-}

applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H

@@ -1,19 +0,0 @@
-surfaceScalarField::GeometricBoundaryField& phivPatches =
-    phiv.boundaryField();
-
-const volVectorField::GeometricBoundaryField& Upatches =
-    U.boundaryField();
-
-const surfaceVectorField::GeometricBoundaryField& SfPatches =
-    mesh.Sf().boundaryField();
-
-forAll(phivPatches, patchI)
-{
-    if (phiv.boundaryField().types()[patchI] == "calculated")
-    {
-        calculatedFvsPatchScalarField& phivPatch =
-            refCast<calculatedFvsPatchScalarField>(phivPatches[patchI]);
-
-        phivPatch == (Upatches[patchI] & SfPatches[patchI]);
-    }
-}

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C

@@ -41,10 +41,10 @@ void Foam::conformalVoronoiMesh::conformToSurface()
 {
     reconformationMode reconfMode = reconformationControl();
 
-    if (Pstream::parRun())
+//    if (Pstream::parRun())
-    {
+//    {
-        seedProcessorBoundarySurfaces(true);
+//        seedProcessorBoundarySurfaces(true);
-    }
+//    }
 
     if (reconfMode == rmNone)
     {
@@ -74,15 +74,15 @@ void Foam::conformalVoronoiMesh::conformToSurface()
         storeSurfaceConformation();
     }
 
-    if (Pstream::parRun())
+//    if (Pstream::parRun())
-    {
+//    {
-        label nFarPoints = removeProcessorBoundarySeeds(true);
+//        label nFarPoints = removeProcessorBoundarySeeds(true);
-
+//
-        reduce(nFarPoints, sumOp<label>());
+//        reduce(nFarPoints, sumOp<label>());
-
+//
-        Info<< "    Removed " << nFarPoints
+//        Info<< "    Removed " << nFarPoints
-            << " far points from the mesh." << endl;
+//            << " far points from the mesh." << endl;
-    }
+//    }
 
     // reportSurfaceConformationQuality();
 }
@@ -915,7 +915,7 @@ void Foam::conformalVoronoiMesh::seedProcessorBoundarySurfaces
 
     label nFarPoints = 0;
 
-    const scalar normalDistance = 5.0;
+    const scalar normalDistance = 2.0;
     const scalar pert = 0.1*(rndGen_.scalar01() - 0.5);
 
     forAll(bMesh, patchI)
@@ -1465,7 +1465,7 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
     label cIInner = 0;
     label cIOuter = 0;
 
-//    seedProcessorBoundarySurfaces(true);
+    seedProcessorBoundarySurfaces(true);
 
     label cellIndexCount = 0;
     for
@@ -1530,7 +1530,7 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
         if
         (
             (testCellInfluence[cit->cellIndex()] == 0)
-         //&& (cit->real() || cit->hasFarPoint())
+         && (cit->real() || cit->hasFarPoint())
         )
         {
             const Foam::point& cc = cit->dual();
@@ -1553,15 +1553,18 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
 
     timeCheck("End of testing cell influence");
 
-    Pout<< "Number of quick rejections = " << nQuickRejections << endl;
+//    Pout<< "Number of quick rejections = " << nQuickRejections << endl;
     Pout<< "Number of influences = " << circumcentre.size() << endl;
 
     // Increasing the circumspheres to increase the overlaps and compensate for
     // floating point errors missing some referrals
     labelListList circumsphereOverlaps
-    (
+        = decomposition_().overlapsProcessors
-        overlapsProc(circumcentre, sqr(1.01)*circumradiusSqr)
+          (
-    );
+              circumcentre,
+              sqr(1.01)*circumradiusSqr,
+              false
+          );
 
     timeCheck("End of increasing overlaps");
 
@@ -1624,13 +1627,12 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
         cIOuter++;
     }
 
+    label nFarPoints = removeProcessorBoundarySeeds(true);
 
-//    label nFarPoints = removeProcessorBoundarySeeds(true);
+    reduce(nFarPoints, sumOp<label>());
-//
+
-//    reduce(nFarPoints, sumOp<label>());
+    Info<< "    Removed " << nFarPoints
-//
+        << " far points from the mesh." << endl;
-//    Info<< "    Removed " << nFarPoints
-//        << " far points from the mesh." << endl;
 
 
 //    seedProcessorBoundarySurfaces(false);
@@ -1745,6 +1747,8 @@ void Foam::conformalVoronoiMesh::referVertices
 
     timeCheck("Start of referVertices " + stageName + " insertion.");
 
+    label inserted = 0;
+
     for (label procI = 0; procI < Pstream::nProcs(); procI++)
     {
         const labelList& constructMap = pointMap.constructMap()[procI];
@@ -1774,12 +1778,18 @@ void Foam::conformalVoronoiMesh::referVertices
                         encodedProcI
                     );
 
+                    inserted++;
+
                     receivedVertices[procI].insert(origIndex);
                 }
             }
         }
     }
 
+    reduce(inserted, sumOp<label>());
+
+    Info<< "    Inserted " << stageName << " vertices " << inserted << endl;
+
     Info<< "    Total " << stageName << " vertices " << totalVertices << endl;
 
     timeCheck("End of referVertices " + stageName);

src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C

@@ -50,12 +50,13 @@ namespace Foam
     const char* Foam::NamedEnum
     <
         Foam::fieldValues::cellSource::operationType,
-        8
+        9
     >::names[] =
     {
         "none",
         "sum",
         "average",
+        "weightedAverage",
         "volAverage",
         "volIntegrate",
         "min",
@@ -68,7 +69,7 @@ namespace Foam
 const Foam::NamedEnum<Foam::fieldValues::cellSource::sourceType, 2>
     Foam::fieldValues::cellSource::sourceTypeNames_;
 
-const Foam::NamedEnum<Foam::fieldValues::cellSource::operationType, 8>
+const Foam::NamedEnum<Foam::fieldValues::cellSource::operationType, 9>
     Foam::fieldValues::cellSource::operationTypeNames_;
 
 

src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.H

@@ -51,9 +51,12 @@ Description
       - none
       - sum
       - average
+      - weightedAverage
       - volAverage
       - volIntegrate
       - CoV (Coefficient of variation: standard deviation/mean)
+      - min
+      - max
 
 SourceFiles
     cellSource.C
@@ -105,6 +108,7 @@ public:
             opNone,
             opSum,
             opAverage,
+            opWeightedAverage,
             opVolAverage,
             opVolIntegrate,
             opMin,
@@ -113,7 +117,7 @@ public:
         };
 
         //- Operation type names
-        static const NamedEnum<operationType, 8> operationTypeNames_;
+        static const NamedEnum<operationType, 9> operationTypeNames_;
 
 
 private:
@@ -169,7 +173,8 @@ protected:
         Type processValues
         (
             const Field<Type>& values,
-            const scalarField& V
+            const scalarField& V,
+            const scalarField& weightField
         ) const;
 
         //- Output file header information

src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceTemplates.C

@@ -78,7 +78,8 @@ template<class Type>
 Type Foam::fieldValues::cellSource::processValues
 (
     const Field<Type>& values,
-    const scalarField& V
+    const scalarField& V,
+    const scalarField& weightField
 ) const
 {
     Type result = pTraits<Type>::zero;
@@ -94,6 +95,11 @@ Type Foam::fieldValues::cellSource::processValues
             result = sum(values)/values.size();
             break;
         }
+        case opWeightedAverage:
+        {
+            result = sum(values)/sum(weightField);
+            break;
+        }
         case opVolAverage:
         {
             result = sum(values*V)/sum(V);
@@ -169,7 +175,7 @@ bool Foam::fieldValues::cellSource::writeValues(const word& fieldName)
 
         if (Pstream::master())
         {
-            Type result = processValues(values, V);
+            Type result = processValues(values, V, weightField);
 
             if (valueOutput_)
             {

src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C

@@ -53,12 +53,13 @@ namespace Foam
     const char* Foam::NamedEnum
     <
         Foam::fieldValues::faceSource::operationType,
-        8
+        9
     >::names[] =
     {
         "none",
         "sum",
         "average",
+        "weightedAverage",
         "areaAverage",
         "areaIntegrate",
         "min",
@@ -72,7 +73,7 @@ namespace Foam
 const Foam::NamedEnum<Foam::fieldValues::faceSource::sourceType, 3>
     Foam::fieldValues::faceSource::sourceTypeNames_;
 
-const Foam::NamedEnum<Foam::fieldValues::faceSource::operationType, 8>
+const Foam::NamedEnum<Foam::fieldValues::faceSource::operationType, 9>
     Foam::fieldValues::faceSource::operationTypeNames_;
 
 

src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H

@@ -59,6 +59,7 @@ Description
       - none
       - sum
       - average
+      - weightedAverage
       - areaAverage
       - areaIntegrate
       - min
@@ -132,6 +133,7 @@ public:
             opNone,
             opSum,
             opAverage,
+            opWeightedAverage,
             opAreaAverage,
             opAreaIntegrate,
             opMin,
@@ -140,7 +142,7 @@ public:
         };
 
         //- Operation type names
-        static const NamedEnum<operationType, 8> operationTypeNames_;
+        static const NamedEnum<operationType, 9> operationTypeNames_;
 
 
 private:
@@ -215,7 +217,8 @@ protected:
         Type processValues
         (
             const Field<Type>& values,
-            const scalarField& magSf
+            const scalarField& magSf,
+            const scalarField& weightField
         ) const;
 
         //- Output file header information

src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C

@@ -97,7 +97,9 @@ template<class Type>
 Type Foam::fieldValues::faceSource::processValues
 (
     const Field<Type>& values,
-    const scalarField& magSf
+    const scalarField& magSf,
+    const scalarField& weightField
+
 ) const
 {
     Type result = pTraits<Type>::zero;
@@ -113,6 +115,11 @@ Type Foam::fieldValues::faceSource::processValues
             result = sum(values)/values.size();
             break;
         }
+        case opWeightedAverage:
+        {
+            result = sum(values)/sum(weightField);
+            break;
+        }
         case opAreaAverage:
         {
             result = sum(values*magSf)/sum(magSf);
@@ -203,7 +210,7 @@ bool Foam::fieldValues::faceSource::writeValues(const word& fieldName)
 
         if (Pstream::master())
         {
-            Type result = processValues(values, magSf);
+            Type result = processValues(values, magSf, weightField);
 
             if (valueOutput_)
             {

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/CH4

@@ -27,22 +27,26 @@ boundaryField
         inletValue      uniform 0;
         value           uniform 0;
     }
+
     sides
     {
         type            inletOutlet;
         inletValue      uniform 0;
         value           uniform 0;
     }
+
     base
     {
         type            zeroGradient;
     }
+
     inlet
     {
         type            fixedValue;
         value           uniform 1.0;
     }
-    frontBack
+
+    frontAndBack
     {
         type            empty;
     }

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/G

@@ -28,6 +28,11 @@ boundaryField
         emissivity      uniform 1.0;
         value           uniform 0;
     }
+
+    frontAndBack
+    {
+        type            empty;
+    }
 }
 
 // ************************************************************************* //

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/IDefault

@@ -28,6 +28,11 @@ boundaryField
         emissivity      uniform 1.0;
         value           uniform 0;
     }
+
+    frontAndBack
+    {
+        type            empty;
+    }
 }
 
 // ************************************************************************* //

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/N2

@@ -25,19 +25,23 @@ boundaryField
     {
         type            calculated;
     }
+
     sides
     {
         type            calculated;
     }
+
     base
     {
         type            calculated;
     }
+
     inlet
     {
         type            calculated;
     }
-    frontBack
+
+    frontAndBack
     {
         type            empty;
     }

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/O2

@@ -27,22 +27,26 @@ boundaryField
         inletValue      $internalField;
         value           $internalField;
     }
+
     sides
     {
         type            inletOutlet;
         inletValue      $internalField;
         value           $internalField;
     }
+
     base
     {
         type            zeroGradient;
     }
+
     inlet
     {
         type            fixedValue;
         value           uniform 0;
     }
-    frontBack
+
+    frontAndBack
     {
         type            empty;
     }

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/T

@@ -27,22 +27,26 @@ boundaryField
         inletValue      uniform 300;
         value           uniform 300;
     }
+
     sides
     {
         type            inletOutlet;
         inletValue      uniform 300;
         value           uniform 300;
     }
+
     base
     {
         type            zeroGradient;
     }
+
     inlet
     {
         type            fixedValue;
         value           uniform 300;
     }
-    frontBack
+
+    frontAndBack
     {
         type            empty;
     }

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/U

@@ -29,23 +29,27 @@ boundaryField
         value           uniform (0 0 0);
 
     }
+
     sides
     {
         type            pressureInletOutletVelocity;
         outletValue     uniform (0 0 0);
         value           uniform (0 0 0);
     }
+
     base
     {
         type            fixedValue;
         value           uniform (0 0 0);
     }
+
     inlet
     {
         type            fixedValue;
         value           uniform (0 0.05 0);
     }
-    frontBack
+
+    frontAndBack
     {
         type        empty;
     }

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/Ydefault

@@ -45,7 +45,8 @@ boundaryField
         inletValue      $internalField;
         value           $internalField;
     }
-    frontBack
+
+    frontAndBack
     {
         type            empty;
     }

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/alphaSgs

@@ -25,19 +25,23 @@ boundaryField
     {
         type            zeroGradient;
     }
+
     sides
     {
         type            zeroGradient;
     }
+
     base
     {
         type            zeroGradient;
     }
+
     inlet
     {
         type            zeroGradient;
     }
-    frontBack
+
+    frontAndBack
     {
         type            empty;
     }

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/k

@@ -27,23 +27,27 @@ boundaryField
         inletValue      uniform 1e-4;
         value           uniform 1e-4;
     }
+
     sides
     {
         type            inletOutlet;
         inletValue      uniform 1e-4;
         value           uniform 1e-4;
     }
+
     base
     {
         type            fixedValue;
         value           uniform 1e-4;
     }
+
     inlet
     {
         type            fixedValue;
         value           uniform 1e-4;
     }
-    frontBack
+
+    frontAndBack
     {
         type            empty;
     }

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/muSgs

@@ -25,19 +25,23 @@ boundaryField
     {
         type            zeroGradient;
     }
+
     sides
     {
         type            zeroGradient;
     }
+
     base
     {
         type            zeroGradient;
     }
+
     inlet
     {
         type            zeroGradient;
     }
-    frontBack
+
+    frontAndBack
     {
         type            empty;
     }

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/omega

@@ -28,26 +28,29 @@ boundaryField
         inletValue      uniform 0;
         value           uniform 0;
     }
+
     sides
     {
         type            inletOutlet;
         inletValue      uniform 0;
         value           uniform 0;
     }
+
     base
     {
         type            zeroGradient;
     }
+
     inlet
     {
         type            fixedValue;
         value           uniform 0.0;
     }
-    frontBack
+
+    frontAndBack
     {
         type            empty;
     }
-
 }
 
 // ************************************************************************* //

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p

@@ -32,17 +32,20 @@ boundaryField
         type            calculated;
         value           $internalField;
     }
+
     base
     {
         type            calculated;
         value           $internalField;
     }
+
     inlet
     {
         type            calculated;
         value           $internalField;
     }
-    frontBack
+
+    frontAndBack
     {
         type            empty;
     }

tutorials/combustion/fireFoam/les/smallPoolFire2D/0/p_rgh

@@ -26,6 +26,7 @@ boundaryField
         type            buoyantPressure;
         value           $internalField;
     }
+
     sides
     {
         type            totalPressure;
@@ -37,17 +38,20 @@ boundaryField
         gamma           1;
         value           $internalField;
     }
+
     base
     {
         type            buoyantPressure;
         value           $internalField;
     }
+
     inlet
     {
         type            buoyantPressure;
         value           $internalField;
     }
-    frontBack
+
+    frontAndBack
     {
         type            empty;
     }