From e68480913f2c2792ed40d4662ada7e099c090555 Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Thu, 4 Jun 2009 19:11:12 +0100
Subject: [PATCH] pd => p

---
 .../compressibleInterDyMFoam/UEqn.H           |  6 +--
 .../compressibleInterDyMFoam/correctPhi.H     |  2 +-
 .../compressibleInterDyMFoam/createFields.H   | 30 +++---------
 .../compressibleInterDyMFoam/pEqn.H           | 43 ++++++++---------
 .../multiphase/compressibleInterFoam/UEqn.H   |  6 +--
 .../compressibleInterFoam/createFields.H      | 29 ++----------
 .../multiphase/compressibleInterFoam/pEqn.H   | 36 +++++++--------
 .../multiphase/interDyMFoam/correctPhi.H      |  4 +-
 .../multiphase/interDyMFoam/createFields.H    | 46 ++++---------------
 .../multiphase/interDyMFoam/interDyMFoam.C    | 12 -----
 .../solvers/multiphase/interDyMFoam/pEqn.H    | 22 ++++-----
 .../multiphase/interPhaseChangeFoam/UEqn.H    |  6 +--
 .../interPhaseChangeFoam/correctPhi.H         | 10 ++--
 .../interPhaseChangeFoam/createFields.H       | 29 +++---------
 .../multiphase/interPhaseChangeFoam/pEqn.H    | 24 +++++-----
 .../multiphase/multiphaseInterFoam/pEqn.H     |  4 +-
 .../solvers/multiphase/settlingFoam/UEqn.H    |  6 +--
 .../multiphase/settlingFoam/createFields.H    |  3 --
 .../solvers/multiphase/settlingFoam/pEqn.H    |  2 +-
 19 files changed, 102 insertions(+), 218 deletions(-)

diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H b/applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H
index 1576f6ba787..138e58fc7f7 100644
--- a/applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H
+++ b/applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H
@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
+                  - fvc::snGrad(p)
                 ) * mesh.magSf()
             )
         );
diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H b/applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H
index 41521057846..d82a03edb59 100644
--- a/applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H
@@ -12,7 +12,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
         pcorrTypes
     );
 
diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H b/applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H
index dc04fb45468..6b8e67cc492 100644
--- a/applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H
@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -88,24 +88,6 @@
 
     dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
 
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        )
-    );
-
     volScalarField rho1 = rho10 + psi1*p;
     volScalarField rho2 = rho20 + psi2*p;
 
@@ -152,11 +134,11 @@
     );
 
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes(p.boundaryField().types());
 
-    for (label i=0; i<pd.boundaryField().size(); i++)
+    for (label i=0; i<p.boundaryField().size(); i++)
     {
-        if (pd.boundaryField()[i].fixesValue())
+        if (p.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }
diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H
index 7e4b37061fd..e6004eb9de9 100644
--- a/applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H
@@ -2,17 +2,17 @@
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rUAf = fvc::interpolate(rUA);
 
-    tmp<fvScalarMatrix> pdEqnComp;
+    tmp<fvScalarMatrix> pEqnComp;
 
     if (transonic)
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
     }
     else
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
     }
 
 
@@ -26,16 +26,16 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-          - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+          + fvc::interpolate(rho)*(g & mesh.Sf())
+        )*rUAf;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqnIncomp
+        fvScalarMatrix pEqnIncomp
         (
             fvc::div(phi)
-          - fvm::laplacian(rUAf, pd)
+          - fvm::laplacian(rUAf, p)
         );
 
         if
@@ -51,9 +51,9 @@
                     max(alpha1, scalar(0))*(psi1/rho1)
                   + max(alpha2, scalar(0))*(psi2/rho2)
                 )
-               *pdEqnComp()
-              + pdEqnIncomp,
-                mesh.solver(pd.name() + "Final")
+               *pEqnComp()
+              + pEqnIncomp,
+                mesh.solver(p.name() + "Final")
             );
         }
         else
@@ -64,8 +64,8 @@
                     max(alpha1, scalar(0))*(psi1/rho1)
                   + max(alpha2, scalar(0))*(psi2/rho2)
                 )
-               *pdEqnComp()
-              + pdEqnIncomp
+               *pEqnComp()
+              + pEqnIncomp
             );
         }
 
@@ -73,26 +73,21 @@
         {
             dgdt =
                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
-               *(pdEqnComp & pd);
-            phi += pdEqnIncomp.flux();
+               *(pEqnComp & p);
+            phi += pEqnIncomp.flux();
         }
     }
 
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 
-    p = max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        );
+    p.max(pMin);
 
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(pd) " << min(pd).value() << endl;
+    Info<< "min(p) " << min(p).value() << endl;
 
     // Make the fluxes relative to the mesh motion
     fvc::makeRelative(phi, U);
