Simulation Class: Advection

Advection

In our previous post, we concluded our exploration of iterative matrix methods, and the Jacobi method in particular, for solving the 2D wave equation. We'll come back to that method a bit later when we explore different applications of The Poisson Equation.

We've now built up enough machinery to start building the 2D "smoke" solver promised at the outset of the class. The system we'd like to model is a fluid in a rectangular area, with a coloration and impulse that are added via external forces (the mouse, in this simple case). The coloration will decay over time, and we will use it to depict smoke, or alternatively it could be interpreted as ink in water.

Before concerning ourselves with fluid pressure, incompressibility, viscosity, vorticity, and other nuanced fluid behaviors, we will first turn our attention to Advection. Advection is just a fancy fluid-dynamics term for "the movement of stuff through a fluid due to velocity". There's another way that stuff (highly technical term) is moved through a fluid, and that is via diffusion, which we'll consider later. When the movement of stuff via fluid velocity and diffusion are considered together, it is called Convection.

Examples of advection are ink or soot particles moving through a dynamic fluid such as water or air.

While advection is outwardly one of the simpler parts of fluid dynamics to visualize and understand, it is very difficult to simulate without sacrificing either stability or detail. A very popular approach to advection, which we'll use and discuss is called Semi-Lagrangian Advection. This scheme is unconditionally stable, but sacrifices sharp details of the fluid and its properties, as we'll demonstrate.

Fluid Velocity

Firstly, we must consider what we're simulating and how to represent it in our State. For a homogeneous, incompressible fluid without an interface to a second fluid (called a free surface), the main property of the fluid that we are simulating is its velocity. This makes intuitive sense - when we think of "fluidity", we think of smooth, stream-like motion that changes in curvy ways over time. (Or at least I do). Velocity is a vector field, represented in 2D by two scalar fields, VelocityU (the fluid velocity in the x-direction) and VelocityV (the fluid velocity in the y-direction). These Scalar fields will be represented in our State exactly like the 2D scalar "height" property was represented in the previous wave equation example.

Additionally, we'll need a scalar field to represent our coloration, or soot amount. Before we even bother with looking at how fluid velocity or the soot amount changes over time, we'll first consider how to represent these fields in our simulation, and how to visualize them. (Once again returning to our main strategy). Let's start with a simple state that's reminiscent of the previous 2D Wave Equation:

float WorldSize = 10.0;
int NX = 64;
int NY = 64;
int ArraySize = NX * NY;
float DXY = WorldSize / NX;

float LX = WorldSize;
float LY = WorldSize;

int StateSize = 3;
float[][] State = new float[StateSize][ArraySize];
int StateVelU = 0;
int StateVelV = 1;
int StateSoot = 2;

To set the initial state of velocity, we'll start with a noise function just like we used with the height field in the previous example. We'll soon see that this isn't a good choice, but we need to start somewhere. In previous examples we've been writing a small noise fractal inline, but since we'll reuse that in multiple places, we'll separate it out into its own function. We'll also initialize the soot value to hatch pattern of lines, which will be easy to see move and deform.

float Fractal(float x, float y) {
  float octave_gain = 0.5;
  float lacunarity = 1.67;
  float f = 0.0;
  float amp = 1.0;
  float px = x;
  float py = y;
  for (int oct = 0; oct < 4; ++oct) {
    f += amp * snoise(px, py);
    amp *= octave_gain;
    px *= lacunarity;
    py *= lacunarity;
  }
  return f;
}

void SetFieldToFractal(int field, float gain,
    float offset_x, float offset_y,
    float pattern_size) {
  noiseSeed( 0 );
  float spacing = DXY / pattern_size;
  for (int j = 0; j < NY; ++j) {
    float noise_pos_y = offset_y + spacing * float(j);
    for (int i = 0; i < NX; ++i) {
      float noise_pos_x = offset_x + spacing * float(i);
      float d = gain * Fractal(noise_pos_x, noise_pos_y);
      State[field][IX(i,j)] = d;
    }
  }
}

// Draws a hatched grid of lines.
void SetFieldToHatch(int field, int line_width, int cell_width) {
  for (int j = 0; j < NY; ++j) {
    int j_line = (j + (line_width/2)) % cell_width;
    boolean j_in_line = j_line < line_width;
    for (int i = 0; i < NX; ++i) {
      int i_line = (i + (line_width/2)) % cell_width;
      boolean i_in_line = i_line < line_width;

      State[field][IX(i,j)] = i_in_line || j_in_line ? 1.0 : 0.0;
    }
  }
}

float VelocityInitAmplitude = WorldSize * 0.125;
void SetInitialState() {
  SetFieldToFractal(StateVelV, VelocityInitAmplitude,
      2341.17, 9911.44, WorldSize * 0.617);
  SetFieldToFractal(StateVelU, VelocityInitAmplitude,
      -8181.13, -1881.81, WorldSize * 0.617);
  SetFieldToHatch(StateSoot, 2, NX/8);
}

