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

/*Currently, this simulation freaks out drawing weird green, yellow, and red patterns.

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;

}

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

Let me know what code you changed.

/*Again, let me know what code you changed.

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;

}

"derivative of u with respect to x" = "delta u delta x" = du / dx = ∂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;

}

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).

∂

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:

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

∂

If you calculate ∂

Turn in your modified code snippets on Blackboard.