Assignment 3: Smoke Simulation

**Posted:** Thursday, 13 September 2018

**Updated:** Sunday, 16 September 2018

**Due:** Thursday, 27 September 2018

In this assignment you will implement a basic 2D incompressible fluid simulator, and extend it with various improvements to advection, pressure projection, and boundary handling.

Implement a 2D fluid solver using a staggered grid with semi-Lagrangian advection, a Gauss-Seidel pressure solver, and buoyancy forces for smoke. As you can see, there are several different components that need to work correctly even for this base requirement, so make sure you get started on this assignment early!

To help you get started, I’ve provided a `Grid<T>`

class that represents a regular \(m\times n\) grid with spacing \(\it dx\), the origin \(\mathbf x_0\) specifying the position of the \((0,0)\)th grid point, and values of type `T`

at each point. Note that if you want to divide the domain, say \([0,1]\times[0,1]\), into \(n\times n\) grid cells, the pressure grid points should be located at grid cell *centers*, \(((i+\frac12)/n,(j+\frac12)/n)\); choose the grid origin accordingly.

The **first** thing you should do is implement a staggered grid to store velocity values. One way to do this is to simply have the staggered grid contain two regular grids, each with a slightly different size and an appropriately offset origin. Implement a function to compute the interpolated velocity at any given point \(\mathbf x\) by performing interpolation on the underlying grids.

At this point, it is useful to be able to visualize the data in a given scalar or vector field. For scalar fields, I suggest creating a function that draws a \({\it dx}\times{\it dx}\) square centered at the location of each grid point \((i,j)\), coloured according to the corresponding value \(c_{ij}\). Some easy colour maps you could use are \(c \mapsto (c,c,c)\), which goes from black to white over \([0,1]\), and \(c \mapsto (1+c, 1-c^2, 1-c)\), which goes from blue to white to red over \([-1,1]\). For vector fields, drawing the vector components on the grid faces is quite unilluminating. Instead, draw an arrow at the center of each grid cell by averaging the velocity components on the adjacent faces. If drawing full 3D arrows becomes slow for large grids, just draw a line segment.

The **second** step is to implement semi-Lagrangian advection. I suggest doing this in two steps: first advecting particles forward through the velocity field, then advecting scalar and vector fields using backtracking. (i) Advecting a particle \(\mathbf x\) through a velocity field \(\mathbf u\) essentially involves integrating the ODE \(\dot{\mathbf x} = \mathbf u(\mathbf x)\). Write a method to do this over a single time step \(\Delta t\), using a simple second-order scheme like explicit midpoint or Heun’s method. (ii) To do semi-Lagrangian advection on a scalar field, you will require two grids: one which gives the initial state, and one into which you will write the advected output. For each point in the output grid, use the particle advection routine for time \(-\Delta t\) to find the position to interpolate from. To do semi-Lagrangian advection on a vector field, simply perform scalar field advection on each component.

Test your code on the following velocity field, which is a simple divergence-free vector field on \([0,1]\times[0,1]\): \[\begin{align} \mathbf u(x,y) = \begin{bmatrix}\sin(\pi x)\cos(\pi y) \\ -\cos(\pi x)\sin(\pi y)\end{bmatrix}. \end{align}\] Create a staggered grid on \([0,1]\times[0,1]\) and initialize it with this velocity field. Also create a collection of particles distributed in \([0,1]\times[0,1]\). At each frame, advect the particles through the velocity field, so that they circulate smoothly around the square. Finally, create a scalar field with some interesting initial values (e.g. \(c=0\) if \(x<\frac12 \mathrel{\rm xor} y<\frac12\), \(c=1\) otherwise), and perform semi-Lagrangian advection on it as well. Verify that the motion of the values on the scalar field is consistent with the motion of the particles (apart from the numerical diffusion, which will be a lot).

The **third** step is to add pressure projection. As a warm-up, you may want to first implement a function to compute the gradient of a scalar field and write it into a vector field, and another to compute the divergence of a vector field and write it into a scalar field. You will need to be careful about which entries in the staggered grid are adjacent to the chosen cell in the regular grid, and vice versa. At grid boundaries, set the normal component of the gradient to \(0\).

