Position Based Fluids

13 Dec 2022

Notes Computer Animation Position Based Fluids

§ Description

For my final project in CSCE-451, I reimplemented Position Based Fluids (PBF) by Macklin and Müller (2013)^[1].

Essentially, PBF is just an extension of smoothed-particle hydrodynamics (SPH) into the Position Based Dynamics framework (Müller 2007).

In contrast to SPH, the major benefit that PBF gives is stability, allowing us to take larger physics time-steps. Typically SPH requires many physics substeps for a single animation frame; whereas PBF can remain stable with a single physics step per frame.

In practice, PBF with $\sim 4$ physics substeps per animation frame gives reasonable results for real-time applications.

§ Algorithm

The PBF algorithm is very simple, and is easily parallelisable.

§ Euler Integration

Integration

for all particles $i$ $i$ do
- $\mathbf{v}_i$ $\leftarrow$ $\mathbf{v}_i + \Delta t\ \mathbf{F}(\mathbf{x}_i)$ # apply forces
- $\mathbf{x}^{\star}_i$ $\leftarrow$ $\mathbf{x}_i + \Delta t \mathbf{v}_i$ # predict position

First, the algorithm does an unconstrained step by applying external forces and particle velocities; but the particles do not move quite like a liquid.

§ Particle Sorting

A lot of our computations are dependent on a particle's neighbours, e.g.

Neighbour Computations

for all particles $i$ $i$ do
- neighbours $\leftarrow$ $N_i(\mathbf{x}^{\star}_i)$ # find neighbouring particles
- for all particles $j$ $j$ in neighbours do
  - ComputeWeightedSum( $i$ , $j$ )

Naively, for each particle we look at every other particle, this ends up being $\mathcal{O}\left(n^2\right).$

A better way is to discretise the simulation space into a uniform grid of bins, insert the particles to the bins, then sort.

Using counting sort, this step becomes $\mathcal{O}(nk)$ ^[2]. And with a prefix-sum implementation^[3], we can sort in parallel.

§ Jacobi Optimizer

Optimizer

while iter < solverIterations do
- for all particles $i$ $i$ do
  - $\rho_i$ $\leftarrow$ ComputeDensity( $i$ )
- for all particles $i$ $i$ do
  - $\lambda_i$ $\leftarrow$ LagrangeMultiplier( $i$ )
- for all particles $i$ $i$ do
  - Compute $\Delta\mathbf{p}_i$
- for all particles $i$ $i$ do
  - $\mathbf{x}^\star_i$ $\leftarrow$ $\mathbf{x}^\star$

This step projects the particles' positions back into their expected fluid position using a constraint function.

Constraint Function

Next, define the constraint function by

C_i(\mathbf{p}_i) = \frac{\rho_i}{\rho_0} - 1.

Due to the integration step we have that $C(\mathbf{p}_1, \dots, \mathbf{p}_n) \neq 0.$ As a result, we optimize for a positional correction $\Delta\mathbf{p} \approx \nabla C_i(\mathbf{p})\lambda$ such that $C(\mathbf{p} + \Delta\mathbf{p}) = 0,$ where $\lambda$ is a Lagrange multiplier given by

\lambda_i = -\frac{C_i(\mathbf{p})}{\sum_k\|\nabla_{\mathbf{p}_k}C_i\|^2 + \varepsilon}.

$\Delta\mathbf{p}$

The $\Delta\mathbf{p}$ computation is straightforward,

\Delta\mathbf{p}_i = \frac{1}{\rho_0}\sum_j\left(\lambda_i + \lambda_j + s\right)\nabla W(\mathbf{p}_i - \mathbf{p}_j, h)

where $s$ is the tensile instability, and $\nabla W(\cdot)$ is the spiky kernel (Müller 2003)[^3]; the gradient of the poly6 kernel .

The tensile instability adds an artificial pressure that encourages particles to clump together to mimic surface tension.

Typically, the optimizer is not run until convergence -- rather, it is run for a few (4-6) steps which gives results that look good visually.

§ Viscosity and Vorticity

Viscosity and Vorticity

for all particles $i$ $i$ do
- $\mathbf{v}_i$ $\leftarrow$ $\frac{1}{\Delta t}\left(\mathbf{x}^{\star}_i - \mathbf{x}_i\right)$ # update velocities
- $\mathbf{v}_i$ $\leftarrow$ $\mathbf{v}_i + \frac{c}{\rho_i} \nu_i$ # apply XSPH viscosity
- $\mathbf{v}_i$ $\leftarrow$ $\mathbf{v}_i + \varepsilon\Delta t\,\left(\hat{\eta}_i \times \omega_i\right)$ # apply vorticity confinement
- $\mathbf{x}_i$ $\leftarrow$ $\mathbf{x}^{\star}_i$ # update position

To make the results look more visually appealing, XSPH viscosity and vorticity confinement is applied.

Vorticity confinement reintroduces energy that was lost in the integration and optimization steps making the simulation look nicer, however, this technique is not physically based.

§ Implementation

I used CUDA to implement the entire algorithm on the GPU to exploit its highly parallel nature.

Originally, I had a lot of issues with performance in simulating ~150K particles in real time.

I had a lot of CPU $\longleftrightarrow$ $⟷$ GPU data transfers which introduced some latency.
- Instead, I used the CUDA/OpenGL interop features to directly map the OpenGL VBOs to the PBF simulation buffer.
GPUs are optimised for single precision.
- A lot of the library functions I was using were causing type promotion from floats to doubles, preventing this gave me a 10x speed up.

My renderer simply reuses the code from the particle simulation labs, but I replaced the shaders with one that did Blinn-Phong shading on point-sprite spheres based on Green (2010)^[4].

Initially, I wanted to implement screen space fluid rendering^[4:1] with an anisotropic kernel^[5]. Unfortunately, I only managed to implement the kernel and not the shaders.

Macklin and Müller. (SIGGRAPH 2013). Position Based Fluids ↩︎
R. Hoetzlin. (GDC 2014). Fast Fixed-Radius Nearest Neighbors: Interactive Million-Particle Fluids ↩︎
(Fall 2022). UPenn CIS565 - GPU Programming, Project 2: Stream Compaction ↩︎
S. Green. (GDC 2010). Screen Space Fluid Rendering for Games ↩︎ ↩︎
Yu and Turk. (SIGGRAPH 2010). Reconstructing Surfaces of Particle-Based Fluids Using Anisotropic Kernels ↩︎