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

Simulation Class: The Wave Equation in 2D

The Wave Equation in 2D

The 1D wave equation solution from the previous post is fun to interact with, and the logical next step is to extend the solver to 2D. It turns out that this is almost trivially simple, with most of the work going into making adjustments to display and interaction with the state arrays. As we'll see, we'll make almost no changes to the physics whatsoever - this will be an ongoing theme in this class, that simulation programming is mostly about input, output, iteration over arrays, and display, and involves very little actual physics (in terms of lines of code).

The first thing we'll do is change the state array from a 1D array to a 2D array. We'll also change notation so that X & Y represent the dimensions of the array, and Z represents the vertical dimension expressed by heights in the array. Our state initialization now looks like this (showing only the lines that have changed):

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

float LX = WorldSize;
float LY = WorldSize;
float LZ = WorldSize / 2.0;

We'll keep the State array declaration exactly the same, which means that we're still declaring the individual fields themselves as one-dimensional arrays. In order to index these 1D arrays from 2D coordinates, we'll need an indexing function:

int IX( int i, int j )
{
    return i + NX*j;
}

Next we implement the initialization of the state arrays, which is almost exactly the same as in the 1D case, except using a 2D iteration and 2D noise patterns. Otherwise, it's basically identical:

void SetInitialState() {
    noiseSeed( 0 );
    for (int j = 0; j < NY; ++j) {
        for (int i = 0; i < NX; ++i) {
            float worldX = 2341.17 + DXY * ( float )i;
            float worldY = 9911.44 + DXY * ( float )j;

            float n = 0.5 * snoise( worldX * 0.0625, worldY * 0.0625 ) +
                0.4 * snoise( worldX * 0.125, worldY * 0.125 ) +
                0.3 * snoise( worldX * 0.25, worldY * 0.25 ) +
                0.2 * snoise( worldX * 0.5, worldY * 0.5 );
            n = 1.0 - abs(n);
            n = (2.0 * n) - 1.0;

            State[StateHeight][IX(i,j)] = n;

            State[StateVel][IX(i,j)] = 0.0;
        }
    }
    EnforceHeightBoundaryConditions( StateHeight );
    EnforceNeumannBoundaryConditions( StateVel );
    StateCurrentTime = 0.0;

    CopyArray( StateHeight, StateHeightPrev );
    CopyArray( StateVel, StateVelPrev );
    InputActive = false;
}

Boundary Conditions in 2D

Our boundary conditions used to only be concerned with the first and last points of the linear 1D Array. Now we have to be concerned with the entire boundary, on all four edges. The Dirichlet and Neumann boundary functions are changed thus:

void EnforceDirichletBoundaryConditions( int io_a ) {
    for (int j = 0; j < NY; ++j) {
        if (j == 0 || j == (NY-1)) {
            for (int i = 0; i < NX; ++i) {
                State[io_a][IX(i,j)] = 0.0;
            }
        } else {
            State[io_a][IX(0,j)] = 0.0;
            State[io_a][IX(NX-1,j)] = 0.0;
        }
    }
}

void EnforceNeumannBoundaryConditions( int io_v ) {
    for (int j = 0; j < NY; ++j) {
        if (j == 0) {
            for (int i = 0; i < NX; ++i) {
                State[io_v][IX(i,0)] = State[io_v][IX(i,1)];
            }
        } else if (j == (NY-1)) {
            for (int i = 0; i < NX; ++i) {
                State[io_v][IX(i,NY-1)] = State[io_v][IX(i,NY-2)];
            }
        }

        State[io_v][IX(0,j)] = State[io_v][IX(1,j)];
        State[io_v][IX(NX-1,j)] = State[io_v][IX(NX-2,j)];
    }
}

Mouse Input in 2D

Input gathering is a fairly simple process as well - we didn't explicitly go over this in the previous example, so we'll provide a quick explanation here. In the top of the code, we introduce a small amount of input state:

boolean InputActive = false;
int InputIndexX = 0;
int InputIndexY = 0;
float InputHeight = 0;

And then we've got a function that gathers input when the mouse button is pressed:

void GetInput() {
    InputActive = false;
    if ( mousePressed && mouseButton == LEFT ) {
        float mouseCellX = mouseX / ( ( float )PixelsPerCell );
        float mouseCellY = mouseY / ( ( float )PixelsPerCell );

        int iX = ( int )floor( mouseCellX + 0.5 );
        int iY = ( int )floor( mouseCellY + 0.5 );
        if ( iX > 0 && iX < NX-1 && iY > 0 && iY < NY-1 ) {
            InputIndexX = iX;
            InputIndexY = iY;
            InputHeight = 1.5;
            InputActive = true;
        }
    }
}

This input data is used in only one place - when the height field boundary conditions are enforced, we set the input height explicitly (a Dirichlet condition). This is implemented via an explicit boundary handler for height:

void EnforceHeightBoundaryConditions( int io_h ) {
    EnforceNeumannBoundaryConditions(io_h);
    if ( InputActive ) {
        State[io_h][IX(InputIndexX, InputIndexY)] = 1.0;
    }
}

Displaying the 2D Results

Before we worry about any of the physics - we'll just work on getting the drawing of the data working. This is a revisiting of our basic philosophy and workflow - define your State, figure out how to display your State, and then work on your physics!

Processing offers a 2D image type that we can use to write to the pixels directly. Near the top of the file, we'll globally create this image:

PImage StateImage = createImage( NX, NY, RGB );

And then in the display function, we'll convert our height data into pixel colors. This is a somewhat aesthetic decision - our simulation data will range from -1.0 to 1.0 in height, generally, so we'll remap those values so that smaller numbers are drawn in dark blue and bigger numbers are drawn in white. Here's the draw-image code. We'll try this out temporarily without a time step just to make sure that the display works.

// Draw height field into the image.
void DrawHeightField( int i_field ) {
    float pixr, pixg, pixb;
    float d;
    color dark_blue = color(0.01, 0.01, 0.2);
    color light_blue = color(0.9, 0.9, 1.0);
    StateImage.loadPixels();
    for (int j = 0; j < NY; ++j) {
        for (int i = 0; i < NX; ++i) {
            float d = State[i_field][IX(i,j)];

            // Renormalize
            d = constrain( (d + 1.0) / 2.0, 0.0, 1.0 );

            // Add contrast
            d = pow( d, 8 );

            // Interpolate color.

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

void draw() {
    background( 0.5 );

    // Simulate
    //GetInput();
    //TimeStepRK4( 1.0 / 24.0 );
    DrawHeightField(StateHeight);

    // Label.
    fill( 1.0 );
    text( "2D Wave Equation Init Only", 10, 30 );
}

Let's look at that in Processing, with just the drawing of the initial state, and nothing else:

Discretizing the 2D Wave Equation

In the 1D solver, we implemented our time-step by a somewhat lengthy discretization of the 1D wave equation, which we massaged and rearranged until it was suitable for use with our simple Jacobi iterative matrix solver. We'll follow exactly the same steps again. We'll take some bigger leaps this time around, but the method is exactly the same as in the 1D case, which you can review to understand each step broken down.

The 2D Wave Equation describes the motion of waves in a height field, and can be written simply as:

\[ \begin{equation} \frac{\partial^2 h}{\partial t^2} = c^2 \nabla^2 h \end{equation} \]

Which is the same as:

\[ \begin{equation} \frac{\partial^2 h}{\partial t^2} = c^2 \Delta h \end{equation} \]

Which are both just different shorthand notations for the following differential equation:

\[ \begin{equation} \frac{\partial^2 h}{\partial t^2} = c^2 (\frac{\partial^2 h}{\partial x^2} + \frac{\partial^2 h}{\partial y^2}) \end{equation} \]

The form of this equation simply adds the second partial derivative of height with respect to y to the right hand side of the 1D wave equation. The resulting discretization is simple and looks very familiar, we follow the same series of steps as we did with the 1D equation to arrive at:

\[ \begin{equation} \begin{split} a_{{i,j}_{t+\Delta t}} = &c^2 \frac{h^\star_{{i+1,j}_{t+\Delta t}} + \ h^\star_{{i-1,j}_{t+\Delta t}} + \ h^\star_{{i,j+1}_{t+\Delta t}} + \ h^\star_{{i,j-1}_{t+\Delta t}} - \ 4 h^\star_{{i,j}_{t+\Delta t}} } {\Delta x^2} + \\ &c^2 \Delta t^2 \frac{a_{{i+1,j}_{t+\Delta t}} + \ a_{{i-1,j}_{t+\Delta t}} + \ a_{{i,j+1}_{t+\Delta t}} + \ a_{{i,j-1}_{t+\Delta t}} - \ 4 a_{{i,j}_{t+\Delta t}} } {\Delta x^2} \label{eq:AccelDisc1Subscripted} \end{split} \end{equation} \]

The above equation assumes that \(\Delta x = \Delta y\), as does the derivation below. If working with non-square cells, the discretization gets a little longer, left as an exercise for the reader! (It rarely occurs in practice). Following the same process as in the 1D case, we rearrange our terms, substituting intermediate constants \(\kappa\) and \(\gamma\):

\[ \begin{equation} \kappa = \frac{c^2 \Delta t^2}{\Delta x^2}, \ \gamma = \frac{c^2}{\Delta x^2} \label{eq:Constants} \end{equation} \]

Giving the following equation which expresses the relationship between the acceleration at the index \((i,j)\) and its spatially adjacent neighbors, dropping the time subscript because it is unnecessary:

\[ \begin{equation} \begin{split} (1 + 4 \kappa) a_{i,j} + &( -\kappa ) ( a_{i-1,j} + a_{i+1,j} + a_{i,j-1} + a_{i,j+1} ) = \\ &\gamma ( h^\star_{i-1,j} + \ h^\star_{i+1,j} + \ h^\star_{i,j-1} + \ h^\star_{i,j+1} - \ 4 h^\star_{i,j} ) \\ \label{eq:AccelOneMatrixRow} \end{split} \end{equation} \]

The above equation \(\eqref{eq:AccelOneMatrixRow}\) represents a single row in the linear matrix we are solving. Once again, we'll use a jacobi solver, computing a new guess for acceleration \(a^{\star k+1}\) in each iteration \(k+1\) from the acceleration \(a^{\star k}\) of the previous iteration, which looks like this:

\[ \begin{eqnarray} b_{i,j} &=& \gamma ( h^\star_{i-1,j} + \ h^\star_{i+1,j} + \ h^\star_{i,j-1} + \ h^\star_{i,j+1} - \ 4 h^\star_{i,j} ) \\ c^{\star k}_{i,j} &=& \kappa ( a^{\star k}_{i-1,j} + \ a^{\star k}_{i+1,j} + \ a^{\star k}_{i,j-1} + \ a^{\star k}_{i,j+1} ) \\ a^{\star k+1}_{i,j} &=& \frac{b_{i,j} + c^{\star k}_{i,j}}{1 + 4 \kappa}\\ \label{eq:JacobiIteration} \end{eqnarray} \]

This translates very directly into code:

// Jacobi iteration to get temp acceleration
void JacobiIterationAccel( int i_aOld, int o_aNew, int i_hStar, float i_dt ) {
    float kappa = sq( WaveSpeed ) * sq( i_dt ) / sq( DXY );
    float gamma = sq( WaveSpeed ) / sq( DXY );

    for (int j = 1; j < NY-1; ++j) {
        for (int i = 1; i < NX-1; ++i) {
            float a_left = State[i_aOld][IX(i-1,j)];
            float a_right = State[i_aOld][IX(i+1,j)];
            float a_down = State[i_aOld][IX(i,j-1)];
            float a_up = State[i_aOld][IX(i,j+1)];

            float h_star_left = State[i_hStar][IX(i-1,j)];
            float h_star_right = State[i_hStar][IX(i+1,j)];
            float h_star_down = State[i_hStar][IX(i,j-1)];
            float h_star_up = State[i_hStar][IX(i,j+1)];
            float h_star_cen = State[i_hStar][IX(i,j)];

            float b = gamma *
                (h_star_left + h_star_right + h_star_down + h_star_up -
                 (4.0 * h_star_cen));

            float c = kappa * (a_left + a_right + a_down + a_up);

            State[o_aNew][IX(i,j)] = (b + c) / (1.0 + kappa);
        }
    }

    EnforceDirichletBoundaryConditions( o_aNew );
}

Arduously Refactoring What's Left

We've made these changes so far in our port from 1D to 2D:

• Updated the initial conditions to work in 2D
• Updated boundary condition enforcement to work in 2D
• Added 2D drawing code
• Updated the jacobi iteration to work in 2D

Now we need to undertake the work of porting the remainder of the solver, which is still the bulk of the code, by length. Maybe go get a coffee and have a stretch before this undertaking. Ready?

And... we're done. (We didn't have to do anything at all)

By organizing our state vector into generic arrays, and sticking with an underlying one-dimensional layout of the arrays themselves, most of the actual solver machinery - the estimation of new fields, the time integration - works without modification. Additionally, it's modular. We can swap out RK4 for RK2, or try a first order time step to see its effect on stability. This modularity and flexibility is the payoff for the good design practices we established at the outset, and it translates into the design and implementation of bigger systems, even (especially) production-quality solvers. It also tends to abstract out parallelism very effectively.

I've added a little extra feature to this final solver - hit the 't' key to cycle between first-order, RK2, and RK4 time steps. It's interesting to notice the increase in stability with each one, though the gap between first-order and RK2 is much bigger than the gap between RK2 and RK4.

There's a bit of aliasing around the mouse input point, and the waves created have a small amount of aliasing to them as well. As an exercise for the reader, consider these two approaches towards fixing these problems. Firstly, instead of just drawing a single one-pixel dot as the mouse input boundary condition, draw a smooth shape instead. This is challenging because the application of boundary conditions needs to be idempotent - meaning once they're applied, applying them a second time should have no effect. Secondly, we're using a first-order discretization of space in the discretization of our laplacian. If you used a higher-order spatial discretization, what would that look like? How would it change the Jacobi iteration?

Download

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