Arvutigraafika - Kursused - Arvutiteaduse instituut

Interactive 2D Fluid Simulation

Henrik Peegel, Jaakob-Jaan Avvo, Ander Pavlov

itch.io link: https://vapsik.itch.io/fluid-simulation

Repository link: https://github.com/vapsik/Interactive_Fluid_Simulations

Introduction

In this project we created an interactive 2D fluid simulation that enacts the Navier-Stokes' equations for incompressible fluids. The core algorithm used to simulate the fluids is Stable Fluids based on Jos Stam's seminal work, a 2001 article in the field of real-time fluid simulation that has been extensively used in video games since its inception. In this algorithm, a simple incompressible (constant density) fluid is represented in a rectangular grid as a field of coupled velocity vectors and pressure scalars that will dictate the movement of the dye field representing the color of the fluid or smoke being simulated. The method relies on semi-Lagrangian advection, which allows for stable, large time-step simulation without numerical divergence, combined with a projection step that enforces incompressibility, which makes the velocity field divergence-free by iteratively solving the Poisson's equation. Diffusion, external forces, and boundary interactions are incorporated in a computationally efficient time evolution scheme such that larger time steps can be used without losing physical accuracy or stability.

Navier-Stokes' Equations

The N-S equations are the following {$$\begin{align} & \frac{\partial \vec{u}}{\partial t} = - \frac{\nabla p}{\rho} - (\vec{u} \cdot \nabla)\vec{u} + \nu \nabla^{2} \vec{u} + \frac{\vec{F}_{\text{external}}}{\rho}, \\&\nabla \cdot \vec{u} = 0, \end{align}$$} where the terms are the following:

{$\vec{x}$} - let that be the position vector which in 2D is just {$\vec{x} = (x,y)$}.
{$\rho$} - constant density of the fluid,
{$\nu$} - the viscosity constant,
{$\vec{u} = \vec{u}(\vec{x}, t)$} - the vector field of velocities,
{$p = p(\vec{x},t)$} - the scalar-valued pressure field (implicit in the solver),
{$\vec{F}_{\text{external}}$} - vector field of external forces such as gravity, buoyant forces or other forces.

Because solving for viscosity is computationally expensive and it can be cheaply imitated by applying a time dependent decay on the dye color, this project will simulate a 0 viscosity fluid and apply exponential decay on the dye colours.

The zero-divergence requirement {$\nabla \cdot \vec{u} = 0$} means that in the final solution of each iteration step the sum of fluid components entering each discrete cell will be equal to the components exiting the same cell. This kind of requirement will produce the swirliness of incompressible fluids in the simulation

Discrete Representation of the Vector Fields

Let's say we want to discretely evaluate (solve) in a rectangular spatial domain {$R$} and time {$t>0$} some type of function {$f(x,y,t)$} which is governed by some partial differential equation. Let {$R$} be defined as {$$\begin{align} R = \{(x,y) \in \mathbb{R}| \ 0<x<a, 0<y<a\}, \quad a = \text{const.} >0. \end{align}$$} In finite-difference methods for solving PDE-s in 2D, the spatial domain is usually discretized by some fixed smallest spatial step for all coordinates as {$\Delta x = h_{x} >0, \Delta y = h_{y} >0$}. Time is also discretized to consist of smallest possible timesteps {$\Delta t >0$}. To assure the integer number of smallest spaces in a domain {$R$}, the smallest steps are defined by the integer numbers of cells in each direction. Let that be an integer constant {$N>0$} for both axes. Then the following can be stated about the spatial discretization:

{$$\begin{align} &\Delta x = \Delta y = \frac{a}{N} \end{align}$$} and the position vector {$\vec{x} =(x,y)$} coordinates in the domain can be expressed for some integer {$0\leq i\leq N$} as {$$\begin{align} &x_{i} = i\Delta x, \\ &y_{i} = i\Delta y, \end{align}$$} therefore there are {$N+1$} nodes from {$n=0$} to {$n=N$} where a scalar field {$f(x,y) =f(x_{i}, y_{j})$} can be evaluated therefore all the {$(N+1)(N+1)$} values of {$f(x_{i},y_{i})$} can be represented as a matrix {$$\begin{align} f_{ij}=f(x_{j},y_{i}) = \begin{pmatrix} f(x_{0},y_{0}) & f(x_{1},y_{0}) & \dots&f(x_{n},y_{0}) \\ f(x_{0},y_{1}) & f(x_{1},y_{1}) & \dots & f(x_{n},y_{0}) \\ \vdots&\vdots&\ddots&\vdots \\ f(x_{0},y_{n}) & f(x_{1},y_{n}) & \dots & f(x_{n},y_{n}) \end{pmatrix}. \end{align}$$}

Staggered Discretization of Velocities

In computational fluid dynamics, the pressure is (implicitly) discretized. However, to ensure the numerical ability to create 0-divergence, a staggered grid will be used to map velocity field components on the edges of cells that surround some pressure values. Therefore, if pressure for example takes on a matrix {$(p)_{i,j}$} with dimensions {$(N_{y}) \times (N_{x})$} in the form {$$(p)_{i,j} = \begin{pmatrix} p_{00}&p_{01}&\dots&p_{0,N_{x}-2}&p_{0,N_{x}-1} \\ p_{10}&p_{11}&\dots&p_{1,N_{x}-2}&p_{1,N_{x}-1} \\ \dots&\dots&\ddots&\vdots \\ p_{N_{y},0}&p_{N_{y},1}&\dots&p_{N_{y}-2,N_{x}-2}&p_{N_{y}-2,N_{x}-1} \\ p_{N_{y}-1,0}&p_{N_{y}-1,1}&\dots&p_{N_{y}-1,N_{x}-2}&p_{N_{y}-1,N_{x}-1} \end{pmatrix} \in \mathbb{R}^{N_{y}\times N_{x}},$$} then the corresponding velocity fields {$(u_{x})_{k,l}$} and {$(u_{y})_{m,n}$} will be created such that each pressure value {$p_{ij} = (p)_{i,j}$} will have two vertical neighbors that correspond to the vertical velocities {$(u_{y})_{i,j}$} its cell and {$(u_{y})_{i+1,j}$} below, and two horizontal neighbors that correspond to the horizontal velocities {$(u_{x})_{i,j}$} from the left and {$(u_{x})_{i,j+1}$} from the right.

