//SIMPLE correction
volScalarField DU = 1.0/( UEqn.A() );

const surfaceVectorField ed = mesh.delta()()/mag(mesh.delta()());

rho = thermo->rho();

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_prevIter_f = linearInterpolate(U.prevIter());
surfaceVectorField U_avg_old_f = linearInterpolate(U.oldTime());
surfaceVectorField U_avg_f = linearInterpolate(U);

surfaceScalarField rhofd = upwind<scalar>(mesh,mdotf).interpolate(rho);
surfaceScalarField rhof("srho",fvc::interpolate(rho));
surfaceScalarField DUf("srUA",fvc::interpolate(DU,"interpolate((1|A(U)))"));

volScalarField dt 
(
  	IOobject
  	(
    	 "dt",
    	 runTime.timeName(),
    	 mesh,
    	 IOobject::NO_READ,
    	 IOobject::NO_WRITE
  	),
	mesh ,
	dimensionedScalar("dt" ,runTime.deltaT().dimensions(), runTime.deltaT().value()),
  	zeroGradientFvPatchScalarField::typeName
); 
surfaceScalarField dt_f = linearInterpolate(dt);
surfaceScalarField drhodp = linearInterpolate(thermo->psi());

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

// Rhie-Chow interplation
mdotf = rhof*(
              (U_avg_f & mesh.Sf()) - ( (DUf*( gradp_f - gradp_avg_f)) & mesh.Sf() ) 
            );
           + (scalar(1) - UURF)*(mdotf.prevIter() - ( (rhof*U_avg_prevIter_f) & mesh.Sf() ) )
           +    rhof*(DUf/dt_f)*(mdotf.oldTime() - ( (rhof*U_avg_old_f) & mesh.Sf() ) );

surfaceScalarField phid
(
  	IOobject
  	(
    	 "phid",
    	 runTime.timeName(),
    	 mesh,
    	 IOobject::NO_READ,
    	 IOobject::NO_WRITE
  	),
    mdotf*drhodp/(rhofd)
); 

surfaceScalarField pDiff = rhof*DUf;

pp = dimensionedScalar("pp", p.dimensions(), scalar(0.0));
pp.correctBoundaryConditions();

fvMatrix<scalar> ppCompEqn
(
    fvm::ddt(thermo->psi(),pp)  
  + fv::gaussConvectionScheme<scalar>
     (
         mesh,
         phid,
         tmp<surfaceInterpolationScheme<scalar> >
         (
             new upwind<scalar>(mesh,phid)
         )
     ).fvmDiv(phid, pp)
  + fvm::laplacian(-pDiff, pp)
  + fvc::div(mdotf) + fvc::ddt(rho)
);

ppCompEqn.solve();

// correct Fluxes
mdotf += ppCompEqn.flux();

p += alphaP*pp;
p.correctBoundaryConditions();

rho += alphaP*pp*thermo->psi();

boundMinMax(p, pMin, pMax);
boundMinMax(rho, rhoMin, rhoMax);

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


