HW3: Gridded Simulations and PDEs

CS 493/693 Homework, Dr. Lawlor

Prepare a single plain-text file (.txt) giving me your results from this homework.

Problem 1: Basic Gridded Simulations

Here's a nice little 1D simulation of the shallow water wave equations: x is velocity, y is height.  It can be run on NetRun, or on your own machine using the osl/vec4.h present in every example program (view "out.ppm" using the GIMP).
/*
493/693 hw3.1: 1D shallow-water wave equation on the CPU

Creates an output PPM image with time increasing down the vertical axis.
*/
#include <iostream>
#include <fstream> /* for ofstream */
#include <vector> /* data storage */
#include "osl/vec4.h" /* data format */

int main(void)
{
// Figure out how big an image we should render:
int wid=400, ht=400;

// Create a PPM output image file header:
std::ofstream out("out.ppm",std::ios_base::binary);
out<<"P6\n"
<<wid<<" "<<ht<<"\n"
<<"255\n";

std::vector<vec4> state(wid,vec4(0.0)); /* read here */
std::vector<vec4> newstate(wid,vec4(0.0)); /* write here */
for (int i=wid/2-20;i<wid/2;i++)
{ /* prepare the initial conditions */
state[i-100]=vec4(0.4,0.4,0.5,0.0);
state[i]=vec4(0.0,-0.5,0.5,0.0);
state[i+100]=vec4(0.4,-0.4,0.5,0.0);
}

int substep=8; /* simulation steps per image output row */
for (int time=0;time<substep*ht;time++)
{
for (int x=0;x<wid;x++)
{ /* Run physics */
float dt=0.1;
/* Left and right boundary conditions */
if (x==0) newstate[x]=vec4(0.0,state[x+1].y,0.0,0.0);
if (x==wid-1) newstate[x]=vec4(0.0,state[x-1].y,0.0,0.0);
/* Interior */
if (x-1>=0 && x+1<wid) {
vec4 L=state[x-1]; /* left, center, and right */
vec4 C=state[x];
vec4 R=state[x+1];

/* Run physics */
C.y+=dt*(L.x-R.x);
C.x+=dt*(L.y-R.y);

newstate[x]=C;
}
}
if (time%substep==0) /* do output this step */
for (int x=0;x<wid;x++)
{ /* Figure out the output pixel color */
unsigned char r,g,b;
vec4 outcolor=newstate[x];
outcolor=clamp(outcolor+vec4(0.1),0.0,1.0);
r=(unsigned char)255*outcolor.x;
g=(unsigned char)255*outcolor.y;
b=(unsigned char)255*outcolor.z;
out<<r<<g<<b;
}

/* After each timestep, ping-pong the buffers */
std::swap(state,newstate);
}

return 0;
}

(Try this in NetRun now!)

Currently, this simulation freaks out drawing weird green, yellow, and red patterns. 
Initial simulation result.

Fix it (a small change) so that it draws a nice simulation that looks like this.
Smooth gridded simulation result.

Let me know what code you changed.

Problem 2: GPGPU Simulations

This GPGPU code has the same bug.  But luckily, the same fix should work!
/*
493/693 hw3.1: 1D shallow-water wave equation on the CPU

Creates an output PPM image with time increasing down the vertical axis.
*/
#include <iostream>
#include <fstream> /* for ofstream */
#include <vector> /* data storage */
#include "ogl/gpgpu.h" /* all GPU stuff below */
#include "ogl/glsl.cpp" /* to simplify linking */
/* will need GLEW at runtime: either #include "ogl/glew.c", or link in yourself. */

int main(void)
{
// Figure out how big an image we should render:
int wid=400, ht=400;

// Create a PPM output image file header:
std::ofstream out("out.ppm",std::ios_base::binary);
out<<"P6\n"
<<wid<<" "<<ht<<"\n"
<<"255\n";

std::vector<vec4> state(wid,vec4(0.0)); /* read here */
std::vector<vec4> newstate(wid,vec4(0.0)); /* write here */
for (int i=wid/2-20;i<wid/2;i++)
{ /* prepare the initial conditions */
state[i-100]=vec4(0.4,0.4,0.5,0.0);
state[i]=vec4(0.0,-0.5,0.5,0.0);
state[i+100]=vec4(0.4,-0.4,0.5,0.0);
}

/* GPGPU stuff: */
gpu_env e;
gpu_array State(e,"state",wid,1,state[0]);
gpu_array Newstate(e,"newstate",wid,1,newstate[0]);

int substep=8; /* simulation steps per image output row */
for (int time=0;time<substep*ht;time++)
{
GPU_RUN(Newstate,
/* Run physics */
float dt=0.1;

vec2 delx=vec2(stateScale.x,0.0); /* one pixel */
vec4 L=texture2D(state,location-delx);
vec4 C=texture2D(state,location);
vec4 R=texture2D(state,location+delx);

/* Left and right boundary conditions */
if (location.x<stateScale.x) /* left edge */
gl_FragColor=vec4(0.0,R.y,0.0,0.0);
else if (location.x>1.0-stateScale.x) /* right edge */
gl_FragColor=vec4(0.0,L.y,0.0,0.0);
else /* Interior */
{
/* Run physics */
C.y+=dt*(L.x-R.x);
C.x+=dt*(L.y-R.y);

gl_FragColor=C;
}
)
if (time%substep==0) /* do output this step */
{
Newstate.read(newstate[0],wid,1);
for (int x=0;x<wid;x++)
{ /* Figure out the output pixel color */
unsigned char r,g,b;
vec4 outcolor=newstate[x];
outcolor=clamp(outcolor+vec4(0.1),0.0,1.0);
r=(unsigned char)255*outcolor.x;
g=(unsigned char)255*outcolor.y;
b=(unsigned char)255*outcolor.z;
out<<r<<g<<b;
}
}
/* After each timestep, ping-pong the buffers */
State.swapwith(Newstate);
}

return 0;
}

