// Relative momentum predictor

    volVectorField divPhiUrel
    (
  	IOobject
  	(
    	 "divPhiUrel",
    	 runTime.timeName(),
    	 mesh,
    	 IOobject::NO_READ,
    	 IOobject::NO_WRITE
  	),
	fvc::div(phi, Urel)
    );
    divPhiUrel.write();

    volVectorField gradP
    (
  	IOobject
  	(
    	 "gradP",
    	 runTime.timeName(),
    	 mesh,
    	 IOobject::NO_READ,
    	 IOobject::NO_WRITE
  	),
	fvc::grad(p)
    );
    gradP.write();

    tmp<fvVectorMatrix> UrelEqn
    (
        fvm::div(phi, Urel)
      + turbulence->divDevReff(Urel)
      + SRF->Su()
     ==
        sources(Urel)
    );

    UrelEqn().relax();

    sources.constrain(UrelEqn());

    solve(UrelEqn() == -fvc::grad(p));