Simulation  Stable fluids
Path of the code: [11_stable_fluids]
The objective of this scene is to simulate the "stable fluid" model (in 2D) proposed by Jos Stam in [Stable Fluids. Jos Stam. SIGGRAPH 1999]. And later described with more algorithmic details in [RealTime Fluid Dynamics for Games. Jos Stam. Game Developers Conference, 2003].
In the current state, the code encodes a velocity grid (the velocity of the fluid) and a density grid (the color/image which is visualized). The values of these grid are only advected through time, and your objective is to complete the three steps of the stable fluids method, namely
  Diffusion (function diffuse)
  Advection (function advect)
  Projection of the velocity to a divergence free vector field (function project)
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 unconditional stability: \((1+4\Delta t\,\mu)\,f ^ {t} _ {x,y}=f ^ {t\Delta t} _ {x,y}+ \mu\,\Delta t\,\left(f ^ {t} _ {x+1,y}+f ^ {t} _ {x1,y}+f ^ {t} _ {x,y+1}+f ^ {t} _ {x,y1}\right)\) > Implement the numerical solution of this equation using Gauss Seidel iterations (consider around 15 iterations). The function to complete use template arguments and is defined in the filesimulation.hpp
.
Note:
  For each Gauss Seidel iteration, you may set the values of the grid boundary using the precoded function set_boundary.
  The Gauss Seidel method consists in iterating over

 \(\forall (x,y), \;\displaystyle (f_{x,y}^{t})^{\text{new}} = \frac{1}{1+4\,\mu\,\Delta t}\left(f ^ {t\Delta t} _ {x,y}+\mu\,\Delta t\left(f ^ {t} _ {x+1,y}+f ^ {t} _ {x1,y}+f ^ {t} _ {x,y+1}+f ^ {t} _ {x,y1}\right)\right)\)
 \(\displaystyle \forall (x,y),\;f_{x,y}^{t} \leftarrow (f_{x,y}^{t})^{\text{new}}\)
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. 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 substeps. 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\)