Wave Simulation with CUDA

Two things that I've been wanting to do for a long time are:

1. Write my own physics simulation that runs on my laptop's graphics card, and
2. Mess around with wave visualizations.

The Concept

I wanted code that could help me think about spacial dynamics rather than just be another complicated thing I have to deal with. I also have a bit to learn before attempting really relevant problems anyway, so I decided to at least begin to build something that can function as a sort of general-purpose thinking, computing and visualization tool that can grow with me as I grow.

The idea behind the code goes something like this: imagine a point in space, just a single point without any spacial extent of its own. That point has certain properties that we might care about in a given situation: x, y, and z coordinates; mass; amplitude; whatever. The main thing is just to clearly define the properties associated with this point that are important for the problem and to ignore everything else.

At any given moment, those properties will have certain values, but at the next moment those values may have changed. So there's a function that goes along with our point that determines how its property values change from one moment to the next. Figuring out what that function is is called Physics, and yes, it's the hard part. I like to imagine myself as the point, and then think about what forces might be pulling on me, and in what directions. I guess the whole "point" of this computer program is to get away from all the extraneous complexity and clear a space where it's easy to think about just one thing.

The Code

I started by defining a "point" data type and giving it the properties it needs. Since I was doing a pretty simple 2D wave simulaton for this project all I needed were x, y, and Ψ.

typedef struct {

    float x;
    float y;
    float psi;

} point;

Now I could envision my "point" as a zero-dimensional "container" that carries all of its relevant data around with it wherever it goes, make as many copies of the point as I like, and scatter them around space by assigning unique values to the properties of each and storing them all in a giant array that gets passed to the GPU for number-crunching.

Back in the mindset of a single point, I wrote a GPU function (a.k.a. a "kernel") that takes the properties of the point at a single instant in time and replaces them with the values they will have at the next instant. A simplified version of it looks like this:

__device__ float t = 0;

__global__ void incrementPoint( point *p ) {

    // CUDA's index for this "point":
    int	i = blockIdx.x * blockDim.x
            + threadIdx.x;

    // A disbursive wave packet:
    float r = sqrt( p[i].x*p[i].x
              + p[i].y*p[i].y );
    float alpha = r - c * t;
    p[i].psi = 2 * A * dk * sin( dk * alpha )
          * cos( k * alpha ) / ( dk * alpha );

    // Incrementing the time:
    t += dt;
}

The Rest

was really just filling in the details, which of course in practice was less than trivial; learning to program in CUDA is not for the feint of heart! In fact I probably spent 99% of my time on these "details." But it was a worthwhile effort, because now I can plug any equation I like into incrementPoint() and get a totally different simulation (with only a few other minor tweaks.)

2D mesh interpolated over 32,768 points