I am modifying OpenFVM to deal with the neutron diffusion equation. The equation consists of a second order elliptic operator, in a way identical to the viscous diffusion term in the momentum equation. I noticed that when I perform the non-orthogonal correction, the code "blows up" if I divide the cell-center gradient by the volume of the mesh (Vp). If I don't divide by it, everything seems to be fine, meaning I get convergence and the results "make sense" compared to a case where I do not apply any non-orthogonal correction.

Any clues?