Ashuthosh Shridhar ~

    // Momentum predictor
 
    MRF.correctBoundaryVelocity(U);
 
    tmp<fvVectorMatrix> tUEqn       // tUEqn -> temporary UEqn
    (
        fvm::div(phi, U)            // Convection term
      + MRF.DDt(U)
      + turbulence->divDevReff(U)   //Diffusion Term, Phi=Face-Flux, Jasak Thesis 
     ==                             //(Pg.80)
        fvOptions(U)                
    );

By equating the Momentum Equation to fvOptions instead of $\nabla P$ , we can add even more source terms.

	fvVectorMatrix& UEqn = tUEqn.ref();
 
    UEqn.relax();
 
    fvOptions.constrain(UEqn);
 
    if (simple.momentumPredictor())
    {
        solve(UEqn == -fvc::grad(p));    
 
        fvOptions.correct(U);
    }

Velocity field is calculated by taking a guess value of pressure.

References