Note that we've removed the boundary conditions temporarily, because they'll require a bit of discussion. Given this initial stab at a fluid velocity initialization, let's work on drawing it.

Visualizing Fluid Velocity and Soot Simultaneously

How should we draw velocity and soot? There's a few different ways we could consider. Velocity is a vector field - in this case it has two components, U (the velocity of the fluid in the x direction) and V (the velocity of the fluid in the y direction). It will vary from positive to negative, and so we could try just mapping it into the red and green channels of a color image. Then we can use the blue channel to display the soot value. Here's what that would look like, genericized slightly to include a range to map the display to. We've also made some small improvements to the drawing of the label.

// Display Gamma
float DisplayGamma = 2.2;

// Draw vector field into the image. The from_low and from_high
// represent the numerical range that will be mapped into the 0-1
// space of red and green. A display gamma will then be applied.
void DrawVectorField(int field_r, int field_g, int field_b,
    float from_low_r, float from_low_g, float from_low_b,
    float from_high_r, float from_high_g, float from_high_b) {
  float r, g, b;
  StateImage.loadPixels();
  for (int j = 0; j < NY; ++j) {
    for (int i = 0; i < NX; ++i) {
      r = field_r < 0 ? 0 : State[field_r][IX(i,j)];
      g = field_g < 0 ? 0 : State[field_g][IX(i,j)];
      b = field_b < 0 ? 0 : State[field_b][IX(i,j)];

      // remap.
      r = (r - from_low_r) / (from_high_r - from_low_r);
      g = (g - from_low_g) / (from_high_g - from_low_g);
      b = (b - from_low_b) / (from_high_b - from_low_b);

      // constrain.
      r = constrain(r, 0.0, 1.0);
      g = constrain(g, 0.0, 1.0);
      b = constrain(b, 0.0, 1.0);

      // display gamma.
      r = pow(r, DisplayGamma);
      g = pow(g, DisplayGamma);
      b = pow(b, DisplayGamma);

      // Set
      StateImage.pixels[IX(i,j)] = color(r, g, b);
    }
  }
  StateImage.updatePixels();
  image(StateImage, 0, 0, width, height);
}

void draw() {
  // Draw Velocity
  DrawVectorField(StateVelU, StateVelV, StateSoot,
    -VelocityInitAmplitude, -VelocityInitAmplitude, 0.0,
    VelocityInitAmplitude, VelocityInitAmplitude, 1.0);

  // Label.
  noStroke();
  fill(1.0, 1.0, 1.0, 0.75);
  rect(0.0, 0.0, width, 40.0);
  fill( 0.0 );
  text("Fluid Velocity & Soot : RG, B", 10, 30);
}

Let's see what that looks like in Processing:

Hmmm. That's not very easy to interpret, visually. If you work with fluid solvers frequently, you can develop an ability to look at velocity in this way, but it's always non-intuitive. Plus, the simultaneous display of velocity and soot is very visually confusing. Instead, let's try drawing line segments to show the direction of the fluid at each point. This is easy to do using the "LINES" shape tools in Processing. We'll create a line starting at each cell and pointing in the direction of the velocity, with the line lengths proportional to the velocity magnitude, scaled by an arbitrary constant for ease of comprehension. In code,

// Draw velocity as line segments. We can optionally skip
// some of the points, and we can scale the length of the vectors
// accordingly.
void DrawVelocityAsLines(int U, int V,
    int step_u, int step_v,
    float dt) {
  beginShape(LINES);
  for (int j = step_v / 2; j < NY; j += step_v) {
    float pixel_start_y = PixelsPerCell * float(j);
    for (int i = step_u / 2; i < NX; i += step_u) {
      float pixel_start_x = PixelsPerCell * float(i);

      float vel_u = State[U][IX(i,j)];
      float vel_v = State[V][IX(i,j)];
      float pixel_vel_u = PixelsPerCell * State[U][IX(i,j)] / DXY;
      float pixel_vel_v = PixelsPerCell * State[V][IX(i,j)] / DXY;

      float pixel_end_x = pixel_start_x + (dt * pixel_vel_u);
      float pixel_end_y = pixel_start_y + (dt * pixel_vel_v);
      vertex(pixel_start_x, pixel_start_y);
      vertex(pixel_end_x, pixel_end_y);
    }
  }
  endShape();
}

void DrawScalarField(int field,
    float from_low, float from_high, float exponent,
    color color_low, color color_high) {
  float d;
  StateImage.loadPixels();
  for (int j = 0; j < NY; ++j) {
    for (int i = 0; i < NX; ++i) {
      d = State[field][IX(i,j)];

      d = (d - from_low)/(from_high - from_low);
      d = constrain(d, 0.0, 1.0);
      d = pow(d, exponent);

      StateImage.pixels[IX(i,j)] = lerpColor(color_low, color_high, d);
    }
  }
  StateImage.updatePixels();
  image(StateImage, 0, 0, width, height);
}