(Try this in NetRun now!)

Again, let me know what code you changed.

Problem 3: Partial Differential Equations

First, a note on the partial derivative syntax I'm using below.  There are a bunch of equivalent notations to write the derivative of some value 'u' with respect to x:
    "derivative of u with respect to x" = "delta u delta x" = du / dx = ∂u / ∂x = ∂x u = ux = u_x

Only that last version is a valid C++ (or GLSL) variable name.  Thus a velocity (a derivative with respect to time) can be written as "u_t" in that last notation.  A second derivative, like acceleration (the second derivative with respect to time), can be written as "u_tt" in that last notation. 

In the code below, I'm using a centered difference approximation to calculate the second spatial derivative of our main value 'u' (the red channel): that is, u_xx.  I'm using the usual leapfrog integration we used for particles to integrate the second time derivative u_tt.  Evidently the 1D wave equation is simply "u_tt = u_xx", which seems to generate a moving diamond pattern that bounces around somewhat like waves (see Stephen Wolfram, A New Kind of Science, page 163).  Here's the code:
/*
493/693 hw3.3: 1D PDE-style wave equation on the CPU

Creates an output PPM image with time increasing down the vertical axis.
*/
#include <iostream>
#include <fstream> /* for ofstream */
#include <vector> /* data storage */
#include "osl/vec4.h" /* data format */

int main(void)
{
// Figure out how big an image we should render:
int wid=400, ht=400;

// Create a PPM output image file header:
std::ofstream out("out.ppm",std::ios_base::binary);
out<<"P6\n"
<<wid<<" "<<ht<<"\n"
<<"255\n";

std::vector<vec4> state(wid,vec4(0.0)); /* read here */
std::vector<vec4> newstate(wid,vec4(0.0)); /* write here */
for (int i=0;i<wid;i++)
{ /* prepare the initial conditions */
float x=(i-wid/2.0)/wid;
float gaussian=exp(-x*x*10000);
state[i]=vec4(1.0*gaussian,0.0,0.5,0.0);
}

int substep=8; /* simulation steps per image output row */
for (int time=0;time<substep*ht;time++)
{
for (int x=0;x<wid;x++)
{ /* Run physics */
float dt=0.1;
/* Interior */
if (x-1>=0 && x+1<wid) {
float L=state[x-1].x; /* left, center, and right */
float u=state[x].x;
float R=state[x+1].x;
float u_t=state[x].y; /* velocity (at center) */

/* u_xx: second derivative in space (centered difference) */
float u_xx=L-2*u+R;

/* u_tt: second derivative in time (acceleration)
Setting this to u_xx makes the wave equation.
*/
float u_tt=u_xx;

/* u_t: first derivative in time (velocity) */
u_t += dt*u_tt;

/* Leapfrog update for u */
u += dt*u_t;

newstate[x]=vec4(u,u_t,u_tt,0.0);
}
}
if (time%substep==0) /* do output this step */
for (int x=0;x<wid;x++)
{ /* Figure out the output pixel color */
unsigned char r,g,b;
vec4 outcolor=newstate[x];
outcolor=clamp(outcolor+vec4(0.1),0.0,1.0);
r=(unsigned char)255*outcolor.x;
g=(unsigned char)255*outcolor.y;
b=(unsigned char)255*outcolor.z;
out<<r<<g<<b;
}

/* After each timestep, ping-pong the buffers */
std::swap(state,newstate);
}

return 0;
}

(Try this in NetRun now!)

This initially produces a diamond pattern as the initial disturbance propagates sideways in both directions, until it bounces off the walls.  There's a little tiny bit of dispersion (ringing).
Diamond pattern from 1D wave equation

I'd like to you modify this simulation to solve the first PDE listed in Wolfram, page 165:

    ∂tt u[t,x] = ∂xx u[t,x] + (1 - u[t,x]2)(1 + u[t,x])

Keep in mind that his "u[t,x]" is just plain "u" above.  The whole point is this is easier than it sounds.

Your modified simulation should look like this:
Wolfram's page 165 PDE solution.

This actually isn't unstable, it's just a colorful pattern.  Note the weirdly organic interacting patterns!

Problem 4: More complex PDEs (for 693 students only)

You can use the same code as problem 3 to simulate this PDE, which is a nonlinear modification of the wave equation:

tt u = (∂xx u) / (0.5 + ∂x u)

If you calculate ∂x u using a left or right handed difference, the simulation breaks down and produces an ugly yellow blotch, but a centered difference approximation produces this nice skewed diamond:

Skewed wave equation

Turn in your modified code snippets on Blackboard.