diff --git a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
index 528e0aaafd8..0b1a9ac029d 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
+                  - fvc::snGrad(p)
                 ) * mesh.magSf()
             )
         );
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index 1f579d245bd..3e6904d383e 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -46,11 +46,6 @@
     #include "createPhi.H"
 
 
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi);
 
@@ -88,24 +83,6 @@
 
     dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
 
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        )
-    );
-
     volScalarField rho1 = rho10 + psi1*p;
     volScalarField rho2 = rho20 + psi2*p;
 
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index ebf24498ade..9d2dc23916b 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -2,17 +2,17 @@
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rUAf = fvc::interpolate(rUA);
 
-    tmp<fvScalarMatrix> pdEqnComp;
+    tmp<fvScalarMatrix> pEqnComp;
 
     if (transonic)
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
     }
     else
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
     }
 
 
@@ -26,16 +26,16 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-          - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+          + fvc::interpolate(rho)*(g & mesh.Sf())
+        )*rUAf;
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqnIncomp
+        fvScalarMatrix pEqnIncomp
         (
             fvc::div(phi)
-          - fvm::laplacian(rUAf, pd)
+          - fvm::laplacian(rUAf, p)
         );
 
         solve
@@ -44,31 +44,27 @@
                 max(alpha1, scalar(0))*(psi1/rho1)
               + max(alpha2, scalar(0))*(psi2/rho2)
             )
-           *pdEqnComp()
-          + pdEqnIncomp
+           *pEqnComp()
+          + pEqnIncomp
         );
 
         if (nonOrth == nNonOrthCorr)
         {
             dgdt =
                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
-               *(pdEqnComp & pd);
-            phi += pdEqnIncomp.flux();
+               *(pEqnComp & p);
+            phi += pEqnIncomp.flux();
         }
     }
 
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 
-    p = max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        );
+    p.max(pMin);
 
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(pd) " << min(pd).value() << endl;
+    Info<< "min(p) " << min(p).value() << endl;
 }
diff --git a/applications/solvers/multiphase/interDyMFoam/correctPhi.H b/applications/solvers/multiphase/interDyMFoam/correctPhi.H
index c975c9b3741..1f7845f347d 100644
--- a/applications/solvers/multiphase/interDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interDyMFoam/correctPhi.H
@@ -12,7 +12,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
         pcorrTypes
     );
 
@@ -27,7 +27,7 @@
             fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
         );
 
-        pcorrEqn.setReference(pdRefCell, pdRefValue);
+        pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
         if (nonOrth == nNonOrthCorr)
diff --git a/applications/solvers/multiphase/interDyMFoam/createFields.H b/applications/solvers/multiphase/interDyMFoam/createFields.H
index bcceb9d7481..26644e17f47 100644
--- a/applications/solvers/multiphase/interDyMFoam/createFields.H
+++ b/applications/solvers/multiphase/interDyMFoam/createFields.H
@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -92,47 +92,17 @@
         incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
     );
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes(p.boundaryField().types());
 
-    for (label i=0; i<pd.boundaryField().size(); i++)
+    for (label i=0; i<p.boundaryField().size(); i++)
     {
-        if (pd.boundaryField()[i].fixesValue())
+        if (p.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }
     }
 
 
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rho*(g & mesh.C())
-    );
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
-
+    label pRefCell = 0;
     scalar pRefValue = 0.0;
-
-    if (pd.needReference())
-    {
-        pRefValue = readScalar
-        (
-            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
-        );
-
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, pdRefCell)
-        );
-    }
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
diff --git a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
index 046503a8c45..cf8ec8eb323 100644
--- a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
@@ -114,18 +114,6 @@ int main(int argc, char *argv[])
             #include "pEqn.H"
         }
 
-        p = pd + rho*gh;
-
-        if (pd.needReference())
-        {
-            p += dimensionedScalar
-            (
-                "p",
-                p.dimensions(),
-                pRefValue - getRefCellValue(p, pdRefCell)
-            );
-        }
-
         turbulence->correct();
 
         runTime.write();
diff --git a/applications/solvers/multiphase/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interDyMFoam/pEqn.H
index f8ff18f6063..8825661c830 100644
--- a/applications/solvers/multiphase/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interDyMFoam/pEqn.H
@@ -7,38 +7,38 @@
 
     phi = phiU +
     (
-        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-      - ghf*fvc::snGrad(rho)
-    )*rAUf*mesh.magSf();
+        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+      + fvc::interpolate(rho)*(g & mesh.Sf())
+    )*rAUf;
 
-    if (pd.needReference())
+    if (p.needReference())
     {
         fvc::makeRelative(phi, U);
-        adjustPhi(phi, U, pd);
+        adjustPhi(phi, U, p);
         fvc::makeAbsolute(phi, U);
     }
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::laplacian(rAUf, pd) == fvc::div(phi)
+            fvm::laplacian(rAUf, p) == fvc::div(phi)
         );
 
