[randolph]

In the buoyantBoussinesqPimpleFoam, the gradient of pressure and the density term is calculated as the following. In my understanding, the surface normal gradient is first calculated and then used to reconstruct the gradient.

   solve
         (
             UEqn
              ==
             fvc::reconstruct
             (
                 (
                   - ghf*fvc::snGrad(rhok)
                   - fvc::snGrad(p_rgh)
                 )*mesh.magSf()
             )
         );

My question is why not just simply calculating the gradient as:

-fvc::grad(p_rgh)

Are these two approaches equivalent? If not, why is the first approach better than the second?

Thanks,
Rdf