CUDA 4.1 Particle Update

Question

I'm using CUDA 4.1 to parse in the update of my Particle system that I've made with DirectX 10. So far, my update method for the particle systems is 1 line of code within a for loop that makes each particle fall down the y axis to simulate a waterfall:

m_particleList[i].positionY = m_particleList[i].positionY - (m_particleList[i].velocity * frameTime * 0.001f);

In my .cu class I've created a struct which I copied from my particle class and is as follows:

struct ParticleType
{
    float positionX, positionY, positionZ;
    float red, green, blue;
    float velocity;
    bool active;
};

Then I have an UpdateParticle method in the .cu as well. This encompass the 3 main parameters my particles need to update themselves based off the initial line of code. :

   __global__ void UpdateParticle(float* position, float* velocity, float frameTime)
{

}

This is my first CUDA program and I'm at a loss to what to do next. I've tried to simply put the particleList line in the UpdateParticle method, but then the particles don't fall down as they should.

I believe it is because I am not calling something that I need to in the class where the particle fall code use to be. Could someone please tell me what it is I am missing to get it working as it should?

If I am doing this completely wrong in general, the please inform me as well.

Answer 1

You're thinking is slightly off, but you're on the right track. When you use CUDA, and update your particles on the GPU, depending on how many particles you have, you will roughly be dedicating one processor on the GPU to each particle. So your line of code in UpdateParticle will look very similar to your CPU call listed at the top of your question. The key is how to get the data to GPU, and how to access it in the device function.

I'll do my best to give an overview of all the parts you'll need, but for a more detailed reference, I recommend this book. If you happen to be a student somewhere, or have an account with Safari Books online, then you could start reading online for free.

So, the first thing you need to do, after you've created your particles on the CPU (host), is allocate memory for them on the GPU (device), and copy over the data. CUDA makes this very easy:

ParticleType* dev_particles; // The array of particles on the GPU.
// Assume that 'particles' is the pointer to your array of CPU particles
cudaMalloc((void**)&dev_particles, num_particles*sizeof(ParticleType));
cudaMemcpy(dev_particles, particles, num_particles*sizeof(ParticleType), cudaMemcpyHostToDevice);

So dev_particles is your pointer to memory on the GPU. It's important that you DO NOT try to access data from this via host code (code that runs on the CPU). Now you can write UpdateParticle as something like this:

__global__ void UpdateParticle(ParticleType* particles) {
    int index = blockIdx.x;
    particles[index].positionY = // Update the position here;
}

blockIdx is a CUDA built-in variable that stores the information about which block of the device is currently executing. What's happening, conceptually, is you will be launching N blocks on the GPU, where N is your number of particles, and each block will update only one particle. The 'x' value of the block determines which particle in the list that block will update. Note that you could mix this up however you want, so long as every particle is addressed, but this assignment makes the most sense.

Finally, to launch the GPU simulation of your particles, you would make a call like this:

UpdateParticles<<<num_particles,1>>>(dev_particles);

This launches num_particles blocks, and one thread per block. More information about threads can be found in that book I mentioned. Just note that the number of blocks you can launch is limited, and using threads can significantly increase the amount of particles you can operate on. For starters, I suggest just getting familiar with blocks.

If you want to operate on the particles in host code, after you've updated them on the device, you'll need to copy the data back over:

cudaMemcpy(particles, dev_particles, num_particles*sizeof(ParticleType), cudaMemcpyDeviceToHost);

and then particles will contain the updated data. Note that you can set up graphics interoperability with DirectX, allowing you not to have to copy data back and forth from the GPU each frame. However, that's beyond the scope of your question.

Also, if you need to send additional information to the device function, like frameTime, you might need to allocate space for that as well on the GPU, using a similar process as I used with the particles. It might be possible to just pass a single variable by value, but a pointer will not work, as a device cannot access host memory, and a host cannot access device memory.

If you have any questions, feel free to ask in comments.

Hope this helps :)

Answer 2

Wow, this ended up being a pretty long answer but I've explained it as thoroughly as I could, so it might be worth the read (or at least I hope so :D) I explain as much as I can about the theory, the practice and little optimizations you can make. Enjoy!

The theory

A kernel in CUDA is actually what gives you extreme parallel processing. CUDA follows a SPMD (Single Program-Multiple Data) approach: a single program (your kernel) is executed along several threads, manipulating multiple data in parallel. Bear in mind it's SPMD and not SIMD since your program code can branch (e.g. in if clauses, which you should avoid, since branching is heavy for GPUs) so it's not just a single instruction executed in multiple data (which would be a vector-processor, which GPUs are not! I hope I made this clear.)

