What are projection methods

Question

Quoting from Solenthaler et. al. Predictive-Corrective Incompressible SPH (ACM Transactions on Graphics, Vol. 28, No. 3, Article 40, Publication date: August 2009) (PDF link here)

These incompressible SPH (ISPH) methods first integrate the velocity field in time without enforcing incompressibility. Then, either the intermediate velocity field, the resulting variation in particle density, or both are projected onto a divergence-free space to satisfy incompressibility through a pressure Poisson equation.

What does the author mean by "project?? Is there a simpler way of understanding this operation? I tried reading other articles, but they go even more deep saying solenoidal vectors etc. I am looking for a simple explanation at first to understand the concept.

Answer 1

Projection methods typically split up the solution of transient Stokes or Navier-Stokes methods into the solution of two separate problems - one for the velocity, one for the pressure. For Stokes, the simplest version (by Chorin/Temam) goes as follows for a timestep $k$:

1. Compute an intermediate velocity $u^*$ by solving $$\frac{(u^*-u^k)}{dt} - \nu\Delta u^* = f(t^k)$$
2. Correct the velocity by solving $$\begin{cases} \frac{(u^{k+1}-u^*)}{dt} + \nabla p^k &= 0\\ \nabla \cdot u^{k+1} &= 0\\ u^{k+1}\cdot n &= 0 \quad \text{on the boundary} \end{cases}$$

This above step is the reason these splitting schemes are referred to as projection methods - the weak form of (2) can be interpreted as the scaled $L^2$ projection of $u^*$ onto a divergence-free $u^{k+1}$ (after multiplying by a test function, $p$ can then be viewed as a Lagrange multiplier to enforce the divergence free constraint).

A note: since pressure is undetermined in the second step, we need to determine it. Based on the assumption that $\nabla \cdot u^{k+1}$ must be $0$, we can take the divergence of the second equation to get $$\nabla \cdot \frac{u^*}{dt} + \nabla\cdot \nabla p^k = 0,$$ which yields a Poisson equation for $p$ given $u^*$. The overall scheme is then

1. Solve for $u^*$ through $$\frac{(u^*-u^k)}{dt} - \nu\Delta u^* = f(t^k)$$
2. Solve for $p$ with $$\Delta p = \frac{1}{dt}\nabla \cdot u^*$$
3. Solve for $u^{k+1}$ through $$\frac{(u^{k+1}-u^*)}{dt} + \nabla p^k = 0.$$

See this paper for a great overview of available projection methods.

Answer 2

To give an easier example, consider the ODE for rotation

\begin{align} \dot x=-ay\\ \dot y=ax \end{align}

If one solves that with the Euler method, the next point is found in the direction of the tangent of the circle, increasing the radius. \begin{align} x(t+h)^2+y(t+h)^2&=(x(t)-ahy(t))^2+(y(t)+ahx(t))^2\\ &=(1+a^2h^2)(x(t)^2+y(t)^2) \end{align}

However, since one knows that the exact solution of this ODE preserves the initial radius

$$x(t)^2+y(t)^2=r^2=x(0)^2+y(0)^2$$

one can correct this drifting away by projecting the numerical solution back to that circle.

$$(\tilde x(t+h),\tilde y(t+h)) =\frac{r}{\sqrt{x(t+h)^2+y(t+h)^2}}(x(t+h),y(t+h))$$

The angular velocity will still be different, but the approximation will still be better.

In general, if you have a system of first integrals for the problem, you may want to return, after a certain number of steps, to the set satisfying the first integrals. Using Newton's method with the Moore-Penrose pseudo-inverse will usually find a rather short path to that set.

Finding the closest point on some set usually is called "projection". This is conform to the characteristic of orthogonal projections in linear algebra.

Answer 3

If you're confused why this operation is called "projection", note that the equation $$ \text{div} \, \mathbf u = 0 $$ means that we only consider velocity fields that are divergence free. This is no different than asking to find a solution of some problem subject to a constraint of the form $Bx=0$ where $B$ is a matrix with fewer rows than columns. The set of all vectors $x$ that satisfy this form a subspace (think of a plane in 3-space, or a hyperplane in $n$-space). What projection algorithms do is to find some approximation for a simpler subproblem, and then "project" back to the hyperplane of functions that are divergence free, where the projection is really just to be considered what one does in regular, finite dimensional algebra or geometry.