Perform pressure projection by computing the divergence of the input vector field, computing a scalar field \(p\) whose Laplacian is the said divergence, and subtracting its gradient from the velocity field. Compute \(p\) by implementing a Gauss-Seidel solver that starts with a given initial guess (e.g. all zeroes) and updating one cell at a time. Be careful to use the appropriate stencil when next to grid boundaries. For the update order, use a red-black scheme as mentioned in class: first update all the grid cells \((i,j)\) for which \(i+j\) is even, then all the cells for which it is odd. If you know OpenMP, you can parallelize this step very easily. Run the Gauss-Seidel solver until the maximum absolute divergence on any cell is sufficiently small, or a maximum iteration count is reached.

Test your code on some simple input velocities, for example (i) a constant nonzero velocity everywhere (except on boundary faces, of course), and (ii) a velocity that is zero everywhere except in a region in the interior, like the example seen in class. Verify that the divergence of the output velocity field is nearly zero. Try this on small grids first, since you will need a lot of iterations if the grid is large. It may be helpful to visualize the computed pressure field while debugging.

**Finally**, add smoke as a scalar field advected using the current velocity, and apply a buoyancy force proportional to the amount of smoke. A standard example is to define a smoke source near the bottom of the domain, which is a region in which you reset the smoke amount to \(1\) on each time step. For buoyancy, since smoke is defined on cell centers, you will have to do some interpolation to apply forces on cell faces. Then, you should be able to simulate a plume of smoke rising over time and creating an interesting flow.

Note that the entire domain may fill up with smoke if you run the simulation for a long time. To avoid this, you could add a decay term to the evolution of the smoke field to prevent this (which will amount to just decreasing the smoke density in each grid cell a little on each time step).

In your submission, include the following results:

a video of a collection of particles and a scalar field being advected by the vector field \(\mathbf u\) defined in the advection section;

images of an example non-divergence-free vector field, the corresponding pressure field computed in the projection step, and the final divergence-free vector field; and

a video of a smoke plume rising from a source.

Completing the base requirement is worth 60% of the points for this assignment. Each of the following four extensions is worth an additional 20%, so you need to do at least two of them to have the possibility of getting 100% points. Note that the total points for the assignment are capped at 100%. The extensions are mostly independent and do not have to be done in sequential order, i.e. you can choose any two or more to implement in your submission, not necessarily the first two. Your implemented extensions also do not have to be mutually compatible.

Many of the components below are also described in Bridson and Müller-Fischer fluids notes, if you need more detail.

There are a couple of simple improvements you can make to your smoke simulator to make it more effective and useful.

First, replace the linear interpolation scheme used in the semi-Lagrangian advection lookup with a monotone cubic interpolation, as described in Appendix B of Fedkiw et al.’s “Visual Simulation of Smoke” (2001). Note that their modification is a *necessary* condition but not a *sufficient* one, so you may see a very small amount of overshooting, though much less than with traditional cubic interpolation. If you want to fix this, a sufficient (but no longer necessary) condition is to force the slopes to have the same sign as \(\Delta_k\) and with absolute value no greater than \(3|\Delta_k|\) (in their notation).

Run the scalar field advection example from the base requirement again, this time with monotone cubic interpolation. You should observe much less numerical diffusion than before, while still having negligible overshooting (i.e. the values of the scalar field remain within the range of the original data).

Second, modify your simulator to support solid boundaries in the interior, not just at the ends of the grid. This requires attaching to each grid cell a flag that labels it as solid or fluid, and then modifying both the advection step (to not sample values inside solids) and the pressure step (to set appropriate boundary conditions at fluid/solid interfaces). Refer to Fedkiw et al. (2001) Sec. 4 and/or the lecture slides for more details.

Create an animation of the rising smoke plume, but with a solid obstacle (e.g. a voxelized disk or a rectangle) in its path, so that the plume is deflected realistically by the presence of the obstacle. Draw the obstacle cells with a different colour so that the obstacle itself is visible.

