//SIMPLE correction 

volScalarField DU = 1.0/UEqn.A();
surfaceScalarField DUf("DUf",linearInterpolate(DU));


//- Grabbing necessary data for Rhie-Chow interplation

dimensionedScalar deltaT ("deltaT",dimensionSet(0,0,1,0,0,0,0),scalar(runTime.deltaTValue()));
surfaceScalarField DTf("DTf",DUf/deltaT);

const surfaceVectorField ed = mesh.delta()()/mag(mesh.delta()());
Foam::fv::orthogonalSnGrad<scalar> faceGradient(mesh);

surfaceVectorField gradp_avg_f = linearInterpolate(fvc::grad(p));
surfaceVectorField gradp_f = gradp_avg_f - (gradp_avg_f & ed)*ed + (faceGradient.snGrad(p))*ed;

surfaceVectorField U_avg_f = linearInterpolate(U);
surfaceVectorField U_avg_prevIter_f = linearInterpolate(U.prevIter());
surfaceVectorField U_old_f = linearInterpolate(U.oldTime());

surfaceScalarField mdotf_old = mdotf.oldTime();

scalar URFU = mesh.equationRelaxationFactor("U");

// Rhie-Chow interplation
mdotf =   (U_avg_f & mesh.Sf()) - ( (DUf*( gradp_f - gradp_avg_f)) & mesh.Sf() )
        +(scalar(1) - URFU)*(mdotf.prevIter() - (U_avg_prevIter_f & mesh.Sf()))
        + DTf*( mdotf_old - (U_old_f& mesh.Sf()));

//- resetting pressure correction
pp.internalField() = scalar(0.0);
pp.correctBoundaryConditions();

//- Assembling
fvScalarMatrix ppEqn
(
    - fvm::laplacian(DUf, pp, "laplacian(pDiff,pp)")
    + fvc::div(mdotf)
);

ppEqn.setReference(pRefCell, pRefValue);

//- Solution
ppEqn.solve();

//- update of the fluxes based on the flux correction
mdotf += ppEqn.flux();

//- update of the variable pressure and velocity based on the pressure correction solution and pressure correction relaxation factor
scalar URFPP = mesh.equationRelaxationFactor("pp");
p += URFPP*pp;
p.correctBoundaryConditions();

U -= fvc::grad(pp)*DU;
U.correctBoundaryConditions();


