CorrectPhi: Added pressureReference argument to set the reference cell

so that the same reference cell is used for pcorr and p or p_rgh to improve consistency between flux and flux correction.
2021-04-30 21:07:15 +01:00
parent 007658b768
commit e8ff92cd67
72 changed files with 172 additions and 122 deletions
--- a/applications/solvers/combustion/reactingFoam/buoyantReactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/buoyantReactingFoam/createFields.H
@ -96,7 +96,6 @@ pressureReference pressureReference
 (
    p,
    p_rgh,
-    rho,
    pimple.dict(),
    thermo.incompressible()
 );
--- a/applications/solvers/combustion/reactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/createFields.H
@ -38,7 +38,7 @@ volScalarField& p = thermo.p();

 #include "compressibleCreatePhi.H"

-pressureReference pressureReference(p, rho, pimple.dict(), false);
+pressureReference pressureReference(p, pimple.dict(), false);

 mesh.setFluxRequired(p.name());

--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
@ -43,7 +43,6 @@ volVectorField U
 pressureReference pressureReference
 (
    p,
-    rho,
    pimple.dict(),
    thermo.incompressible()
 );
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@ -41,7 +41,6 @@ volVectorField U
 pressureReference pressureReference
 (
    p,
-    rho,
    simple.dict(),
    thermo.incompressible()
 );
--- a/applications/solvers/electromagnetics/mhdFoam/createFields.H
+++ b/applications/solvers/electromagnetics/mhdFoam/createFields.H
@ -109,10 +109,7 @@ DB.name() = "DB";
 dimensionedScalar DBU = 1.0/(2.0*mu*rho);
 DBU.name() = "DBU";

-
-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, piso.dict(), pRefCell, pRefValue);
+pressureReference pressureReference(p, piso.dict());

 mesh.setFluxRequired(p.name());
 mesh.setFluxRequired(pB.name());
--- a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
+++ b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
@ -2,7 +2,7 @@
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -50,6 +50,7 @@ Description

 #include "fvCFD.H"
 #include "pisoControl.H"
+#include "pressureReference.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

@ -113,7 +114,12 @@ int main(int argc, char *argv[])
                        fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
                    );

-                    pEqn.setReference(pRefCell, pRefValue);
+                    pEqn.setReference
+                    (
+                        pressureReference.refCell(),
+                        pressureReference.refValue()
+                    );
+
                    pEqn.solve();

                    if (piso.finalNonOrthogonalIter())
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
@ -83,7 +83,6 @@ pressureReference pressureReference
 (
    p,
    p_rgh,
-    rho,
    pimple.dict(),
    thermo.incompressible()
 );
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
@ -81,7 +81,6 @@ pressureReference pressureReference
 (
    p,
    p_rgh,
-    rho,
    simple.dict(),
    thermo.incompressible()
 );
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidPressureControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidPressureControls.H
@ -8,7 +8,6 @@ forAll(fluidRegions, i)
        new pressureReference
        (
            p_rghFluid[i],
-            rhoFluid[i],
            pimples.pimple(i).dict(),
            false
        )
--- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
+++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C
@ -49,6 +49,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "kinematicMomentumTransportModel.H"
 #include "simpleControl.H"
+#include "pressureReference.H"
 #include "fvModels.H"
 #include "fvConstraints.H"

@ -143,7 +144,12 @@ int main(int argc, char *argv[])
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                );

-                pEqn.setReference(pRefCell, pRefValue);
+                pEqn.setReference
+                (
+                    pressureReference.refCell(),
+                    pressureReference.refValue()
+                );
+
                pEqn.solve();

                if (simple.finalNonOrthogonalIter())
--- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/createFields.H
+++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/createFields.H
@ -29,9 +29,8 @@ volVectorField U
 #include "createPhi.H"


-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, simple.dict(), pRefCell, pRefValue);
+pressureReference pressureReference(p, simple.dict());
+
 mesh.setFluxRequired(p.name());


--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/SRFPimpleFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/SRFPimpleFoam.C
@ -36,6 +36,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "kinematicMomentumTransportModel.H"
 #include "pimpleControl.H"
+#include "pressureReference.H"
 #include "SRFModel.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/createFields.H
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/createFields.H
@ -40,9 +40,7 @@ surfaceScalarField phi
    linearInterpolate(Urel) & mesh.Sf()
 );

