Physics-Based Animation: PIC-FLIP Fluids

In this assignment you will learn to implement a standard hybrid particle-grid fluid simulator that solves on the GPU.

Rendered assignment output played back at high-speed

WARNING: Do not create public repos or forks of this assignment or your solution. Do not post code to your answers online or in the class discussion board. Doing so will result in a 20% deduction from your final grade.

Checking out the code and setting up the python environment

These instructions use Conda for virtual environment. If you do not have it installed, follow the instructions at the preceeding link for your operating system

Checkout the code git clone [email protected]:panuelosj/pba-assignment-fluids.git {ROOT_DIR}, where {ROOT_DIR}* is a directory you specify for the source code.

Next create a virtual environment and install relevant python dependencies install.

cd {ROOT_DIR}
conda create -n csc417  python=3.12 -c conda-forge
conda activate csc417
pip install -e . 

Optionally, if you have an NVIDIA GPU you might need to install CUDA if you want to use the GPU settings

conda install cuda -c nvidia/label/cuda-12.1.0

Assignment code templates are stored in the {ROOT_DIR}/assignment directory.

WINDOWS NOTE: If you want to run the assignments using your GPU you may have to force install torch with CUDA support using

pip install torch torchvision torchaudio --index-url https://download.pytorch.org/whl/cu124

Installation without conda

If you are having too many problems with Conda or prefer to use another package manager, we recommend using UV. If you do not have it installed, follow the instructions at the preceeding link for your operating system

Next, create a virtual environment and install relevant python dependencies:

cd {ROOT_DIR}
uv venv
uv pip install -e . 

If you opted to use UV, you can run the examples using:

uv run python main.py <arguments-for-tests>

Tools You Will Use

1. NVIDIA Warp -- python library for kernel programming
2. PyTorch -- python library for array management, deep learning etc ...

Running the Assignment Code

cd {ROOT_DIR}
python main.py --scene=tests/{SCENE_JSON_FILE}.json

By default the assignment code runs on the cpu, you can run it using your GPU via:

python main.py --scene=tests/{SCENE_JSON_FILE}.json --device=cuda

Finally, the code runs, headless and can write results to a USD file which can be viewed in Blender:

python main.py --scene=tests/{SCENE_JSON_FILE}.json --usd_output={FULL_PATH_AND_NAME}.usd --num_steps={Number of steps to run}

Occasionaly on windows the assignment will fail to run the first time, but subsequent attempts should work fine

Assignment Structure and Instructions

1. You are responsible for implementing all functions found in the assignments subdirectory.
2. The tests subdirectory contains the scenes, specified as python files, we will validate your code against.
3. This Google Drive link contains output from the solution code that you can use to validate your code. The output consists of USD (Universal Scene Description) files which contain simulated results. These can be played back in any USD viewer. I use Blender. You can output your own simulations as USD files, load both files in blender and examine the simulations side-by-side. A blender file is available at the blender folder if you want to render out with the same material as the gif above.

Background and Resources

PIC-FLIP was introduced into graphics in this paper by Zhu and Bridson. The SIGGRAPH Coursenotes is also a good resource for implementing this technique, as is a recent video by Sebastian Lague.

PIC-FLIP Pipeline

In this assignment, you will be implementing a fluid simulator according to the following PIC-FLIP Pipeline:

This method involves doing operations in two data structures: particles and a grid. As you learned in class, and as presented in the pipeline above, certain operations are more amenable to either one of the data structures. You will be responsible for implementing the operations for advection, body force (gravity), pressure projection, and the grid-to-particle and particle-to-grid interpolations that transfers data between the particle and the grid data structures.

Advection

Advection moves particles through the velocity field. For each particle, update its position using the simple forward Euler timestepping:

Advection is simply physically moving the particle forward by shifting its position based on its current velocity. The positions of these particles are markers that dictate which portion of the domain is considered as fluid (any grid cell that contains a particle is considered as being fluid, whose velocity is then solved for in subsequent steps).

Body Force (Gravity)

Body forces are external forces that act uniformly throughout the fluid volume. In this assignment, we primarily consider gravity, which accelerates the fluid downwards. Body forces are applied with simple forward Euler again, this time operating on velocities and, taking g to be the gravitational acceleration vector (remember that acceleration is the time-derivative of velocity):

