# 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 CPUCreates 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.

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

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

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:

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:

Turn in your modified code snippets on Blackboard.