Monday, June 2, 2014

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

Monday, May 26, 2014

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

Monday, June 3, 2013

Simulation Class: Iterative Matrix Solvers and the Jacobi Method

Iterative Methods for Large Linear Systems

In our previous post, we carefully discretized the equation of motion for a 1D wave equation into a large set of simple linear equations, which we then converted into a sparse matrix notation. It's important to understand that the matrix & vector notation that we reached at the end of our derivation, after creating a list of simultaneous linear equations, was just a change of notation, nothing more. Matrices and the notation surrounding linear algebra can sometimes be obfuscating and generate confusion - whenever that is the case, it's helpful to just remember that the matrix notation is a short-hand, and that we can always revert to the system of equations that the matrix expression represents.

So, armed with this fantastically simple linear matrix equation, we are tasked with solving the system for the unknown solution vector, \(x\):

\[ \begin{equation} A x = b \label{eq:AxEqualsB} \end{equation} \]
Surely we can just construct our matrix \(A\), invert the matrix to produce \(A^{-1}\), and then our solution would be:
\[ \begin{equation} x = A^{-1} b \label{eq:XequalsAinvB} \end{equation} \]
This simple solution is comforting, but sadly, totally unusable in any practical setting. Why?

The first and most important problem we run into with an explicit matrix inversion solution is that the matrices themselves would be incredibly massive in size for any real-world problem. Imagine we have a fluid simulation in 2D that is \(1024 \times 1024\) in size. In any such simulation, our state array size \(N\) is equal to the total number of points being solved for, so in this case, \(N = 1024 \times 1024 = 1048576\). The matrix \(A\) is, as we've seen, size \(N \times N\), which would in that simple case be \(N^4 = 109951162776\), or just slightly larger than \(10^{12}\) points, which would take just over 4 TERABYTES of RAM to hold, for single precision floating point.

Sparse Matrices

The obvious immediate solution to the above storage problem is to take advantage of the fact that our matrices are sparse. A Sparse Matrix is one in which most of the matrix elements are zero. In our example, we notice that there is at most three non-zero elements per-row, one of which is always the diagonal element. The other non-zero elements are always one column before or after our diagonal element. We could be clever and store our matrix \(A\) as an array of \(N\) "Sparse Rows", where a Sparse Row has a structure something like this:


class SparseMatrixRow
{
    float indices[3];
    float values[3];
    int numIndices;
    SparseMatrixRow()
    {
        numIndices = 0;
    }
    void addElement( int i_index, float i_val )
    {
        if ( numIndices < 3 )
        {
            values[numIndices] = i_val;
            indices[numIndices] = i_index;
            ++numIndices;
        }
    }
};

Indeed, structures similar to the one above are the basis of sparse-matrix linear algebra libraries such as Sparse BLAS. With the above representation, our matrix in the 2D height field case would take about 28 megabytes to store, a huge improvement over 4 terabytes.

However, even with the sparse representation of \(A\), how do we invert a huge matrix like this? If you've ever worked to invert a \(3 \times 3\) or \(4 \times 4\) matrix, you know how complex that could be - how on earth would you go about inverting a square matrix with a million rows, and how would you store your intermediate results, which may not enjoy the simple sparse form that the original matrix \(A\) had?

Problems like the above are what supercomputers were essentially built to solve, and there are excellent libraries and tools for doing this. However, such systems are arcane in nature, usually written in FORTRAN, and in our case, drastically more firepower than we actually need.

Iterative Methods

Our matrix has some nice properties that will allow us to solve the problem much more simply. Firstly, it is symmetric. A Symmetric Matrix is one which is identical under a transposition of the rows and columns. If you take any element to the upper right of the diagonal, and look at the mirror of that element in the same location on the lower left of the diagonal, you'll see the same value. A symmetric matrix as a representation of a physical equation of motion is an expression of Newton's Third Law, that each action has an equal but opposite reaction. It means that the effect point \(i\) has on point \(j\) is the same as the effect point \(j\) has on point \(i\). Conservation of energy emerges from this law.

