Navier Stokes Fluid Dynamics

CS 482 Lecture, Dr. Lawlor
Oftentimes we don't want to deal with watching pressure waves slowly expand and bounce around--when we push on the fluid, we want it to flow, not to just bounce around! If the pressure is defined to be uniform, we better make sure that the fluid velocities don't ever try to change the pressure. Recall that pressure change is driven from velocity divergence, often written as "div", and defined with the del operator as "-dp / dt = ∇ . v" (the "." here stands for dot like dot product). Divergence, ∇ . v, takes a vector field like velocity, and returns a scalar field like pressure difference.

From the definition of ∇,

-dp / dt = div v = ∇ . v = dot(vec3(∂ _ / ∂x, ∂ _ / ∂y, ∂ _ / ∂z),v) = ∂v/∂x + ∂v/∂y + ∂v/dz

Recall that in our little simple wave simulator,

C.z+=dt*(L.x - R.x + B.y - T.y); /* velocity difference -> height */

Again, this is just a centered-difference approximation for:

dp/dt = -(∂v/∂x + ∂v/∂y) = -∇ . v = -div v

Of course, our goal with incompressible fluids is to make dp/dt equal zero; to do this, we've got to ensure the divergence of the velocity always equals zero. One frustrating problem: by changing the velocity at one pixel we can easily force divergence to be zero there, but this will create divergence at the neighboring pixels. This means divergence cannot be eliminated locally, but must act across wide areas of the mesh.

The problem of extracting a divergence-free version of a vector field is called the Helmholtz problem. Mathematically, the cleanest way to exactly solve the Helmholtz problem is to project the vector field to frequency space with the FFT, drop the component of the resulting vectors that points toward the origin, and then inverse-FFT. This works nicely (at least on a rectangular 2D grid with power of two dimensions) but it's a lot of code to import, not terribly fast, and I haven't found a reliable fast WebGL FFT implementation.

We can more easily eliminate divergence by just computing it, and gently penalizing it. For example, our standard shallow-water velocity update dv / dt = - ∇ p pushes against pressure differences, so we can reuse this to eliminate divergence. Actually... that's what happens in the real world--at places where velocities are converging, pressure increases, pushing back on the velocities, until equilibrium is restored. The only trouble is this process takes thousands of steps to redistribute velocities across a big mesh.

We can speed up pressure's normal divergence correction process by just computing the divergence at several different spatial scales (or pixel offsets). For example, we can easily just use the texture2D bias factor to grab a velocity sample from a higher mipmap level:

        float div=0.0; /* magnitude of divergence of our vector field */
        for (float lvl=9.0;lvl>=0.0;lvl-=1.0) {
                float dist=pow(2.0,lvl)*(1.0/texture_size);
		
                L=texture2D(tex,P+vec2(-dist,0.0),lvl);
                R=texture2D(tex,P+vec2(+dist,0.0),lvl);
                T=texture2D(tex,P+vec2(0.0,+dist),lvl);
                B=texture2D(tex,P+vec2(0.0,-dist),lvl);
		
		div += L.x - R.x + B.y - T.y; /* velocity difference -> divergence */
	}

We can then change our velocities per-pixel to try to eliminate this divergence, just like it was pressure (it's often called "pseudopressure" in simulations). This goes a little bit faster if we correct the divergence as we go, like so:

        vec2 corr=vec2(0.0);
        float div=0.0;

        for (float lvl=9.0;lvl>=0.0;lvl-=1.0) {
                float dist=pow(2.0,lvl)*(1.0/texture_size);
		
                L=texture2D(tex,P+vec2(-dist,0.0),lvl);
                R=texture2D(tex,P+vec2(+dist,0.0),lvl);
                T=texture2D(tex,P+vec2(0.0,+dist),lvl);
                B=texture2D(tex,P+vec2(0.0,-dist),lvl);
        
        /* figure out how to fix velocity divergence */
                corr+=vec2(R.z-L.z,T.z-B.z);
                div+=R.x-L.x + T.y-B.y;
        }
        C.xy+=divcorr*corr; /* repair divergence using our velocity */
        C.z=divscale*div; /* write out new divergence estimate */

Note that we *don't* do a C.z += operation, like a partial differential equation, since this results in pressure waves carried by z. Instead, if there's no divergence, z has done its job and should be zero.

This seems to work quite well, despite being an Euler-style update!
Try the multigrid Navier-Stokes fluids WebGL demo.