-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, pimple.dict(), pRefCell, pRefValue);
+pressureReference pressureReference(p, pimple.dict());
 mesh.setFluxRequired(p.name());

 Info<< "Creating SRF model\n" << endl;
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/pEqn.H
@ -38,7 +38,11 @@ while (pimple.correctNonOrthogonal())
        fvm::laplacian(rAtUrel(), p) == fvc::div(phiHbyA)
    );

-    pEqn.setReference(pRefCell, pRefValue);
+    pEqn.setReference
+    (
+        pressureReference.refCell(),
+        pressureReference.refValue()
+    );

    pEqn.solve();

--- a/applications/solvers/incompressible/pimpleFoam/correctPhi.H
+++ b/applications/solvers/incompressible/pimpleFoam/correctPhi.H
@ -11,6 +11,7 @@ CorrectPhi
    p,
    dimensionedScalar("rAUf", dimTime, 1),
    geometricZeroField(),
+    pressureReference,
    pimple
 );

--- a/applications/solvers/incompressible/pimpleFoam/createFields.H
+++ b/applications/solvers/incompressible/pimpleFoam/createFields.H
@ -30,10 +30,8 @@ volVectorField U

 #include "createPhi.H"

+pressureReference pressureReference(p, pimple.dict());

-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, pimple.dict(), pRefCell, pRefValue);
 mesh.setFluxRequired(p.name());


--- a/applications/solvers/incompressible/pimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/pimpleFoam/pEqn.H
@ -42,7 +42,11 @@ while (pimple.correctNonOrthogonal())
        fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
    );

-    pEqn.setReference(pRefCell, pRefValue);
+    pEqn.setReference
+    (
+        pressureReference.refCell(),
+        pressureReference.refValue()
+    );

    pEqn.solve();

--- a/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
+++ b/applications/solvers/incompressible/pimpleFoam/pimpleFoam.C
@ -37,6 +37,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "kinematicMomentumTransportModel.H"
 #include "pimpleControl.H"
+#include "pressureReference.H"
 #include "CorrectPhi.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
--- a/applications/solvers/incompressible/pisoFoam/createFields.H
+++ b/applications/solvers/incompressible/pisoFoam/createFields.H
@ -29,9 +29,8 @@ volVectorField U
 #include "createPhi.H"


-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, piso.dict(), pRefCell, pRefValue);
+pressureReference pressureReference(p, piso.dict());
+
 mesh.setFluxRequired(p.name());


--- a/applications/solvers/incompressible/pisoFoam/pEqn.H
+++ b/applications/solvers/incompressible/pisoFoam/pEqn.H
@ -24,7 +24,11 @@ while (piso.correctNonOrthogonal())
        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
    );

-    pEqn.setReference(pRefCell, pRefValue);
+    pEqn.setReference
+    (
+        pressureReference.refCell(),
+        pressureReference.refValue()
+    );

    pEqn.solve();

--- a/applications/solvers/incompressible/pisoFoam/pisoFoam.C
+++ b/applications/solvers/incompressible/pisoFoam/pisoFoam.C
@ -38,6 +38,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "kinematicMomentumTransportModel.H"
 #include "pisoControl.H"
+#include "pressureReference.H"
 #include "fvModels.H"
 #include "fvConstraints.H"

--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/SRFSimpleFoam.C
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/SRFSimpleFoam.C
@ -35,6 +35,7 @@ Description
 #include "kinematicMomentumTransportModel.H"
 #include "SRFModel.H"
 #include "simpleControl.H"
+#include "pressureReference.H"
 #include "fvModels.H"
 #include "fvConstraints.H"

--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/createFields.H
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/createFields.H
@ -40,9 +40,8 @@ surfaceScalarField phi
    linearInterpolate(Urel) & mesh.Sf()
 );

-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, simple.dict(), pRefCell, pRefValue);
+pressureReference pressureReference(p, simple.dict());
+
 mesh.setFluxRequired(p.name());

 Info<< "Creating SRF model\n" << endl;
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H
@ -29,7 +29,11 @@
            fvm::laplacian(rAtUrel(), p) == fvc::div(phiHbyA)
        );

-        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.setReference
+        (
+            pressureReference.refCell(),
+            pressureReference.refValue()
+        );

        pEqn.solve();

--- a/applications/solvers/incompressible/simpleFoam/createFields.H
+++ b/applications/solvers/incompressible/simpleFoam/createFields.H
@ -29,9 +29,8 @@ volVectorField U
 #include "createPhi.H"


