Fluid Advection

CS 482 Lecture, Dr. Lawlor

In the Navier-Stokes equations, "v · ∇ v" is the term that represents moving fluid--this is called various things like convection (when driven by heat), advection (when driven by anything else), or just "wind" or "bulk transport".

Now, "∇ v" is the gradient of the velocity: vec3( ∂v/∂x, ∂v/∂y, ∂v/∂z). This is a little odd, because velocity is already a vector, so taking the gradient gives us a vector of vectors (a matrix, or tensor). Dotting this vector of vectors in "v · ∇ v" gives us a vector again (this is good because it gets added to dv/dt, which is also a vector). The bottom line is:

v · ∇ v = dot(v,vec3( ∂v/∂x, ∂v/∂y, ∂v/∂z)) = v.x * ∂v/∂x + v.y * ∂v/∂y + v.z * ∂v/∂z

In English, we're dotting our vector field with its own gradient. This tells us how similar the vectors are to their direction of greatest change, which in turn says how much the value of the vector will change when moving in that direction. This is how "v · ∇ v" simulates fluid transport.

Now, we could implement this by actually taking some finite difference approximation of ∂v/∂x and actually computing the vector v · ∇ v at each pixel:

   // Advection physics via explicit dot(V,grad(V)):
   vec2 gradVx=vec2(R.x-L.x,T.x-B.x);
   vec2 gradVy=vec2(R.y-L.y,T.y-B.y);
   vec2 VgV=vec2(dot(V,gradVx),dot(V,gradVy));
   V-=50.0*VgV*dt;
   Demo: fluid advection via V dot grad V

The problem with this approach is it tends to break down and go unstable if the velocity gets too big. The problem with the differential approximation to transport is that it fits a linear model to the local velocity neighborhood; moving by more than one pixel is essentially using linear extrapolation, which will amplify small waves in the mesh and explode first in zebra stripes and then NaN. (Curious? Read about the Courant-Fredrichs-Lewy (CFL) stability condition.)

Coordinate Shift Advection

On the graphics hardware, it's actually a lot easier to move stuff around onscreen by just changing your texture coordinates--normally, you look "upwind" against the flow to see where your values are coming from:

vec4 C = texture2D(srcTex,texCoords); /* center pixel, for vel estimate */
vec2 dir = C.xy; /* my velocity */
vec2 tc = texCoords - vel*dir; /* move "upwind" */
C = texture2D(srcTex,tc); /* advected source center pixel */
Demo: fluid advection via upwind coordinate shift

This also gives advection just like explicitly computing v · ∇ v, but it's more error-tolerant than v · ∇ v, because we can advect by multiple pixels in a single step. Jos Stam calls this technique "stable fluids". The hardest part about this technique is knowing where to look--the problem is your calm still pixel might be right next to an impending wave, but you won't know to advect yourself from there unless you look carefully at all your neighbors beforehand.

Worse yet, it's possible for several waves to come crashing into a single pixel at once, in which case no single source point is entirely correct. This isn't a problem if the velocity field is sufficiently smooth and well resolved, but around shocks the V dot grad V technique can be better behaved, because it explicitly includes information from your neighbors.