Inside your kernel you should update your particles. Each thread executing your kernel would correspond to a particle (or a set of particles if you wish, which may help coalescing memory accesses.) That way you should IMMEDIATELY get rid of the for loop, wherever it's right now. Your kernel is practically your new loop (although it's not a loop per se.)

Since your particles don't affect each other you can safely update themselves inside the kernel, not needing a double buffer (which you would need to avoid race conditions.) This means you can safely overwrite each particle data because it doesn't matter to other particles' simulation.

I can see you barely got it, but it's weird how, although you do have a struct ParticleType, you're passing pointers to position and velocity to your kernel. Why are you not passing the struct array in the first place?

I'm going to continue this answer without using the struct, but if your data is laid as an array of ParticleType in GPU memory you should pass this array to the kernel, because if you didn't, you'd have to separate the data in different arrays (which would defeat all purposes of paralellization, since it would take a LONG time of memory accesses for each frame.)

Some more theory on kernels

Kernels are executed in batches of blocks (actually in batches of warps, which are 32 threads, but you don't need to know this right now since it's hardware-stuff.) These blocks contain up to 768 threads (in G80... the limit is higher in newer architectures) which can communicate with each other, either by shared memory (with the __shared__ qualifier) or by syncing (calling __synchthreads();.) Threads inside blocks are guaranteed to be able to do these tasks, but outside blocks, threads are oblivious to other threads (although not to the non-shared data they change.) Even blocks are oblivious to each other.

Once your kernel is executing you need to address your different data indices. But you're not passing any index as an argument! To address your particles, you have some convenient built-in variables:

1. blockIdx - Your current global block id.
2. blockDim - The global block dimension (threads per block.) Thus blockIdx.x * blockDim.x is the index of the first thread in your block.
3. threadIdx - The current thread id (inside the block.) blockIdx.x * blockDim.x + threadIdx.x would be your current thread index (easy if you think about it, firstThreadInBlockIndex + currentThreadInsideBlockIndex.)

This way you can address different indices using the following statement:

int idx = blockIdx.x * blockDim.x + threadIdx.x;

See the .x after each variable? These built-ins are actually dim3 (can have up to three dimensions) but you don't need this, since you're working on a 1d array so you just need the x dimension (just so you knew it.)

Allocating data in GPU

You first need to allocate and feed your buffers in GPU. If you're not modifying your data in your host code (which you shouldn't) you should just allocate the data once and then re-use the buffer in subsequent calls.

To allocate (and copy from host) a single buffer (if you don't use the ParticleType, you'll need one allocation for the position, one for the velocity buffer, etc.) you have to use these CUDA calls:

float* h_position; // Your host pointer. This holds the data (I assume it's already filled with the data.)
float* d_position; // Your device pointer, we will allocate and fill this
cudaMalloc( (void**)&d_position, numParticles * sizeof(float) );
cudaMemcpy( d_position, h_position, numParticles * sizeof(float), cudaMemcpyHostToDevice);

Pretty easy to understand.

First approach: naive implementation

Summing up, you should call your UpdateParticle kernel like this:

int threads_per_block = 128; // You should play with this value
int blocks = numParticles/threads_per_block + ( (numParticles%threads_per_block)1:0 );

I'm making a stop here to clarify the modulo (%) thing. You need it because if numParticles%threads_per_block is not zero you need one extra block for the remaining threads. E.g. 129 particles/128 part_per_blk = 1 blk, while you actually need two blocks (one for the first 128 threads, and one block for the remaining one.)

There is a problem with this approach. Think of the aforementioned example: since you're allocating 2 blocks with 128 threads per block, your actually spawning 256 threads... but your array has 129 elements! Not only some threads are going to repeat work but they're also going to do the work twice (which would lead to some particles moving twice as fast) or even worse (this is the case here) unintendedly writing to illegal indices.

You can safely overcome that problem checking for the correct indices. This is your complete kernel call:

UpdateParticle <<< blocks, threads_per_block >>> ( d_position, d_velocity, frameTime, numParticles );

Your kernel should look like this:

__global__ void UpdateParticle(float* position, float* velocity, float frameTime, int numParticles)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x; // Compute the index

    if (idx < numParticles) { // Is this index valid?
        position[idx] = position[idx] - velocity[idx] * frameTime * 0.001f;
        ... // some more updates
    }
}

Second approach: re-tessellating