-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, simple.dict(), pRefCell, pRefValue);
+pressureReference pressureReference(p, simple.dict());
+
 mesh.setFluxRequired(p.name());


--- a/applications/solvers/incompressible/simpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/pEqn.H
@ -28,7 +28,11 @@
            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
        );

-        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.setReference
+        (
+            pressureReference.refCell(),
+            pressureReference.refValue()
+        );

        pEqn.solve();

--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H
@ -31,7 +31,11 @@ while (simple.correctNonOrthogonal())

    fvScalarMatrix pEqn = tpEqn.ref();

-    pEqn.setReference(pRefCell, pRefValue);
+    pEqn.setReference
+    (
+        pressureReference.refCell(),
+        pressureReference.refValue()
+    );

    pEqn.solve();

--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C
@ -35,6 +35,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "kinematicMomentumTransportModel.H"
 #include "simpleControl.H"
+#include "pressureReference.H"
 #include "IOporosityModelList.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
--- a/applications/solvers/incompressible/simpleFoam/simpleFoam.C
+++ b/applications/solvers/incompressible/simpleFoam/simpleFoam.C
@ -34,6 +34,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "kinematicMomentumTransportModel.H"
 #include "simpleControl.H"
+#include "pressureReference.H"
 #include "fvModels.H"
 #include "fvConstraints.H"

--- a/applications/solvers/lagrangian/denseParticleFoam/correctPhic.H
+++ b/applications/solvers/lagrangian/denseParticleFoam/correctPhic.H
@ -10,6 +10,7 @@ CorrectPhi
    p,
    dimensionedScalar("rAUf", dimTime, 1),
    geometricZeroField(),
+    pressureReference,
    pimple
 );

--- a/applications/solvers/lagrangian/denseParticleFoam/createFields.H
+++ b/applications/solvers/lagrangian/denseParticleFoam/createFields.H
@ -59,9 +59,8 @@ surfaceScalarField phic
    linearInterpolate(Uc) & mesh.Sf()
 );

-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell(p, pimple.dict(), pRefCell, pRefValue);
+pressureReference pressureReference(p, pimple.dict());
+
 mesh.setFluxRequired(p.name());

 Info<< "Creating turbulence model\n" << endl;
--- a/applications/solvers/lagrangian/denseParticleFoam/denseParticleFoam.C
+++ b/applications/solvers/lagrangian/denseParticleFoam/denseParticleFoam.C
@ -80,6 +80,7 @@ namespace Foam
 #include "singlePhaseTransportModel.H"
 #include "phaseKinematicMomentumTransportModel.H"
 #include "pimpleControl.H"
+#include "pressureReference.H"
 #include "CorrectPhi.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
--- a/applications/solvers/lagrangian/denseParticleFoam/pEqn.H
+++ b/applications/solvers/lagrangian/denseParticleFoam/pEqn.H
@ -34,7 +34,11 @@
          + fvc::div(alphacf*phiHbyASp)
        );

-        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.setReference
+        (
+            pressureReference.refCell(),
+            pressureReference.refValue()
+        );

        pEqn.solve();

--- a/applications/solvers/lagrangian/reactingParticleFoam/buoyantReactingParticleFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParticleFoam/buoyantReactingParticleFoam/createFields.H
@ -109,7 +109,6 @@ pressureReference pressureReference
 (
    p,
    p_rgh,
-    rho,
    pimple.dict(),
    thermo.incompressible()
 );
--- a/applications/solvers/lagrangian/reactingParticleFoam/createFields.H
+++ b/applications/solvers/lagrangian/reactingParticleFoam/createFields.H
@ -90,7 +90,6 @@ volScalarField K("K", 0.5*magSqr(U));
 pressureReference pressureReference
 (
    p,
-    rho,
    pimple.dict(),
    thermo.incompressible()
 );
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C
@ -50,6 +50,7 @@ Description
 #include "compressibleInterPhaseTransportModel.H"
 #include "noPhaseChange.H"
 #include "pimpleControl.H"
+#include "pressureReference.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
 #include "CorrectPhi.H"
--- a/applications/solvers/multiphase/compressibleInterFoam/correctPhi.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/correctPhi.H
@ -11,6 +11,7 @@ CorrectPhi
    p_rgh,
    surfaceScalarField("rAUf", fvc::interpolate(rAU())),
    divU(),
+    pressureReference,
    pimple
 );