Therefore {$$\begin{align} &(u_{x})_{k,l} \in \mathbb{R}^{N_{y}\times (N_{x}+1)}, \\ &(u_{y})_{m,n} \in \mathbb{R}^{(N_{y}+1)\times N_{x}}. \end{align}$$}

Visualization of the staggered grid stencil Pressure scalar field and dye scalar field (matrices) have the same shape.

Computation Implementation

Since the update loop of the simulation includes computations up to time complexity {$\mathcal{O}(kN^2)$}, where {$N$} is the dimension of the grid and {$k \approx 10^1...10^2$} the number of Poisson iterations, the implementation will hold discrete fields info in data textures. The data textures will be passed into compute shaders each frame, where fast parallelized computation will take place. HLSL compute shaders are used alongside with C# scripts for the fluid interactive simulations in Unity game engine.

Update Loop

The simulation will start with some initial conditions {$\vec{u}(\vec{x},0) = \vec{u}_{0}$}. Let the solution of the algorithm after some time {$t$} be {$\vec{u}(\vec{x}, t) = \vec{w}_{0}(\vec{x})$}.Then the update loop of the Stable Fluids simulation algorithm proceeds as follows:

{$$\vec{w}_{0}(\vec{x}) \overbrace{ \to }^{ \text{add force} } \vec{w}_{1}(\vec{x}) \overbrace{ \to }^{ \text{advect} } \vec{w}_{2}(\vec{x})\overbrace{ \to }^{ \text{remove divergence}}\vec{w}_{3}(\vec{x}).$$}

This procedure is calculated each frame/tick of the simulation.

Advection of Velocities Using Bilinear Interpolation

One of the most important techniques introduced in Jos Stam’s Stable Fluids paper is the unconditionally stable semi-Lagrangian advection of the velocity field. For the advection term of the Navier-Stokes' {$$\frac{\partial\vec{u}}{\partial t} = - (\vec{u}\cdot \nabla)\vec{u}$$} instead of calculating the difficult conditionally stable derivatives using finite differences, the Semi-Lagrangian approaches the advection of fluids holistically by finding the new velocity components by bactracking and resampling using only bilinear interpolation (in 3D it would be trilinear interpolation) and an ODE solver step such as Euler or Runge-Kutta.

Visual scheme for advecting a left x-component of an ij-th cell.

Parallelized Divergence Clearance on Compute Shader

The swirliness of non-compressible fluids comes from the 0-divergence constraint {$$\nabla \cdot \vec{u} = 0.$$}

After adding the forces and advecting the velocities, the resulting velocity vector field most probably isn't 0-divergent. Therefore, a pressure field must be found such that its gradient removes the divergence from the velocity field. So far in the update the pressures have not existed (are implicit). Using vector calculus, it can be seen, that the suitable divergence clearing pressure must be the solution to a Poisson's equation.

We are left with the next part of NS that applies the pressure gradient to the field: {$$\begin{align} &\frac{\partial \vec{u}}{\partial t} = - \frac{1}{\rho} \nabla p, \\ & \nabla \cdot \vec{u} = 0. \end{align}$$}

Since, pressure isn't known, it must be algebraically solved for such that the zero-divergence condition is satisfied. Currently the velocity field solution from the earlier step {$\mathbf{u}^{\text{old}}$} most certainly is not divergence free so we want to apply the gradient such that the {$\mathbf{u}^{\text{new}}$} is divergence free. The finite difference version of the first term is:

{$$\begin{align} &\frac{\mathbf{u}^{\text{new}} -\mathbf{u}^{\text{old}}}{\Delta t} = -\frac{1}{\rho} \nabla p, \\ &\nabla \cdot \mathbf{u}^{\text{new}} = 0, \\ \implies& \mathbf{u}^{\text{new}} - \mathbf{u}^{\text{old}} = -\frac{\Delta t}{\rho} \nabla p, \Big| \text{ apply divergence }\nabla \cdot() \\ \implies & \nabla^{2} p = \frac{\rho}{\Delta t} \nabla \cdot \mathbf{u}^{\text{old}}. \end{align}$$} The equation {$$\nabla^{2} p = \frac{\rho}{\Delta t} \nabla \cdot \mathbf{u}^{\text{old}}$$} is the Poisson equation for the implicit pressure that can be solved for pressure with different iterative methods. Most common and the easiest method to implement (which we implement in this project) is the Jacobi Iteration which in this case iteratively sets each {$(p)_{i,j}$} value across the grid to be the sum of the average of its four neighbors minus the value of the RHS of the pressure Poisson. More advanced methods that converge with fewer steps are Gauss-Seidel, Successive Over-Relaxation, Conjugate Gradient and the best but also the most difficult to implement is Multigrid.

Advection of Dye Field & Color Blending

After the velocities are updated, the color dye textures will be updated with an advection procedure similar to the velocities. The center velocities will be used for finding the backtrace location in the texture (must be clamped or reflected, if outside the dye's index space). Using the resulting locations, new colors will be found by bilinearly blending the texture on that position and placed onto the new location.