The project 2DStableFluids.zip contains the following components:

- A simple 2D visualization code using windows GDI API (since we are only dealing with 2D). One region of the client windows draws a 60x60 grid fluid domain, and the rest displays usage information.
- A sparse linear matrix data structure and a simple iterative solver. A.solve(x, b, tolerance, max_iteration) solves Ax=b.
- Simple user interactions are handled as done in previous assignments, including a few pre-defined shortcut keys and mouse event handlers. Left button for adding smoke, and right button for stirring the fluid.
- A fluid solver that updates two things when the timer is set: velocity and smoke density.
- A simple 2D vector data structure to store velocity is also included.

Most of the display functions and interaction functions are in ChildView.cpp, which does not need to be changed except for certain choices of extra credits.

What you need to modify will be mostly within the CFluidSolver class (FluidSolver.cpp and .h)

Run the program and start animation by pressing 'z'. You can use left button to inject smoke anywhere in the fluid domain.

If you want to see the effects of diffusion coefficient, search for diffusion_coef and modify its value. It is located where the diffusion equation is constructed.

Next, we can observe the change of density as a function of spatial location driven by the motion of a constant fluid flow by REMOVING the following lines from density_advection();

for (int i=0; i<size; i++) density[i]=density_source[i]; return;

Now the injected smoke through left button of the mouse will be carried by a velocity field moving towards the lower right direction. Remember to start the animation first.

Overview: There are several step to get the velocity field updated, but one important step (projection) is implemented for you, and a couple of them are omitted (saved for extra credit).

- Advect (the velocity field carried away by the velocity field itself). This is your main task.
- Incorporate user input by adding it to the current velocity.
- Include buoyancy force computed based on density function (extra credit).
- Diffuse (extra credit).
- Project to divergence-free field.

First thing to do: In update() function of the fluid solver (called within the timer), uncomment

updateVelocity();Now the velocity field is being updated together with the density field.

Next: Modify reset() function (used in the beginning or when the user presses 'r' to reset), so that the initial velocity is initialized to 0.

velocity[i] = vec2(0.,0.);

Press 'v' to show the velocity field. Now you can see the velocity changed when you move the mouse while holding right button down: (again, don't forget to press 'z' first).

This certainly doesn't look like fluid motion at all. One of the reasons is that we did not make sure that the fluid is approximate incompressible.

Uncomment the pressure projection step in updateVelocity():

projection();

Now, the user input is always turned into essentially the collection of vortices, the type of velocity field that does not compress the fluid anywhere. You can add some smoke to see the effects of this divergence-free velocity field.

The most important step in fluid simulation is the advection. It is also the part that produces a lot of troubles in numerical simulations.

To create a reasonable simulation as done in Jos Stam's paper Stable Fluids, this step can be implemented by backtracking the grid points:

1. start from coordinates (i, j), travel along the direction -velocity(i,j) by the time step size h.

2. find out the four closest grid points (i - v.x(i,j) *h, j - v.y(i,j) * h) as in texture mapping.

3. interpolate the velocities to get that velocity.

That would be the velocity for (i,j) since after time h, that particle is
supposed to reach (i,j). Note that before the advection function is called,
array *velocity* contains the previous velocity. So it should be used as input,
and the array *advected_velocity* as
output. You are supposed to REPLACE the following lines with your own code in
velocity_advection():

for (int i = 1; i < n-1; i++) for (int j=1; j < n-1; j++) { advected_velocity[i+j*n] = velocity[i+j*n]; }

When it is done, you should be able to observe proper fluid-like behavior, e.g. the two opposite vortices on each side of user mouse stroke (holding right button down).

If you are having trouble, take a look at how density advection is implemented. It follows the same idea.

The main step is fairly simple, so the extra credits can probably be more fun, you may choose any of the follows or come up with something reasonable:

- Add buoyancy force proportional to the density function along negative y direction.
- Add diffusion for the velocity. This should be similar to the density diffusion, but with a different coefficient.
- Create a Pacman world: the top and bottom are connected, and the left side and right side are connected. You will need to change the matrices as well as the advection step. It is hard, and worthy of 2 points.
- Add pseudo-color for any of the DYNAMIC variables.
- Create a Stokes or Darcy flow.
- Use of staggered grid. (hard, so 2 extra points)

When you complete the individual project, proceed to the page on How to submit or email Visual Studio projects.