--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@ -57,6 +57,9 @@ mesh.setFluxRequired(alpha1.name());
 #include "gh.H"


+pressureReference pressureReference(mixture.p(), p_rgh, pimple.dict(), false);
+
+
 // Mass flux
 // Initialisation does not matter because rhoPhi is reset after the
 // alpha1 solution before it is used in the U equation.
--- a/applications/solvers/multiphase/driftFluxFoam/createFields.H
+++ b/applications/solvers/multiphase/driftFluxFoam/createFields.H
@ -99,16 +99,7 @@ volScalarField p
    p_rgh + rho*gh
 );

-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell
-(
-    p,
-    p_rgh,
-    pimple.dict(),
-    pRefCell,
-    pRefValue
-);
+pressureReference pressureReference(p, p_rgh, pimple.dict());

 if (p_rgh.needReference())
 {
@ -116,7 +107,8 @@ if (p_rgh.needReference())
    (
        "p",
        p.dimensions(),
-        pRefValue - getRefCellValue(p, pRefCell)
+        pressureReference.refValue()
+      - getRefCellValue(p, pressureReference.refCell())
    );
    p_rgh = p - rho*gh;
 }
--- a/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
+++ b/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
@ -41,6 +41,7 @@ Description
 #include "momentumTransportModel.H"
 #include "dynamicMomentumTransportModel.H"
 #include "pimpleControl.H"
+#include "pressureReference.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
 #include "gaussLaplacianScheme.H"
--- a/applications/solvers/multiphase/driftFluxFoam/pEqn.H
+++ b/applications/solvers/multiphase/driftFluxFoam/pEqn.H
@ -30,7 +30,11 @@
            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        p_rghEqn.setReference
+        (
+            pressureReference.refCell(),
+            getRefCellValue(p_rgh, pressureReference.refCell())
+        );

        p_rghEqn.solve();

@ -56,7 +60,8 @@
        (
            "p",
            p.dimensions(),
-            pRefValue - getRefCellValue(p, pRefCell)
+            pressureReference.refValue()
+          - getRefCellValue(p, pressureReference.refCell())
        );
        p_rgh = p - rho*gh;
    }
--- a/applications/solvers/multiphase/interFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interFoam/correctPhi.H
@ -13,6 +13,7 @@ if (divU.valid())
        p_rgh,
        surfaceScalarField("rAUf", fvc::interpolate(rAU())),
        divU(),
+        pressureReference,
        pimple
    );
 }
@ -25,6 +26,7 @@ else
        p_rgh,
        surfaceScalarField("rAUf", fvc::interpolate(rAU())),
        geometricZeroField(),
+        pressureReference,
        pimple
    );
 }
--- a/applications/solvers/multiphase/interFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/createFields.H
@ -86,16 +86,7 @@ volScalarField p
    p_rgh + rho*gh
 );

-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell
-(
-    p,
-    p_rgh,
-    pimple.dict(),
-    pRefCell,
-    pRefValue
-);
+pressureReference pressureReference(p, p_rgh, pimple.dict());

 if (p_rgh.needReference())
 {
@ -103,7 +94,8 @@ if (p_rgh.needReference())
    (
        "p",
        p.dimensions(),
-        pRefValue - getRefCellValue(p, pRefCell)
+        pressureReference.refValue()
+      - getRefCellValue(p, pressureReference.refCell())
    );
    p_rgh = p - rho*gh;
 }
--- a/applications/solvers/multiphase/interFoam/initCorrectPhi.H
+++ b/applications/solvers/multiphase/interFoam/initCorrectPhi.H
@ -29,6 +29,7 @@ if
            p_rgh,
            surfaceScalarField("rAUf", fvc::interpolate(rAU())),
            geometricZeroField(),
+            pressureReference,
            pimple
        );
    }
@ -43,6 +44,7 @@ if
            p_rgh,
            dimensionedScalar(dimTime/rho.dimensions(), 1),
            geometricZeroField(),
+            pressureReference,
            pimple
        );
    }
--- a/applications/solvers/multiphase/interFoam/interFoam.C
+++ b/applications/solvers/multiphase/interFoam/interFoam.C
@ -43,6 +43,7 @@ Description
 #include "noPhaseChange.H"
 #include "kinematicMomentumTransportModel.H"
 #include "pimpleControl.H"
+#include "pressureReference.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
 #include "CorrectPhi.H"
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/correctPhi.H
@ -11,6 +11,7 @@ CorrectPhi
    p_rgh,
    surfaceScalarField("rAUf", fvc::interpolate(rAU())),
    geometricZeroField(),