-        pdEqn.setReference(pdRefCell, pdRefValue);
+        pEqn.setReference(pRefCell, pRefValue);
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi -= pdEqn.flux();
+            phi -= pEqn.flux();
         }
     }
 
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
index c59137c7b87..757e37a2cfd 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H
@@ -25,10 +25,10 @@
          ==
             fvc::reconstruct
             (
-                (
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
                     fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
+                  - fvc::snGrad(p)
                 ) * mesh.magSf()
             )
         );
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H b/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H
index 171e1670f47..642aa1c5d2a 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H
@@ -1,11 +1,11 @@
 {
 #   include "continuityErrs.H"
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes(p.boundaryField().types());
 
-    for (label i=0; i<pd.boundaryField().size(); i++)
+    for (label i=0; i<p.boundaryField().size(); i++)
     {
-        if (pd.boundaryField()[i].fixesValue())
+        if (p.boundaryField()[i].fixesValue())
         {
             pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
         }
@@ -22,7 +22,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("pcorr", pd.dimensions(), 0.0),
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
         pcorrTypes
     );
 
@@ -37,7 +37,7 @@
             fvm::laplacian(rUAf, pcorr) == fvc::div(phi)
         );
 
-        pcorrEqn.setReference(pdRefCell, pdRefValue);
+        pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
         if (nonOrth == nNonOrthCorr)
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
index 5e727dea3ec..98d166214b2 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -66,26 +66,9 @@
     rho.oldTime();
 
 
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
-
-    Info<< "Calculating field g.h" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rho*gh
-    );
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
 
 
     // Construct interface from alpha1 distribution
diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
index 0037d71cf52..4f290157f64 100644
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@@ -13,11 +13,11 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-          - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+          + fvc::interpolate(rho)*(g & mesh.Sf())
+        )*rUAf;
 
-    adjustPhi(phi, U, pd);
+    adjustPhi(phi, U, p);
 
     Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
     const volScalarField& vDotcP = vDotP[0]();
@@ -25,31 +25,29 @@
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvc::div(phi) - fvm::laplacian(rUAf, pd)
-          + (vDotvP - vDotcP)*(rho*gh - pSat) + fvm::Sp(vDotvP - vDotcP, pd)
+            fvc::div(phi) - fvm::laplacian(rUAf, p)
+          - (vDotvP - vDotcP)*pSat + fvm::Sp(vDotvP - vDotcP, p)
         );
 
-        pdEqn.setReference(pdRefCell, pdRefValue);
+        pEqn.setReference(pRefCell, pRefValue);
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi += pdEqn.flux();
+            phi += pEqn.flux();
         }
     }
 
-    p = pd + rho*gh;
-
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 }
diff --git a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
index 6efcdbc8b1f..9d7e36dc95a 100644
--- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
@@ -12,9 +12,9 @@
 
     phi = phiU +
         (
-            mixture.surfaceTensionForce()
+            mixture.surfaceTensionForce()*mesh.magSf()
           + fvc::interpolate(rho)*(g & mesh.Sf())
-        )*rUAf*mesh.magSf();
+        )*rUAf;
 
     adjustPhi(phi, U, p);
 
diff --git a/applications/solvers/multiphase/settlingFoam/UEqn.H b/applications/solvers/multiphase/settlingFoam/UEqn.H
index ac4c18ab19d..04e9194363d 100644
--- a/applications/solvers/multiphase/settlingFoam/UEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/UEqn.H
@@ -22,10 +22,8 @@
           ==
             fvc::reconstruct
             (
-                (
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(p)
-                )*mesh.magSf()
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              - fvc::snGrad(p)*mesh.magSf()
             )
         );
     }
diff --git a/applications/solvers/multiphase/settlingFoam/createFields.H b/applications/solvers/multiphase/settlingFoam/createFields.H
index b13649c3124..7c8445a71d4 100644
--- a/applications/solvers/multiphase/settlingFoam/createFields.H
+++ b/applications/solvers/multiphase/settlingFoam/createFields.H
@@ -337,6 +337,3 @@
         ),
         mut + mul
     );
-
-    Info<< "Calculating field (g.h)f\n" << endl;
-    surfaceScalarField ghf = surfaceScalarField("ghf", g & mesh.Cf());
diff --git a/applications/solvers/multiphase/settlingFoam/pEqn.H b/applications/solvers/multiphase/settlingFoam/pEqn.H
index 65b954c6da8..b4fc05b32ab 100644
--- a/applications/solvers/multiphase/settlingFoam/pEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/pEqn.H
@@ -15,7 +15,7 @@ phi =
     );
 
 surfaceScalarField phiU("phiU", phi);
-phi -= ghf*fvc::snGrad(rho)*rUAf*mesh.magSf();
+phi += fvc::interpolate(rho)*(g & mesh.Sf())*rUAf;
 
 for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
-- 
GitLab