CS 567 Assignment #4: Interactive Smoke Control
Professor: Doug James

Due date: End of course (date TBD).

In this fourth and final assignment you defy the laws of physics and perform numerical magic by controlling smoke to match target keyframes of your choosing. The approach we will take is based on the SIGGRAPH 2004 paper by [Fattal and Lischinski 2004]; Graphbib page;  ACM Digital Library link (with SIGGRAPH presentation). Your implementation should run interactively for several nontrivial examples.


Starter Code (cs567.smoke): This project has significant Java starter code, primarily to support smoke simulation based on [Stam 1999; Fedkiw et al. 2001], OpenGL rendering and a skeleton code for smoke control to get you started. It is available here, with online Javadoc documentation here. In this assignment, you will modify this package as needed and however you like. Things you might want to tinker with are simulation parameters, mouse/keyboard controls, exception handling for neural input device, etc.
 Assignment Steps: Your smoke solver works out of the box. Here are the steps and issues you need to address to control smoke---most of which involve modifying the FluidSolver and SmokeControlForces objects:
  1. Driving force: Modify the skeleton code in SmokeKeyframe and SmokeControlForces.getDrivingForce(...) to implement the driving force term.  First, since the ratio of the gradient of the blurred goal density to the blurred goal density, or (GRAD rhoGoalBlur)/rhoGoalBlur, does not change over time, you can compute it once and for all in the constructor of the SmokeKeyframe--fill in this code.  Second, wire in the force in  SmokeControlForces.getDrivingForce(...). The driving force amplitude is controlled by Constants.V_f.  The drag or damping force, Constants.V_d * u, is already implemented for you in FluidSolver.velocitySolver().
  2. Gathering force: Implement the smoke gathering force in SmokeControlForces.getGatheringRate(...) using the discretization discussed in class, i.e., evaluate the effective diffusion coefficient and gradients on cell edges, then evaluate the divergence of the flux (=diffusion coeff * gradient) at the cell center using the edge-based values.
  3. Visualization: Instead of drawing the density field, try drawing other fields using a keyboard control. Drawing the smoke control forces will help you understand what your many parameters, such as blurring, are actually doing.
  4. Smoke preservation: You'll notice that due to numerical dissipation the smoke tends to decay over time, and therefore it can not match keyframes. This was one of the motivations for the particular CFD solver used in [Fattal and Lischinski 2004] but something you'll have to make the best of.  Mass preservation was also discussed in [Treuille et al. 2003]. Experiment with different normalization schemes that help match the overall keyframe density norm. Be careful not to change the normal too quickly or it may introduce artifacts, e.g., when changing keyframes.  Also simple scaling of the density field can produce other artifacts... see what you can do to make it work.
  5. Stability: By this point you should notice that the simulator is not always stable due to the explicit integration of smoke control forces, and associated time-step restriction. Also, reducing the smoke diffusion rate, Constants.SMOKE_DIFFUSION, to zero can lead to artifacts. Experiment with stability for different time steps, grid resolutions, and various parameter settings to understand how to make your simulation more robust.
  6. Preconditioned Conjugate Gradients (PCG) solver:  One of bottlenecks and numerical deficiencies of the starter code is that Gauss-Seidel relaxation is used to solve the linear systems in FluidSolver.linearSolver(...).  This solver is used for both the pressure Poisson equation solve, as well as to diffuse velocity or smoke density if there are nonzero viscosity or smoke diffusion coefficients, respectively.  Use your previous PCG implementation along with a suitable preconditioner to make the simulation faster and better. There are many possible preconditioners including Jacobi, Gauss-Seidel, incomplete LU (ILU), modified incomplete Cholesky (MIC), multigrid, etc. You might also use an FFT-based solver if you are willing to convert the problem to use periodic boundary conditions (see the JGT article [Stam 2001]). For example, this piece code on Robert Bridson's page implements MIC. The simplest choice is probably a Gauss-Seidel relaxation--which you already have! 
    1. Boundary conditions need special care in your linear system solver. The starter code approximates a Neumann BC with the normal derivative of pressure zero at the boundary.  Implement Neumann BCs correctly in your solver.
  7. Blur implementation bottleneck: Another serious bottleneck of the skeleton implementation is that the Gaussian blur (in Blur) is approximated using iterative local-averaging, and is currently O(N^3) for an N-by-N image and O(N) kernel radius.  Although heavy blurring is not required, the support of the effective blurring kernel (Gaussian) must be large to avoid division by zero, and to produce meaningful driving forces.  Implement something better (such as an FFT blur operator) to eliminate this bottleneck.
  8. Try one of the following if you can: This is the final project, and some of you may finish the above steps pretty fast.  In that case, you should try something a little more challenging. Here are some ideas:
    1. 3D Extension:  Given a working 2D implementation, it's only a little more work to convert the smoke solver to a full 3D implementation, e.g., replace I(i,j) with I(i,j,k), etc.  See if you can blow smoke into the shape of a 3D solid model, or make a ship sail through a smoke ring.
    2. MAC grids are better than uniform grids, so give them a try.
    3. Improved advection: Try implementing the improved advection with filtered cubic interpolation from [Fedkiw et al. 2001] and see if that improves the quality of your smoke.
    4. Implement [Fattal and Lischinski 2004]: If you're more ambitious, go ahead and implement the CFD solver described in the original paper, and see how much it improves the quality of your animations.
    5. Internal boundaries can be introduced with some modifications to your pressure solver, and can let you model more interesting domains.
    6. Rigid objects can also be introduced into the fluid with a modified linear solver, and some coupling forces. Some related graphics papers you may find useful are [Carlson et al. 2004], [Baxter and Lin 2004], and [Guendelman et al. 2005].
    7. Multiple smoke fields and colors can be introduced as in [Fattal and Lischinski 2004].  Alternately, see if you can have smoke fields of RGB colors that can transition between colors to match colored images.
    8. Fluid control: See if you can control dyes in fluids with free surfaces similar to smoke. Better still, try to control the shape of the fluid surface using control forces. You can try implementing the FLIP scheme of [Zhu and Bridson 2005] (consider this starter code).
    9. Improved rendering: Consider a ray-marching scheme (such as in [Fedkiw et al. 2001]) for improved smoke rendering.
    10. Something else!  Feel free to talk to me about any other ideas you have.
  9. One Creative Smoke-Control Artifact!: One of the best parts of computer animation is creative use of mathematics and computer programming. Create! Use your imagination to create something really interesting. Modify the starter code in any way you want. 
  10. Include videos of your work, and also please include one still image that represents your most interesting scenario.
On collaboration and academic integrity: You are allowed to collaborate on the assignments to the extent of formulating ideas as a group, and derivation of physical equations. However, you must conduct your programming and write up completely on your own, and understand what you are writing. Please also list the names of everyone that you discussed the assignment with.  (You are expected to maintain the utmost level of academic integrity in the course. Any violation of the code of academic integrity will be penalized severely.)

Hand-in using CMS: Please submit a short write-up (detailing what you did, your findings, and who you discussed the assignment with, etc.), as well as your Java implementation, a creative simulation artifact(s), and videos of anything you want me to see, etc.

References (see course homepage):
Enjoy!!!

Copyright Doug James, April 2007.