+    pressureReference,
    pimple
 );

--- a/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H
@ -99,16 +99,7 @@ volScalarField p
    p_rgh + rho*gh
 );

-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell
-(
-    p,
-    p_rgh,
-    pimple.dict(),
-    pRefCell,
-    pRefValue
-);
+pressureReference pressureReference(p, p_rgh, pimple.dict());

 if (p_rgh.needReference())
 {
@ -116,7 +107,8 @@ if (p_rgh.needReference())
    (
        "p",
        p.dimensions(),
-        pRefValue - getRefCellValue(p, pRefCell)
+        pressureReference.refValue()
+      - getRefCellValue(p, pressureReference.refCell())
    );
    p_rgh = p - rho*gh;
 }
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/initCorrectPhi.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/initCorrectPhi.H
@ -25,6 +25,7 @@ if (correctPhi)
        p_rgh,
        surfaceScalarField("rAUf", fvc::interpolate(rAU())),
        geometricZeroField(),
+        pressureReference,
        pimple
    );

@ -41,6 +42,7 @@ else
        p_rgh,
        dimensionedScalar(dimTime/rho.dimensions(), 1),
        geometricZeroField(),
+        pressureReference,
        pimple
    );

--- a/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C
@ -39,6 +39,7 @@ Description
 #include "immiscibleIncompressibleThreePhaseMixture.H"
 #include "kinematicMomentumTransportModel.H"
 #include "pimpleControl.H"
+#include "pressureReference.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
 #include "CorrectPhi.H"
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/pEqn.H
@ -45,7 +45,11 @@
            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        p_rghEqn.setReference
+        (
+            pressureReference.refCell(),
+            getRefCellValue(p_rgh, pressureReference.refCell())
+        );

        p_rghEqn.solve();

@ -77,7 +81,8 @@
        (
            "p",
            p.dimensions(),
-            pRefValue - getRefCellValue(p, pRefCell)
+            pressureReference.refValue()
+          - getRefCellValue(p, pressureReference.refCell())
        );
        p_rgh = p - rho*gh;
    }
--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@ -49,7 +49,11 @@
         == Sp_rgh
        );

-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        p_rghEqn.setReference
+        (
+            pressureReference.refCell(),
+            getRefCellValue(p_rgh, pressureReference.refCell())
+        );

        p_rghEqn.solve();

@ -81,7 +85,8 @@
        (
            "p",
            p.dimensions(),
-            pRefValue - getRefCellValue(p, pRefCell)
+            pressureReference.refValue()
+          - getRefCellValue(p, pressureReference.refCell())
        );
        p_rgh = p - rho*gh;
    }
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/createFields.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/createFields.H
@ -40,7 +40,6 @@ pressureReference pressureReference
 (
    p,
    p_rgh,
-    fluid.rho(),
    pimple.dict(),
    fluid.incompressible()
 );
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/multiphaseEulerFoam.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/multiphaseEulerFoam.C
@ -158,7 +158,13 @@ int main(int argc, char *argv[])

                        if (correctPhi)
                        {
-                            fluid.correctPhi(p_rgh, divU, pimple);
+                            fluid.correctPhi
+                            (
+                                p_rgh,
+                                divU,
+                                pressureReference,
+                                pimple
+                            );
                        }

                        if (checkMeshCourantNo)
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
@ -40,6 +40,7 @@ License
 #include "BlendedInterfacialModel.H"
 #include "movingWallVelocityFvPatchVectorField.H"
 #include "pimpleControl.H"
+#include "pressureReference.H"

 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

@ -900,6 +901,7 @@ void Foam::phaseSystem::correctPhi
 (
    const volScalarField& p_rgh,
    const tmp<volScalarField>& divU,
+    const pressureReference& pressureReference,
    nonOrthogonalSolutionControl& pimple
 )
 {
@ -955,6 +957,7 @@ void Foam::phaseSystem::correctPhi
            // surfaceScalarField("rAUf", fvc::interpolate(rAU())),
            dimensionedScalar(dimTime/dimDensity, 1),
            divU(),
+            pressureReference,
            pimple
        );

--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
@ -61,6 +61,7 @@ class blendingMethod;
 template<class modelType> class BlendedInterfacialModel;
 class surfaceTensionModel;
 class aspectRatioModel;