void draw() {
  background(0.5);

  // Draw soot as a purplish smoke.
  color soot_zero = color(1.0, 1.0, 1.0);
  color soot_one = color(0.3, 0.0, 0.3);
  DrawScalarField(StateSoot, 0.0, 1.0, 3.0, soot_zero, soot_one);

  // Draw Velocity as lines
  // Length gain is somewhat arbitrary.
  float dt = (1.0 / 24.0);
  stroke(1.0, 0.1, 0.1);
  DrawVelocityAsLines(StateVelU, StateVelV, 2, 2, 10.0 * dt);

  // Label.
  noStroke();
  fill(1.0, 1.0, 1.0, 0.75);
  rect(0.0, 0.0, width, 40.0);
  fill( 0.0 );
  text("Fluid Velocity & Soot : Lines, RGB", 10, 30);
}

Let's see what that looks like in Processing:

That's so much easier to grok! The length scale of the lines is somewhat arbitrary, and can be adjusted on an application-specific basis, or tied to a hotkey. Excellent! Now that we've created a velocity field in our State, and seen how to visualize it, it's time to consider our equations of motion.

The Advection Equation

Ignoring the other aspects of the fluid and just focusing on bulk motion due to velocity, we'll begin with The Advection Equation. This describes how any scalar field, represented by \(\phi\), is changed over time due just to advection by a velocity field \(\mathbf{v}\). The full equation is written in 2D as:

\[ \begin{equation} \frac{\partial \phi}{\partial t} + v_x \frac{\partial \phi}{\partial x} + v_y \frac{\partial \phi}{\partial y} = 0 \label{eq:AdvectionEquationPhi} \end{equation} \]

It is common to use a shorthand notation for the spatial part of \(\eqref{eq:AdvectionEquationPhi}\), called The Advection Operator:

\[ \begin{equation} (\mathbf{v} \cdot \nabla)\phi \equiv v_x \frac{\partial \phi}{\partial x} + v_y \frac{\partial \phi}{\partial y} \label{eq:AdvectionOperator} \end{equation} \]

With this notation, the Advection Equation becomes:

\[ \begin{equation} \frac{\partial \phi}{\partial t} + (\mathbf{v} \cdot \nabla)\phi = 0 \label{eq:AdvectionEquationConcisePhi} \end{equation} \]

This is a difficult equation to solve, and there is a tremendous amount of research and discussion of different ways to solve it. Rather than explore the instabilities of some more obvious choices, we'll jump right to the discussion of a popular, stable method.

Semi-Lagrangian Advection

Semi-Lagrangian Advection uses a technique known as the method of characteristics to provide an unconditionally stable solution to the advection equation. How it works is fairly simple. Given a velocity field, for each location in the field representing the advected quantity to be solved, "Q new", we look up the velocity at that location in the field. We use the velocity to construct a flowline from our position backwards in time, and interpolate the advected quantity, "Q old", from that location.

The above image is from Jos Stam's incredible paper, Stable Fluids, which this course borrows liberally from. Implementing that paper was how I first started to learn about the fluid simulation ideas that I'm focused on in this class.

Implementation

The good news is that implementing Semi-Lagrangian Advection of a scalar field is fairly simple. We will need a way to bi-linearly resample our 2D arrays, which is straightforward to code:

// Bi-Linear Resampling of a field Q at some point fi, fj
float BiLinearResample(int Q, float fi, float fj) {
  int i_lo = floor(fi);
  float s = fi - float(i_lo);
  i_lo = constrain(i_lo, 0, NX-1);
  int i_hi = min(i_lo + 1, NX-1);

  int j_lo = floor(fj);
  float t = fj - float(j_lo);
  j_lo = constrain(j_lo, 0, NY-1);
  int j_hi = min(j_lo + 1, NY-1);

  float q00 = State[Q][IX(i_lo,j_lo)];
  float q10 = State[Q][IX(i_hi,j_lo)];
  float q01 = State[Q][IX(i_lo,j_hi)];
  float q11 = State[Q][IX(i_hi,j_hi)];
  return lerp(lerp(q00, q10, s), lerp(q01, q11, s), t);
}

The above code handles the boundary conditions of the four walls (edges) of our simulation by clamping indices. If we did have additional internal boundary conditions, such as solid collision regions, we'd have to take more care in resampling the fields, something we'll look into in later classes. But for now, the above works great.

At first, we're just going to advect the Soot field by the velocity field, which means we'll be computing a StateSoot from a StatePrevSoot. We already have some Copy & Fill tools from the Wave Equation 2D example, so we'll just need to add an extra field to the State and implement a SwapSoot function. Here are those changes:

int StateSize = 4;
float[][] State = new float[StateSize][ArraySize];
int StateVelU = 0;
int StateVelV = 1;
int StateSoot = 2;
int StatePrevSoot = 3;

...

void SetInitialState() {
  ...
  SetFieldToHatch(StateSoot, 2, NX/8);
  CopyField(StateSoot, StatePrevSoot);
}

...

void SwapSoot() {
  int tmp = StateSootPrev;
  StateSootPrev = StateSoot;
  StateSoot = tmp;
}

A first-order implementation of Semi-Lagrangian Advection is then quite straightforward. At each point, we look up the velocity. We convert the velocity into cell-space, and look backwards along that velocity from the position we're considering by the timestep dt. At this new location, we bilinearly resample the previous field to produce the new, advected field value. In code:

// Semi-Lagrangian Advection of Q by U & V
void SemiLagrangianAdvectFirstOrder(int NEW_Q, int OLD_Q,
    int U, int V, float dt) {
  for (int j = 0; j < NY; ++j) {
    for (int i = 0; i < NX; ++i) {
      // velocity in cell space
      float cell_vel_u = State[U][IX(i,j)] / DXY;
      float cell_vel_v = State[V][IX(i,j)] / DXY;

      // Move point backwards in time by dt
      float fi = float(i) - dt * cell_vel_u;
      float fj = float(j) - dt * cell_vel_v;

      State[NEW_Q][IX(i,j)] = BiLinearResample(OLD_Q, fi, fj);
    }
  }
}

Now, all we need is a time step that advects the soot and swaps the soot fields, and we are done! Here's that, with the corresponding addition to the "draw" function. Note that I've been adding "noLoop" to the setup functions of the static sims so far, to improve performance. Just make sure that, if you've added the noLoop() function, you remove it once you turn on these dynamics.

void TimeStep(float dt) {
  SwapSoot();
  SemiLagrangianAdvect(StateSoot, StateSootPrev,
      StateVelU, StateVelV, dt);
}
...
void draw() {
  TimeStep(1.0/24.0);
  ...
}

I also like having the reset capability, which just flashes things back to the initial state, bound to the 'r' key. Here's code which adds that:

// Reset function. If the key 'r' is released in the display,
// copy the initial state to the state.
void keyReleased() {
    if ( key == 114 ) {
        SetInitialState();
    }
}

That's it! Let's see how it looks in Processing: (Hit 'r' to reset)

Advection of Velocity by Velocity

It's not just the soot that's advected by the velocity of the fluid, though. The fluid mass itself is moving, and so we have to advect the velocity by itself! In equation form, the self advection of the velocity by itself looks like this:

\[ \begin{equation} (\mathbf{v} \cdot \nabla) \mathbf{v} = 0 \label{eq:SelfAdvection} \end{equation} \]

I always found this mildly paradoxical, but fortunately it's easy to do! We'll need Prev versions of the velocity in the State:

int StateSize = 6;
float[][] State = new float[StateSize][ArraySize];
int StateVelU = 0;
int StatePrevVelU = 1;
int StateVelV = 2;
int StatePrevVelV = 3;
int StateSoot = 4;
int StatePrevSoot = 5;

As with soot, we'll need proper initialization, and a swap function:

void SetInitialState() {
  ...
  CopyField(StateVelU, StatePrevVelU);
  CopyField(StateVelV, StatePrevVelV);
}

...

void SwapVel() {
  int tmp = StatePrevVelU;
  StatePrevVelU = StateVelU;
  StateVelU = tmp;

  tmp = StatePrevVelV;
  StatePrevVelV = StateVelV;
  StateVelV = tmp;
}

The time step then just has self-advection added, AFTER the advection of the soot:

void TimeStep(float dt) {
  ++TotalTimeSteps;
  StateCurrentTime += dt;
  SwapSoot();
  SwapVel();
  SemiLagrangianAdvect(StateSoot, StatePrevSoot,
      StatePrevVelU, StatePrevVelV, dt);
  SemiLagrangianAdvect(StateVelU, StatePrevVelU,
      StatePrevVelU, StatePrevVelV, dt);
  SemiLagrangianAdvect(StateVelV, StatePrevVelV,
      StatePrevVelU, StatePrevVelV, dt);
}

And that's all we have to add for self-advection - let's have a look: (Hit 'r' to reset)

It's fairly obvious that our lack of boundary conditions and our questionable initial velocities are producing somewhat strange behavior over time, but we can clearly see that the velocity and soot are being moved along the streamlines, and that we aren't getting any instabilities. In the next classes, we'll consider our fluid boundary conditions and start adding extra fluid behavior such as incompressibility.

Download

All of the files associated with this class are available via a public github repository, found here: https://github.com/blackencino/SimClass