As you will have seen by now, semi-Lagrangian advection is significantly diffusive, and this problem is reduced but far from eliminated by higher-order interpolation. To reduce it further, implement a FLIP advection scheme, as described in Zhu and Bridson’s “Animating Sand as a Fluid” (2005), Sec. 4. You will have to choose how many FLIP particles to initialize at the beginning of the simulation: there shouldn’t be so few particles that some grid cells are left uncovered, nor so many that it becomes a computational bottleneck. A good choice is to create 4 particles in a randomly jittered \(2\times2\) arrangement per grid cell (8 particles in \(2\times2\times2\) in 3D). Also, FLIP by itself is unstable, so you will need to blend it with a small amount of PIC to smooth out the oscillations; typical values are between 1% to 5%.

For testing, it is a good idea to first advect a scalar field, e.g. smoke density, using particles. This way you will be able to check that the particle-to-grid transfer routine, which is new in this component, is behaving as expected. You may want to visualize the values on particles with different colours, so that you can compare them to the resulting values on the grid.

Demonstrate your implementation by simulationg the smoke plume with FLIP advection, which should result in more detailed vortices in the flow. In your animation, visualize the smoke values on the grid as well as the positions of the particles, so that they can be compared.

To permit smooth flow around sloped boundaries, implement the variational pressure solver described in Batty et al.’s “A Fast Variational Framework for Accurate Solid-Fluid Coupling” (2007), Sec. 2.1. This should not actually require a great deal of changes to the pressure projection: only the coefficients of the divergence and Laplacian will change from \(-1\)’s and \(4\)’s to some intermediate values depending on the area covered by the obstacle(s). The area fractions can be quite approximate, as long as they are smooth: you could use the length fraction along the line segment joining the cell centers, as suggested by Batty et al., or simply something based on the signed distance function of the obstacle. In any case, you should precompute them since they are constant over time for static obstacles.

You may have to modify the advection step to account for the more general obstacle shapes, maybe by taking into account the area weights in the interpolation, by projecting the lookup position to the outside of the obstacle, or by extrapolating velocities into the obstacle. If you tried something like this, describe briefly what you did in your report.

Show an animation of the smoke plume with a non-grid-aligned obstacle, such as a disk or a rotated rectangle. The smoke should flow smoothly around the obstacle without being affected by stairstepping artifacts.

Extending the smoke simulator to handle liquids primarily requires the addition of a surface tracking algorithm, and modification of the pressure solver to handle Dirichlet boundary conditions at the free surface.

For surface tracking, the simplest approach is to use marker particles, as described in Foster and Metaxas, “Realistic Animation of Liquids” (1996), Sec. 4.1. If you’ve implemented the FLIP component, the FLIP particles will serve this purpose. Otherwise, initialize a collection of particles in the liquid region at the start of the simulation (see the FLIP section for advice on the particle distribution), and advect them with the velocity field at every time step. Label non-solid cells as fluid or empty depending on whether they contain any marker particles. You may also need to do some velocity extrapolation so that you can sample velocities outside the current liquid region, by marching outward from the current fluid cells. Section 6.3 of the Bridson and Müller-Fischer notes describes the principled way to do this, but something simpler like a Dijkstra-type algorithm will be sufficient for the purposes of this assignment.

For the free surface boundary conditions in the pressure solver, fix the pressure to be zero at all air cells, as described in class.

Two standard examples to demonstrate liquid simulation algorithms are (i) a breaking dam, i.e. a tall, initially static rectangle of liquid that quickly flows under gravity (e.g. Zhu and Bridson (2005), Fig. 5), and (ii) a drop of liquid falling into a pool of static liquid, creating a splash (e.g. Foster and Metaxas (1996), Fig. 4(e)-(h)). Create an animation of one or both of these scenarios, visualizing the grid cell flags (liquid/obstacle/empty) in different colours, as well as the positions of the marker particles.

A submission form will be opened on Moodle for this assignment. There, you should submit a zip file containing the following items:

a report in PDF format describing which extensions you implemented in this assignment, any images required in the base and extensions, and a brief description of your submitted video(s),

a video or set of videos showing the results of your implementation, and

a copy of your implementation code.

Any additional information included in the report about how your implementation works will also be useful for my grading, and will help me give you more detailed feedback on your assignment.