As I said earlier, it might be wise to tessellate further inside your kernel. Your GPU bus is 128 bits, so reading 4 floats (4*32 bits = 128 bits) per thread is just as cheap as reading just one float. You might have to adjust this further if you end up using the ParticleType array. Just to try keep your data 128-bit aligned.

You have to pass a new argument to your kernel:

UpdateParticle <<< blocks, threads_per_block >>> ( d_position, d_velocity, frameTime, numParticles, particlesPerThread );

Bear in mind this would change your number of blocks, so change it accordingly.

Your kernel would almost be the same.

__global__ void UpdateParticle(float* position, float* velocity, float frameTime, int numParticles, int particlesPerThread)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if (idx < numParticles/particlesPerThread) {
        for ( int i = 0; i < particlesPerThread; i++ ) {
            position[idx+i] = position[idx+i] - velocity[idx+i] * frameTime * 0.001f;
           ... // some more updates
       }
    }
}

Also bear in mind your particlesPerThread should be a divisor of your total number of particles (or the for loop is going to access indices outside your array size in the very last thread.)

It might also be wise to unroll the loop (as I said earliear, braching is no good to GPUs) thus you wouldn't need the particlesPerThread parameter (although you still need to adjust blocks accordingly.) I know this might be hell to maintain but hey, if it's worth the performance... why not do it? This would be your new unrolled kernel:

__global__ void UpdateParticle(float* position, float* velocity, float frameTime)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if (idx < numParticles/4) { // the unrolled loop is 4 iterations
        position[idx] = position[idx] - velocity[idx] * frameTime * 0.001f;
        position[idx+1] = position[idx+1] - velocity[idx+1] * frameTime * 0.001f;
        position[idx+2] = position[idx+2] - velocity[idx+2] * frameTime * 0.001f;
        position[idx+3] = position[idx+3] - velocity[idx+3] * frameTime * 0.001f;
        ... // some more updates
    }
}

Conclusion

Unfortunately performace in CUDA involves a lot of tweaking, but I hope with all this background you'll be able to tweak it to fit your needs.

There could be further optimizations such as using __shared__ qualified memory (which is A LOT faster than reading from global-memory) but since you're executing so little operations on your data it would probably not be worth the hassle.

You might also want to precompute the frameTime * 0.001f bit in a register before anything else (just do float realTime = frameTime * 0.001f and use it instead) or even better: pass it already transformed from host code. It won't be a problem for such a small number of operations, but registers are also shared between blocks, so registers (any non-qualified variable inside your kernel, like idx in my examples) can be a bottleneck too. Bear it in mind!

Answer 3

Thank you so much for you in depth answers. These have provided a lot of useful information and I'm currently in the process of implementing the information you gave into my project. Though I currently have one issue I was wondering you could help me with.

Please note, I know that I haven't really modified the code you have given me, I am just trying to get the build first, then modify it proper. If I can get at least a single particle to fall down with this, I'll know I'm on the right track.

Right now, in my kernal class I have the following code:

    __global__ void UpdateParticle(float* position, float* velocity, float frameTime, int numParticles)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x; // Compute the index

    if (idx < numParticles) { // Is this index valid?
        position[idx] = position[idx] - velocity[idx] * frameTime * 0.001f;
      //  ... // some more updates
    }
}

Following that, in the .cpp file I have amended the particle update method as follows:

void ParticleSystemClass::UpdateParticles(float frameTime)
{
    int i;
    float* h_position; // Your host pointer. This holds the data (I assume it's already filled with the data.)
    float* d_position; // Your device pointer, we will allocate and fill this
    float* d_velocity;
    float* d_time;
    int threads_per_block = 128; // You should play with this value
    int blocks = m_maxParticles/threads_per_block + ( (m_maxParticles%threads_per_block)1:0 );
    const int N = 10;
    size_t size = N * sizeof(float);
    cudaMalloc( (void**)&d_position, m_maxParticles * sizeof(float) );
    cudaMemcpy( d_position, h_position, m_maxParticles * sizeof(float), cudaMemcpyHostToDevice);


    UpdateParticle_cuda(d_position, d_velocity,frameTime, m_maxParticles,blocks,  threads_per_block);

    return;
}

When I build the program, the only error I get comes from the following line:

    int blocks = m_maxParticles/threads_per_block + ( (m_maxParticles%threads_per_block)1:0 );

The error behind is being:

error C2059: syntax error: ')'

error C2143: syntax error: missing a ')' before 'constant'

From what I can see, there is no missing bracket in that statement, and no google search I do seems to yield any result for find a solution.

Is there something wrong with my what I have put so far, or I have got it totally wrong?

Please help me fix this bug.

Thanks.