Secondly, assuming that our timestep \(\Delta t\) is relatively small compared to the ratio of wave propagation speed \(c\) to the spatial discretization \(\Delta x\), our off-diagonal elements will be much smaller in absolute magnitude than our diagonal elements. This means that our system is a candidate for The Jacobi Method, which is the simplest iterative \(A x = b\) solver.

Iterative methods, which include Jacobi, but also include more fancy solvers such as The Conjugate Gradient Method, all involve some variation on making an initial guess for the unknown solution vector \(x^{\star0}\), then using that guess in some way to produce a better guess \(x^{\star1}\), and so on, until the change between guesses from one iteration to the next falls below some acceptance threshold, or a maximum number of iterations is reached.

For all iterative methods, we will require storage for the next guess and the previous guess, and we'll swap these temporary storage arrays back and forth as we iterate towards the solution. Our generic state-vector implementation easily allows for this, and indeed we already used this facility when we computed our RK4 approximation.

Jacobi Method

The Jacobi Method is so simple, we don't even need to construct a sparse matrix or a solution vector, we can simply solve for the new estimate of the solution vector \(x\) in place.

Jacobi works by taking each equation and assuming that everything other than the diagonal element is known - which we do by using values from the previous guess of \(x\). The new guess for the diagonal element \(x\) is then, at any row of the matrix:

\[ \begin{equation} x^{\star k+1}_i = \frac{b_i - \sum_{j\neq i}A_{i,j} x^{\star k}_{j}}{A_{i,i}} \label{eq:JacobiIteration} \end{equation} \]
In the above Equation \(\eqref{eq:JacobiIteration}\), the new guess iteration \(k+1\) is computed from the values of the old guess iteration \(k\). This is fairly easy to code, and looks like this:


// Jacobi iteration to get temp acceleration
void JacobiIterationAccel( int i_aOld, int o_aNew, int i_hStar, float i_dt )
{
    float aij = -sq( WaveSpeed ) * sq( i_dt ) / sq( DX );
    float aii = 1.0 + ( 2.0 * sq( WaveSpeed ) * sq( i_dt ) / sq( DX ) );
    float bgain = sq( WaveSpeed ) / sq( DX );
    
    for ( int i = 1; i < ArraySize-1; ++i )
    {
        float aLeft = State[i_aOld][i-1];
        float aRight = State[i_aOld][i+1];
        
        float hStarLeft = State[i_hStar][i-1];
        float hStarCen = State[i_hStar][i];
        float hStarRight = State[i_hStar][i+1];
        
        float bi = bgain * ( hStarLeft + hStarRight - ( 2.0 * hStarCen ) );
        
        float sumOffDiag = ( aij * aLeft ) + ( aij * aRight );
        
        State[o_aNew][i] = ( bi - sumOffDiag ) / aii;
    }
    
    EnforceAccelBoundaryConditions( o_aNew );
}

The code above computes a single new guess for the accelerations from an old guess for the accelerations. We can then write our "acceleration from position" function as we had before, as follows. We're taking a fixed number of iterations for guessing:


// Solve for acceleration.
void JacobiSolveAccel( int i_hStar, float i_dt )
{  
    // Initialize acceleration to zero.
    FillArray( StateAccelStar, 0.0 );
    
    // Solve from StateJacobiTmp into StateAccel
    for ( int iter = 0; iter < 20; ++iter )
    {
        int tmp = StateAccelStar;
        StateAccelStar = StateJacobiTmp;
        StateJacobiTmp = tmp;
        
        JacobiIterationAccel( StateJacobiTmp, StateAccelStar, i_hStar, i_dt );
    }
}

void EstimateAccelStar( float i_dt )
{
    JacobiSolveAccel( StateHeightStar, i_dt );
}

The rest of our implementation looks basically identical. Here's the final timestep, which looks exactly like our previous RK4 implementation:


// Time Step function.
void TimeStep( float i_dt )
{
    // Swap state
    SwapState();
    
    // Initialize estimate. This just amounts to copying
    // The previous values into the current values.
    CopyArray( StateHeightPrev, StateHeight );
    CopyArray( StateVelPrev, StateVel );
    
    // 1
    CopyArray( StateVel, StateVelStar );
    EstimateHeightStar( i_dt );
    EstimateAccelStar( i_dt );
    AccumulateEstimate( i_dt / 6.0 );
    
    // 2
    EstimateVelStar( i_dt / 2.0 );
    EstimateHeightStar( i_dt / 2.0 );
    EstimateAccelStar( i_dt / 2.0 );
    AccumulateEstimate( i_dt / 3.0 );
    
    // 3
    EstimateVelStar( i_dt / 2.0 );
    EstimateHeightStar( i_dt / 2.0 );
    EstimateAccelStar( i_dt / 2.0 );
    AccumulateEstimate( i_dt / 3.0 );
    
    // 4
    EstimateVelStar( i_dt );
    EstimateHeightStar( i_dt );
    EstimateAccelStar( i_dt );
    AccumulateEstimate( i_dt / 6.0 );
    
    // Final boundary conditions on height and vel
    EnforceHeightBoundaryConditions( StateHeight );
    EnforceVelBoundaryConditions( StateVel );

    // Update current time.
    StateCurrentTime += i_dt;
}

And that's it... Let's see how that looks in Processing!

And with that, we're finished... (for now!)

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: Matrix Methods and the Wave Equation

The Wave Equation as a Linear System

In our previous class, we completed our exploration of higher-order explicit time-integration techniques, and moved from a simple spring system to a more complex 1-D Wave Equation. Using the fourth-order Runge-Kutta (RK4) integration scheme, we adapted our spring solver to the wave equation, which seemed stable at first, but as soon as we introduced external forces, the system showed immediate instabilities.

What follows in this post is a step-by-step progression, starting with the original PDE of the 1D Wave Equation, and eventually producing a linear system, represented by the equation \(A x = b\), where \(A\) is a square, mostly sparse matrix of size \(N\times N\), where \(N\) is the number of unknown values in the discretized simulation grid (either heights, or accelerations). This post is a lot of math, but all of it is manipulation of simple terms and equations in a linear algebra setting. I've broken things down to the smallest steps because this is such a core concept to how fluid simulations and grid simulations are performed.

The problem in our previous solution is related to how we calculated our fluid vertical acceleration from the fluid's height. To start with, the 1D Wave Equation is written as:

\[ \begin{equation} \frac{\partial^2 h}{\partial t^2} = c^2 \frac{\partial^2 h}{\partial x^2} \label{eq:Wave1D} \end{equation} \]

Where \(h\) represents height at a given location along an axis, and \(c^2\) represents the square of the wave propagation speed. I think this is my favorite equation, because of the implied equivalency of spatial and temporal dimensions. It says that, as far as waves are concerned, space and time are interchangeable (assuming you choose your units correctly). Anyway, let's examine how we discretized this solution for a given point in our height array, \(h_i\) and its acceleration \(a_i\).

\[ \begin{equation} a_i = c^2 \frac{h_{i+1} -\ 2 h_i + \ h_{i-1}}{\Delta x^2} \label{eq:WaveDisc1} \end{equation} \]
The above equation, \(\eqref{eq:WaveDisc1}\) is just a discretization of Equation \(\eqref{eq:Wave1D}\) above, with \(\frac{\partial^2 h}{\partial t^2}\) taken as acceleration. However, unlike the simple spring system, where the position of the spring varied only in time, our position varies in space as well, which means we have to consider that every point in the array of heights depends on not only its previous positions and velocities, but the heights and velocities and accelerations of neighboring points. This was the oversimplification in our previous attempt, and the source of the instability.

In order to understand this system better, let's replace the acceleration with a finite difference approximation, based on the height and previous heights. Starting with the backwards finite difference approximation of \(a_i\):