Note that this can be done in either the particle or grid regimes. In this assignment, the implementation uses it inside the grid (the warp tid's are given to you and should make that clear).

The Staggered Grid

As previously mentioned, we need to operate on both a particle data structure as well as a grid data structure. The grid data structure, because of the spatial derivatives involved, have to be of a specific form. The PIC-FLIP method uses a staggered grid (also called a MAC grid - Marker-and-Cell grid), which visually have the following schematic:

Here, pressure and velocity components all "live" in different places in space. More specifically:

• Pressures (green squares) are stored at cell centers (i, j, k)
• Velocity components are stored at cell faces:
  • x-component (vx, orange circles) is stored at faces between cells in the x-direction: (i-1/2, j, k)
  • y-component (vy, yellow circles) is stored at faces between cells in the y-direction: (i, j-1/2, k)
  • z-component (vz, not pictured) is stored at faces between cells in the z-direction: (i, j, k-1/2)

Because we cannot have half indexing, we simply shift our indexing by the appropriate amounts for the four different grids (see diagram for integer indexing). Notice that from the point of view of a pressure sample (i, j, k), the left velocity sample is (i, j, k) and the right velocity sample is (i+1, j, k). Likewise, from the point of view of a velocity sample (i, j, k), the left pressure sample is (i-1, j, k) and the right pressure sample is (i, j, k). This specific layout ensures that taking spatial derivatives end up placed where they need to be in space (the divergence of velocities naturally falls on cell centers where it can be added to pressures for example), and avoids nasty nullspaces that causes instabilities.

Pressure Projection

Pressure projection is the step that enforces incompressibility (∇·u = 0) by correcting the velocity field using the pressure gradient. Note that this will be the most involved part of the assignment, and involves solving the following equations:

and

Numerically, this involves three steps:

1. Computing the divergence of the current velocity field at each cell center, then scaling it by ρ / dt (RHS of the top equation).
2. Solving the Poisson equation ∇²p = (ρ/dt)∇·u to find the pressure field (notice that the RHS is exactly what's computed at step 1).
3. Correcting the velocity using the pressure gradient: u_new = u - (dt/ρ) ∇p (the bottom equation).

Note that boundary conditions need to be handled carefully: set velocities to zero at solid boundaries, and set pressures to zero at free surface (air) boundaries. state.grid_occupancy is already computed for you and will give the type of cell at each grid center (valued at 1 for fluid, 0 for air, and -1 for solid).

To perform this solve, you will need to use discrete operators for divergence, gradient, and the Laplacian.

Divergence Operator

The divergence operator computes ∇·v at each cell center from the staggered velocity field, essentially being a measure of how much fluid is entering vs exiting some cell (and thus "converging" vs "spreading out"). Note that this value "lives" at cell centers (so the same spot as pressures).

For a cell at (i, j, k), the divergence is:

where:

• u[i,j] is the x-velocity at face (i-1/2, j, k)
• u[i+1,j] is the x-velocity at face (i+1/2, j, k)
• v[i,j] is the y-velocity at face (i, j-1/2, k)
• v[i,j+1] is the y-velocity at face (i, j+1/2, k)
• Similar for z-component in 3D.

For incompressible flow, we want ∇·u=0, which is enforced by the pressure projection step. Only compute divergence in fluid cells; set it to zero for air and solid cells (since we won't need to update velocities at those cells).

Laplace Operator

The Laplace operator (∇²) is used in the Poisson equation to solve for pressure. Notice that this too must exist at cell centers (since we have ∇²p = (ρ/dt)∇·u and we know that RHS must be defined on cell centers, and the LHS must agree in the domain).

The discrete Laplacian at cell center (i, j, k) is:

This is the standard 5-point stencil in 2D, you can derive a similar 7-point stencil in 3D. The Poisson equation ∇²p = divergence is solved iteratively using either the Jacobi method or Gauss-Seidel with Successive Over-Relaxation (SOR) (see later section). Boundary conditions must be handled:

• Free surface (air cells): p = 0 (Dirichlet boundary condition)
• Solid cells: p_solid = p_fluid (Neumann boundary condition, mirror the fluid pressure)

Gradient Operator

The gradient operator computes ∇p for the pressure projection step. This now must live on cell faces since you need to compute central finite differences of pressures which live on cell centers.

On a staggered grid, the gradient is computed at face locations:

• At x-face (i-1/2, j, k): ∇p_x = (p[i] - p[i-1]) / dx
• At y-face (i, j-1/2, k): ∇p_y = (p[j] - p[j-1]) / dy
• At z-face (i, j, k-1/2): ∇p_z = (p[k] - p[k-1]) / dz

The velocity is then updated using: v_new = v_old - ∇p × dt / ρ. This correction subtracts the divergent component of the velocity, making it incompressible. Handle boundaries carefully: if either adjacent cell is solid, set the velocity component at that face to the solid velocity (zero for purposes of this assignment).

Jacobi Iteration

We want to find a solution for pressures that solves the following equation:

Which if we discretize according to the above Laplacian operator, looks like:

where P is the pressure at our current grid cell, sum(P_nbr) is the sum of pressures at the neighbouring grid cells (i±1, j±1, k±1), and n_nbr is the number of these neighbours (4 in 2D, 6 in 3D). Assuming we have an accurate estimate for the neighbouring pressures, this means we can isolate for the pressure at our current point:

Taking this to be true everywhere, this means that we can start off with some guess for the pressures everywhere, and iteratively compute new pressures. This gives increasingly better estimates for the solution pressure that satisfies this equation, which hopefully converges to the true solution. This iterative method of guessing and recomputing is called Jacobi iteration, which updates all grid cells in parallel using values from the previous iteration.

Key Implementation Details:

1. Read from old values: All cells read from p_old simultaneously, making this method highly parallelizable
2. Boundary conditions:
   • Air cells (free surface): Set p = 0 (Dirichlet boundary condition)
   • Solid cells: Use mirror condition (p_neighbor = p_center) for Neumann boundary condition
   • Fluid neighbors: Use their actual pressure values
3. Relaxation: Apply weighted averaging: p_new = α × p_jacobi + β × p_old, where α + β = 1. Standard Jacobi uses α = 1, β = 0, but relaxation can improve convergence. This is exposed as a parameter in the scene files.
4. Iteration: Repeat this process many times (might take ~500-1000 iterations to have the fluid not lose volume) until convergence

The Jacobi method is simple and parallelizable but converges slowly.

Particle-to-Grid

Particle-to-grid (P2G) transfers particle mass and momentum to the grid. For each particle:

1. Find the base grid cell containing the particle (convert particle position to grid coordinates)
2. Compute interpolation weights using cubic B-spline kernels (or simpler linear weights) for a 3×3×3 neighborhood (see diagram below)
3. Transfer mass and velocity to each grid node in the neighborhood:
   • Mass: m_grid += weight × m_particle
   • Velocity: u_grid += weight × v_particle

Grid-to-Particle

Grid-to-particle (G2P) transfers the updated grid velocity back to particles. This is where PIC-FLIP blending happens:

1. Interpolate the new grid velocity at the particle position using trilinear interpolation (PIC component) ( u_PIC = interp(u_new_grid) )
2. Interpolate the change in grid velocity at the particle position using trilinear interpolation (FLIP component) ( Δv_FLIP = interp(u_new_grid - u_old_grid) )
3. Apply FLIP blending: u_particle = α × u_PIC + (1-α) × (v_particle_old + Δv_FLIP)

where:

• PIC (α = 1.0): Directly use the interpolated new grid velocity
• FLIP (α = 0.0): Add the velocity change to the old particle velocity
• PIC-FLIP (0 < α < 1): Blend both approaches for better stability and detail preservation

Typical values for α are around 0.01-0.1, giving mostly FLIP with a small PIC component for stability.

For interpolation, use trilinear interpolation weights based on the fractional position of the particle within its containing grid cell:

Debugging Hints

1. I suggest starting by implementing advect and apply_gravity, then g2p and p2g. These files can be debugged independently of pressure projection (it should just let your particles fall down). Implementing pressure projection is near impossible without getting these files working.
2. I've included a visual in polyscope for previewing the fluid and solid cells. You can make use of these to visually reason on where grid points are. A lot of the discrete operators (gradient, divergence, laplacian) involve local stencils that simply need to query the nearby grid cells. It is very easy to have an off-by-one error here. Draw yourself a diagram on what each grid stencil should look like.
3. If you are getting memory access issues, it's likely because you're trying to query a grid point that's out of bounds. Make sure to include bounds checks that your indexing doesn't go below 0 or go above the grid limits.
4. You can test all files minus jacobi_pressure_iteration.py by switching the "use_gauss_seidel" flag on the json scene files to true. Instead of using your implementation of jacobi_pressure_iteration, it replaces it with a Gauss-Seidel solve. This method converges faster and will be helpful for checking if the rest of your simulation is working correctly.

Admissable Code and Libraries

You are allowed to use any functions in the warp and warp.sparse packages. You ARE NOT allowed to use code from other warp packages like warp.fem. You cannot use code from any other external simulation library.

Hand-In

We will collect and grade the assignment using MarkUs

Late Penalty

The late penalty is the same as for the course, specified on the main github page.