fvScalarMatrix mfEqn
(
        fvm::ddt(rho,mf) 
	+ mvConvection->fvmDiv(phi, mf)
        - fvm::laplacian(1.47*turbulence->mut(), mf) //let 1/Sc = 1.47
	== parcels.Srho(mf)
);

mfEqn.solve();