## Simulation - Stable fluids

*Path of the code: [09_stable_fluids]*

- - Diffusion (function
*diffuse*) - - Advection (function
*advect*) - - Projection of the velocity to a divergence free vector field (function
*project*)

*simulate*. Note that diffusion and advection steps are applied to both velocity and density field, while the projection is only applied to the velocity. Clicking on the left button of the mouse and dragging it allows to set the speed on the corresponding grid cell (but this speed doesn't evolve yet). We will not consider any other external forces applied on the fluid.

## Diffusion

The function*diffuse*is set to solve one step of the equation \(\frac{\partial f}{\partial t}=\mu\triangle f\), where \(f\) can be either the velocity, or a density. To handle this generic input, the function is encoded to received a 2D grid of template values. The function is also defined in the header file for this reason. The equation is solved numerically using an implicit scheme in order to ensure inconditional stability: \((1+4\Delta t\,\mu)\,f ^ {k+1} _ {x,y}=f ^ {k} _ {x,y}+\Delta t\,\mu\left(f ^ {k+1} _ {x+1,y}+f ^ {k+1} _ {x-1,y}+f ^ {k+1} _ {x,y+1}+f ^ {k+1} _ {x,y-1}\right)\)

**> Implement the numerical solution of this equation using Gauss Seidel iterations (consider around 15 iterations).**

*Note*:

- -
*even if \(f\) is a template parameter, you can use its value as any variable (ex. f(x,y), f(x-1,y), etc.)*

- -
*for each Gauss Seidel iteration, you may set the values of the grid boundary using the pre-coded function set_boundary.*

*Example of diffusion applied to the density with \(\mu=0.02\)*

## Advection

The function*advect*is set to displace a field \(f\) along a prescribed velocity \(v\). \(f\) can itself be a density, or the velocity itself. To ensure stable advection, backtracing can be use in setting the value \(f(p,t)\) where \(p\) is a grid cell position, and \(t\) the current time from the interpolated value of f at previous time step at position \(f(p-\Delta t\,v(p,t), t-\Delta t)\). The advection function is already implemented in the function

*advect*using bilinear interpolation at position \(p-\Delta t\,v(p,t)\). Pay attention to not fetch values outside the boundaries of the grid. Once implemented, you should be able to create visible motion from your mouse trajectory in your density. This motion should be smooth, but doesn't ensure yet incompressible behavior.

## Divergence free projection

The last step ensuring incompressibility is to constraint the velocity field to be divergence free. This step need only to be applied on the velocity field (and not on density), and consists in three sub-steps.- 1. Compute the divergence of the input velocity \(v_0\)
- 2. Compute a gradient field \(q\) satisfying Poisson equation \(\triangle q = div(v)\)
- 3. Compute divergent free vector field as \(v=v_0-\nabla\,q\)

*set_boundary*function.

**> Implement these three steps in the**to reach a divergence free vector field that you can interactively control.

*divergence_free*function