Incompressible Navier-Stokes via Multigrid

CS 493 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 bias=maxbias;bias>=0.0;bias--) {
                float far=pow(2.0,bias);
		
		R = texture2D(srcTex,tc+far*srcDirs[0],bias); /* 4 directions */
		T = texture2D(srcTex,tc+far*srcDirs[1],bias);
		L = texture2D(srcTex,tc+far*srcDirs[2],bias);
		B = texture2D(srcTex,tc+far*srcDirs[3],bias);
		
		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 bias=maxbias;bias>=0.0;bias--) {
                float far=pow(2.0,bias);
                R = texture2D(srcTex,tc+far*srcDirs[0],bias); /* 4 directions */
                T = texture2D(srcTex,tc+far*srcDirs[1],bias);
                L = texture2D(srcTex,tc+far*srcDirs[2],bias);
                B = texture2D(srcTex,tc+far*srcDirs[3],bias);
        
        /* 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 WebGL demo.

Adding boundary conditions to these simulations is usually fairly easy. For example, we can just force the fluid velocity to be zero over land, and the fluid will flow around it.
   Try fluid flowing around disk island demo.

But there's a problem with multigrid and boundary conditions. If we ignore the boundary conditions when sampling the coarser grids, it's easy to get nonphysical behavior where effects like our divergence correction leap across what should be solid boundaries. In this demo, the grid is divided into square zones, so it should act like many disconnected simulations, but unless we modify the multigrid sampling, information propagates between the zones.
   Try multigrid cells demo.

In the general case, with arbitrary inlets surrounded by elongated peninsulas, multigrid can actually be pretty hard to do correctly.