The context of my question is two-phase incompressible solvers such as interFoam in OpenFOAM, but I have seen this trick used elsewhere. To handle gravity these solvers subtract the hydrostatic pressure forming a modified pressure:

$p^* = p - \rho \boldsymbol{x}\cdot\boldsymbol{g}$

and then apply Chorin's projection to Navier-Stokes equation, deriving a Poisson equation. The effects of source terms will appear on the right hand side of the latter.

What is a formal justification for this trick and does it work always or under some circumstances only? The problems should be mathematically identical, but is clearly different numerically. Comparing the resulting Poisson equations a modified pressure will contain a divergence of $\boldsymbol{g} \cdot \boldsymbol{x} \nabla\rho$ and the other a divergence of $\rho\boldsymbol{g}$ with the remaining terms staying the same.

An explanation that I always bought into is that with a modified pressure and with clear separation between incompressible phases a zero solution for pressure is now valid in parts of the domain. Therefore, a zero vector is a reasonable initial guess. With the full pressure, we are asking the linear solvers to find the hydrostatic pressure so it has to work harder. Similarly, BCs can be also set to zero. However, a second derivative of $\rho$ can be quite large around the interface so I think in some circumstances modified pressure may be detrimental and possibly leading to instabilities?

  • $\begingroup$ Have you at least convinced yourself of the usefulness of the trick for the case $\rho=const$? Or is that also in question? $\endgroup$ May 8, 2023 at 22:11
  • $\begingroup$ In $\rho = \textrm{const}$, the trick eliminates effect of gravity completely from the Poisson's solution and the hydrostatic pressure can be easily reconstructed. If there's anything else, I would be interested to hear about it. $\endgroup$ May 9, 2023 at 8:48
  • $\begingroup$ I am quite certain, that this is not the reason, hence I post this as a comment. (Representable) Floating point numbers are spaced more closely around 0, thus doing away with hydrostatic pressure shifts your working range to around 0.0 instead of 1.0e+5. Thus, one might expect less numerical errors stemming from the finite accuracy of floating point numbers. But, maybe this could also be a non-issue, and computers are accurate enough. $\endgroup$
    – Dohn Joe
    May 15, 2023 at 16:32
  • $\begingroup$ Thank you, @DohnJoe - that's already very helpful. $\endgroup$ May 16, 2023 at 19:21

1 Answer 1


It has been a long time since I was involved in the details of developing CFD codes, but here goes:

Writing general purpose CFD solvers is pretty difficult because the range of problems that the code may be applied to is VERY wide. So there are a lot of tricks that get embedded into these codes that help in particular circumstances and do not do (much) harm in other situations.

In this case we are interested in reducing the effect of round-off exacerbated by cancellation. Just to be clear, we all know there can be problems if we compute $fl(X)-fl(Y) = X-Y + \epsilon$.

When $X$ and $Y$ are large and of a similar size then the $\epsilon$ can be significant. Indeed, the relative error: $$ \frac{\epsilon}{X-Y} $$ is unbounded!

When we discretise the pressure gradient term in the momentum equations then $p_i \approx p_{i-1}$ so we could have a problem. Removing the hydrostatic head helps here because we now write: $$(p_i + \rho {\mathbf x}_i \cdot {\mathbf g })- (p_{i-1} + \rho {\mathbf x}_{i-1} \cdot {\mathbf g })= p^\star_i - p^\star_{i-1} + \rho ({\mathbf x}_i - {\mathbf x}_{i-1} ) \cdot {\mathbf g }$$

The right hand side is obviously much less subject to cancellation errors.

Notice also that this problem is much more significant if we are approximating the Poisson operator with terms like: $$ \frac{p_{i+1} - 2 p_i + p_{i-1}}{\Delta x^2} $$

Similar methods are often used when modelling turbulent flows and buoyant flows. In both cases you can derive a simple correction to the pressure which decreases cancellation errors.

We also have to consider implementation details: are we approximating the derivative by dividing by $\Delta x$; have we scaled the equations by $Re$ or some other (potentially big) non-dimensionalization? The right hand side of the Poisson equation probably involves second derivative terms so we have factors of $\Delta x^{-2}$ (at least).

We also have to worry about the condition number $\kappa(A)$ of the matrix we are "inverting." (Of course, we do not invert the matrix but the Wilkinson bound still applies.)

This number can also be large. As you say the pressure for an incompressible flow is only determined to within a constant and any sensible discretisation is going to have a similar property. In other words there is an eigen vector $p_i = 1 ~~ \forall i$ with $\lambda = 0$. People are often uncomfortable with this so they eliminate the problem by constraining the pressure artificially.

(One easy method is to fix $p_i = 0$ for some $i$.) Of course, this changes the position from $\kappa(A) = \infty$ in the singular case to $\kappa(A^\prime)$ very large (because $|| A - A^\prime|| \ll 1$. So this again shows the importance of minimising round-off.

Normally, we do not worry about round-off because discretisation errors and iterative convergence errors dominate. But very occasionally (three times in my 40+ years experience!) round-off is the answer. That is so rare that no-one ever remembers it!

Footnote: For ease of exposition I have written this response using (mostly) 1-D finite differences but the observations are valid for most forms of discretisation and multi-dimensions. Note that cancellation errors can become MORE important with higher accuracy approaches and more space dimensions.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.