+class pressureReference;
 class nonOrthogonalSolutionControl;

 /*---------------------------------------------------------------------------*\
@ -700,6 +701,7 @@ public:
            (
                const volScalarField& p_rgh,
                const tmp<volScalarField>& divU,
+                const pressureReference& pressureReference,
                nonOrthogonalSolutionControl& pimple
            );

--- a/applications/solvers/multiphase/multiphaseInterFoam/correctPhi.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/correctPhi.H
@ -11,6 +11,7 @@ CorrectPhi
    p_rgh,
    surfaceScalarField("rAUf", fvc::interpolate(rAU())),
    geometricZeroField(),
+    pressureReference,
    pimple
 );

--- a/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/createFields.H
@ -69,16 +69,7 @@ volScalarField p
    p_rgh + rho*gh
 );

-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell
-(
-    p,
-    p_rgh,
-    pimple.dict(),
-    pRefCell,
-    pRefValue
-);
+pressureReference pressureReference(p, p_rgh, pimple.dict());

 if (p_rgh.needReference())
 {
@ -86,8 +77,10 @@ if (p_rgh.needReference())
    (
        "p",
        p.dimensions(),
-        pRefValue - getRefCellValue(p, pRefCell)
+        pressureReference.refValue()
+      - getRefCellValue(p, pressureReference.refCell())
    );
+    p_rgh = p - rho*gh;
 }

 mesh.setFluxRequired(p_rgh.name());
--- a/applications/solvers/multiphase/multiphaseInterFoam/initCorrectPhi.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/initCorrectPhi.H
@ -25,6 +25,7 @@ if (correctPhi)
        p_rgh,
        surfaceScalarField("rAUf", fvc::interpolate(rAU())),
        geometricZeroField(),
+        pressureReference,
        pimple
    );

@ -41,6 +42,7 @@ else
        p_rgh,
        dimensionedScalar(dimTime/rho.dimensions(), 1),
        geometricZeroField(),
+        pressureReference,
        pimple
    );

--- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
+++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C
@ -38,6 +38,7 @@ Description
 #include "multiphaseMixture.H"
 #include "kinematicMomentumTransportModel.H"
 #include "pimpleControl.H"
+#include "pressureReference.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
 #include "CorrectPhi.H"
--- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
@ -45,7 +45,11 @@
            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        p_rghEqn.setReference
+        (
+            pressureReference.refCell(),
+            getRefCellValue(p_rgh, pressureReference.refCell())
+        );

        p_rghEqn.solve();

@ -77,7 +81,8 @@
        (
            "p",
            p.dimensions(),
-            pRefValue - getRefCellValue(p, pRefCell)
+            pressureReference.refValue()
+          - getRefCellValue(p, pressureReference.refCell())
        );
        p_rgh = p - rho*gh;
    }
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/correctPhi.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/correctPhi.H
@ -11,6 +11,7 @@ CorrectPhi
    p_gh,
    surfaceScalarField("rAUf", fvc::interpolate(rAU)),
    geometricZeroField(),
+    pressureReference,
    pimple
 );

--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/createFields.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/createFields.H
@ -70,10 +70,8 @@ volScalarField p_gh
 // Height is made relative to field 'refLevel'
 p_gh =  p - (g & mesh.C());

+pressureReference pressureReference(p_gh, pimple.dict());

-label p_ghRefCell = 0;
-scalar p_ghRefValue = 0.0;
-setRefCell(p_gh, pimple.dict(), p_ghRefCell, p_ghRefValue);
 mesh.setFluxRequired(p_gh.name());

 #include "createMRF.H"
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/initCorrectPhi.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/initCorrectPhi.H
@ -9,6 +9,7 @@ if (!runTime.restart() &&correctPhi)
        p_gh,
        surfaceScalarField("rAUf", fvc::interpolate(rAU)),
        geometricZeroField(),
+        pressureReference,
        pimple
    );
 }
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/pEqn.H
@ -34,7 +34,11 @@ while (pimple.correctNonOrthogonal())
        fvm::laplacian(rAUf, p_gh) == fvc::div(phiHbyA)
    );

-    p_ghEqn.setReference(p_ghRefCell, p_ghRefValue);
+    p_ghEqn.setReference
+    (
+        pressureReference.refCell(),
+        pressureReference.refValue()
+    );

    p_ghEqn.solve();

--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceFoam.C
@ -43,6 +43,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "kinematicMomentumTransportModel.H"
 #include "pimpleControl.H"
+#include "pressureReference.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
 #include "CorrectPhi.H"
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H
@ -88,16 +88,7 @@ volScalarField p
    p_rgh + rho*gh
 );

-label pRefCell = 0;
-scalar pRefValue = 0.0;
-setRefCell
-(
-    p,
-    p_rgh,
-    pimple.dict(),
-    pRefCell,
-    pRefValue
-);
+pressureReference pressureReference(p, p_rgh, pimple.dict());

 if (p_rgh.needReference())
 {
@ -105,7 +96,8 @@ if (p_rgh.needReference())
    (
        "p",
        p.dimensions(),
-        pRefValue - getRefCellValue(p, pRefCell)
+        pressureReference.refValue()
+      - getRefCellValue(p, pressureReference.refCell())
    );
    p_rgh = p - rho*gh;
 }
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
@ -27,7 +27,11 @@
            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        p_rghEqn.setReference
+        (
+            pressureReference.refCell(),
+            getRefCellValue(p_rgh, pressureReference.refCell())
+        );

        p_rghEqn.solve();

@ -50,7 +54,8 @@
        (
            "p",
            p.dimensions(),
-            pRefValue - getRefCellValue(p, pRefCell)
+            pressureReference.refValue()
+          - getRefCellValue(p, pressureReference.refCell())
        );
        p_rgh = p - rho*gh;
    }
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C
@ -2,7 +2,7 @@
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
     \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@ -37,6 +37,7 @@ Description
 #include "incompressibleTwoPhaseMixture.H"
 #include "kinematicMomentumTransportModel.H"
 #include "pimpleControl.H"
+#include "pressureReference.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

--- a/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C
+++ b/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C
@ -34,6 +34,7 @@ License
 #include "zeroGradientFvPatchFields.H"
 #include "adjustPhi.H"
 #include "fvcMeshPhi.H"
+#include "pressureReference.H"
 #include "nonOrthogonalSolutionControl.H"

 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -46,6 +47,7 @@ void Foam::CorrectPhi
    const volScalarField& p,
    const RAUfType& rAUf,
    const DivUType& divU,
+    const pressureReference& pressureReference,
    nonOrthogonalSolutionControl& pcorrControl
 )
 {
@ -99,7 +101,7 @@ void Foam::CorrectPhi
            fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divU
        );

-        pcorrEqn.setReference(0, 0);
+        pcorrEqn.setReference(pressureReference.refCell(), 0);

        pcorrEqn.solve();

--- a/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.H
+++ b/src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.H
@ -48,6 +48,7 @@ SourceFiles

 namespace Foam
 {
+    class pressureReference;
    class nonOrthogonalSolutionControl;

    //- If the mesh is moving correct the velocity BCs on the moving walls to
@ -77,6 +78,7 @@ namespace Foam
        const volScalarField& p,
        const RAUfType& rAUf,
        const DivUType& divU,
+        const pressureReference& pressureReference,
        nonOrthogonalSolutionControl& pcorrControl
    );

--- a/src/finiteVolume/cfdTools/general/pressureReference/pressureReference.C
+++ b/src/finiteVolume/cfdTools/general/pressureReference/pressureReference.C
@ -32,7 +32,6 @@ Foam::pressureReference::pressureReference
 (
    const volScalarField& p,
    const volScalarField& pRef,
-    const volScalarField& rho,
    const dictionary& dict,
    const bool pRefRequired
 )
@ -112,12 +111,11 @@ Foam::pressureReference::pressureReference
 Foam::pressureReference::pressureReference
 (
    const volScalarField& p,
-    const volScalarField& rho,
    const dictionary& dict,
    const bool pRefRequired
 )
 :
-    pressureReference(p, p, rho, dict, pRefRequired)
+    pressureReference(p, p, dict, pRefRequired)
 {}


--- a/src/finiteVolume/cfdTools/general/pressureReference/pressureReference.H
+++ b/src/finiteVolume/cfdTools/general/pressureReference/pressureReference.H
@ -69,7 +69,6 @@ public:
        (
            const volScalarField& p,
            const volScalarField& pRef,
-            const volScalarField& rho,
            const dictionary& dict,
            const bool pRefRequired = true
        );
@ -78,7 +77,6 @@ public:
        pressureReference
        (
            const volScalarField& p,
-            const volScalarField& rho,
            const dictionary& dict,
            const bool pRefRequired = true
        );