\[ \begin{equation} a_{i_{t+\Delta t}} \approx \frac{h_{i_{t+\Delta t}} - \ 2 h_{i_{t}} + \ h_{i_{t-\Delta t}}}{\Delta t^2} \label{eq:FiniteDiffA} \end{equation} \]
we can follow the usual convention in height-field solvers, eliminating acceleration by substituting Equation \(\eqref{eq:FiniteDiffA}\) into Equation \(\eqref{eq:WaveDisc1}\) to get the following expression which contains only height values (at different times and positions):
\[ \begin{equation} \frac{h_{i_{t+\Delta t}} - \ 2 h_{i_{t}} + \ h_{i_{t-\Delta t}}}{\Delta t^2} = c^2 \frac{h_{{i+1}_{t+\Delta t}} - \ 2 h_{i_{t+\Delta t}} + \ h_{{i-1}_{t+\Delta t}} } {\Delta x^2} \label{eq:WaveHeightDisc1} \end{equation} \]
However, our RK4 solver and our other time-integration techniques were written to employ a calculation of acceleration based on an estimate of current position and velocity, so it's more helpful if we instead begin with the following implicit relationship between acceleration, velocity, and height:
\[ \begin{equation} h_{i_{t+\Delta t}} = h_{i_t} + \Delta t\ v_{i_{t+\Delta t}} + \Delta t^2\ a_{i_{t+\Delta t}} \label{eq:ImplicitHeightTimeStep} \end{equation} \]
Let's treat the first two terms of the right hand side of Equation \(\eqref{eq:ImplicitHeightTimeStep}\) as an estimate of height at time \(t+\Delta t\), denoted by \(h^\star\), indicated as follows:
\[ \begin{equation} h^\star_{i_{t+\Delta t}} = h_{i_t} + \Delta t\ v_{i_{t+\Delta t}} \label{eq:HeightEstimate} \end{equation} \]
And then we can substitute Equation \(\eqref{eq:HeightEstimate}\) into Equation \(\eqref{eq:ImplicitHeightTimeStep}\) to get:
\[ \begin{equation} h_{i_{t+\Delta t}} = h^\star_{i_{t+\Delta t}} + \Delta t^2\ a_{i_{t+\Delta t}} \label{eq:EstimatedHeightTimeStep} \end{equation} \]
Now, we can subsitute Equation \(\eqref{eq:EstimatedHeightTimeStep}\) into Equation \(\eqref{eq:WaveDisc1}\) to get the following equation, written in terms of the current acceleration and an estimate of the current height.
\[ \begin{equation} a_{i_{t+\Delta t}} = c^2 \frac{h^\star_{{i+1}_{t+\Delta t}} - \ 2 h^\star_{i_{t+\Delta t}} + \ h^\star_{{i-1}_{t+\Delta t}} } {\Delta x^2} + c^2 \Delta t^2 \frac{a_{{i+1}_{t+\Delta t}} - \ 2 a_{i_{t+\Delta t}} + \ a_{{i-1}_{t+\Delta t}} } {\Delta x^2} \label{eq:AccelDisc1Subscripted} \end{equation} \]
All of the terms in Equation \(\eqref{eq:AccelDisc1Subscripted}\) are at the same point in time, so we can drop the time subscript. Assuming that the estimates of height are given by some process (which we are free to tinker with), the equation is a linear relationship between accelerations at different points in space. Let's simplify by defining intermediate constants:
\[ \begin{equation} \kappa = \frac{c^2 \Delta t^2}{\Delta x^2}, \ \gamma = \frac{c^2}{\Delta x^2} \label{eq:Constants} \end{equation} \]
Substituting in our constants, moving all of the acceleration terms to the left-hand side, and gathering coefficients, we have the following equation which expresses the relationship between the acceleration at the index \(i\) and its spatially adjacent neighbors:
\[ \begin{equation} (1 + 2 \kappa) a_i + ( -\kappa ) a_{i-1} + ( -\kappa ) a_{i+1} = \gamma ( h^\star_{i+1} - \ 2 h^\star_i + \ h^\star_{i-1} ) \label{eq:AccelOneMatrixRow} \end{equation} \]
The above Equation \(\eqref{eq:AccelOneMatrixRow}\) is written relative to any position in the array of heights, denoted by the subscript variable \(i\), and this relationship exists at each point in the array. We'll next write out the equations at each index explicitly, to get a system of linear equations. We have to take a bit of care at the boundary points. Our boundary condition is that heights at the boundary are equal to the height at the adjacent, non-boundary position, which, when transformed to a boundary condition on acceleration, says that the acceleration at the boundary is zero. Therefore, we have no equation for indices \(i = 0\) and \(i = N-1\), because we've explicitly defined those points. Furthermore, for indices \(i = 1\) and \(i = N-2\), the relationship is slightly changed, which we'll incorporate into the equations.
\[ \begin{eqnarray*} \ &\ & (1+2\kappa)a_1& +& (-\kappa)a_2& =& \gamma ( - h^\star_1 + h^\star_2 ) \\ (-\kappa)a_1& +& (1+2\kappa)a_2& +& (-\kappa)a_3& =& \gamma ( h^\star_1 - 2 h^\star_2 + h^\star_3 ) \\ (-\kappa)a_2& +& (1+2\kappa)a_3& +& (-\kappa)a_4& =& \gamma ( h^\star_2 - 2 h^\star_3 + h^\star_4 ) \\ ... \\ (-\kappa)a_{n-5}& +& (1+2\kappa)a_{n-4}& +& (-\kappa)a_{n-3}& =& \gamma ( h^\star_{n-5} - 2 h^\star_{n-4} + h^\star_{n-3} ) \\ (-\kappa)a_{n-4}& +& (1+2\kappa)a_{n-3}& +& (-\kappa)a_{n-2}& =& \gamma ( h^\star_{n-4} - 2 h^\star_{n-3} + h^\star_{n-2} ) \\ (-\kappa)a_{n-3}& +& (1+2\kappa)a_{n-2}&\ &\ & =& \gamma ( h^\star_{n-3} - h^\star_{n-2} ) \\ \label{eq:AccelLinearSystem} \end{eqnarray*} \]
This above system of linear equations can be expressed as a sparse matrix equation:
\[ \begin{equation} A x = b \end{equation} \]
Where the sparse, symmetric matrix \(A\) is defined as:
\[ \begin{equation} A = \begin{bmatrix} 1+2\kappa & -\kappa & \ & \ & \ & \ & \cdots \\ -\kappa & 1+2\kappa & -\kappa & \ & \ & \ & \cdots \\ \ & -\kappa & 1+2\kappa & -\kappa & \ & \ & \cdots \\ \ \vdots & \ & \ddots & \ddots & \ddots & \ & \cdots \\ \ \cdots & \ & \ & -\kappa & 1+2\kappa & -\kappa & \ \\ \ \cdots & \ & \ & \ & -\kappa & 1+2\kappa & -\kappa \\ \ \cdots & \ & \ & \ & \ & -\kappa & 1+2\kappa \\ \end{bmatrix} \end{equation} \]
and the column vectors \(x\) and \(b\) are defined as:
\[ \begin{equation} x = \begin{bmatrix} a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-4} \\ a_{n-3} \\ a_{n-2} \\ \end{bmatrix} , b = \begin{bmatrix} \gamma ( - h^\star_1 + h^\star_2 ) \\ \gamma ( h^\star_1 - 2 h^\star_2 + h^\star_3 ) \\ \gamma ( h^\star_2 - 2 h^\star_3 + h^\star_4 ) \\ \vdots \\ \gamma ( h^\star_{n-5} - 2 h^\star_{n-4} + h^\star_{n-3} ) \\ \gamma ( h^\star_{n-4} - 2 h^\star_{n-3} + h^\star_{n-2} ) \\ \gamma ( h^\star_{n-3} - h^\star_{n-2} ) \\ \end{bmatrix} \end{equation} \]
We have thus transformed our acceleration computation problem into an \(A x = b\) matrix problem, a problem which has been carefully studied, and for which many extremely well-understood methods of solution exist. We'll explore simple approaches to solving this system in the next post.

Download

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