What is wrong with my bilinear interpolation in Navier-Stokes model?

Question

I'm modeling density of a fluid with a color map using numerical methods to compute the Navier-Stokes equations. This is the code which computes the advected velocity and density based on the velocities and densities during the previous timestep.

Notes: u and d are the velocity and density grids (each with two time frames: k and k + 1). The dt is a fixed timestep. The delta variable is dy = dx (square grid).

// Computes the advection terms
private void Advect(int i, int j)
{                 
  // Calculate corner indices
  int i0 = Mathf.Clamp(Mathf.FloorToInt(((float)i * delta - u[i, j, k].x * dt) / delta), 0, n - 1);
  int j0 = Mathf.Clamp(Mathf.FloorToInt(((float)j * delta - u[i, j, k].y * dt) / delta), 0, n - 1);
  int i1 = Mathf.Clamp(i0 + 1, 0, n - 1);
  int j1 = Mathf.Clamp(j0 + 1, 0, n - 1);

  // Bilinear interpolation parameters
  Vector2 L   = new Vector2(i0 * delta - delta / 2.0f, j0 * delta - delta / 2);
  Vector2 U   = new Vector2(i0 * delta + delta / 2.0f, j0 * delta + delta / 2);
  Vector2 pos = new Vector2(i * delta - u[i, j, k].x, j * delta - u[i, j, k].y);
  Vector4 a   = new Vector4(Vector2.Distance(pos, L),                     // Bottom left
                            Vector2.Distance(pos, new Vector2(L.x, U.y)), // Bottom right
                            Vector2.Distance(pos, new Vector2(L.y, U.x)), // Top left
                            Vector2.Distance(pos, U));                    // Top right

  // Interpolate velocity
  u1[i, j, k + 1] = a[0] * u[i0, j0, k] +
                    a[1] * u[i1, j0, k] + 
                    a[2] * u[i0, j1, k] +
                    a[3] * u[i1, j1, k] /
                   (a[0] + a[1] + a[2] + a[3]);

  // Interpolate density
  d[i, j, k + 1] = a[0] * d[i0, j0, k] +
                   a[1] * d[i1, j0, k] + 
                   a[2] * d[i0, j1, k] +
                   a[3] * d[i1, j1, k] /
                  (a[0] + a[1] + a[2] + a[3]);
}

Given no external force on the system, and an initial condition with a square of dense fluid in the center, and 0 density everywhere else, my program outputs this model. The expected behavior was for this block of dense fluid to diffuse slowly into the less dense fluid around it.

Now, originally this computation simply averaged the velocities/densities of the neighboring gridpoints instead of using bilinear interpolation. This resulted in what appeared to be correct fluid behavior, but not quite. The problem is that the only force applied here is an upward vertical force in the center. The fluid is leaking toward the bottom left. This, my professor believed was an issue due to truncation caused by averaging. You can see this leak more clearly if I simply don't apply any external forces. Hence, the switch to bilinear interpolation.

So I'm confused why the switch to bilinear interpolation broke the behavior. Due to this, I believe their must be an issue with my implementation of bilinear interpolation.

Answer 1

Is the code you posted exactly what you used? I am wondering a bit about where you interpolated using

d[i, j, k + 1] = a[0] * d[i0, j0, k] +
               a[1] * d[i1, j0, k] + 
               a[2] * d[i0, j1, k] +
               a[3] * d[i1, j1, k] /
              (a[0] + a[1] + a[2] + a[3]);

Are you sure it shouldn't be (note parentheses)

d[i, j, k + 1] = ( a[0] * d[i0, j0, k] +
               a[1] * d[i1, j0, k] + 
               a[2] * d[i0, j1, k] +
               a[3] * d[i1, j1, k] ) /
              (a[0] + a[1] + a[2] + a[3]);