GPU Computing: Particle Simulation

Dipl.-Ing. Jan Novák∗ Dipl.-Inf. Gábor Liktor† Prof. Dr.-Ing. Carsten Dachsbacher‡

Abstract trade speed for quality. We will use a basic integration scheme
called Verlet integration.
In this assignment we will learn how to implement two simple par-
ticle systems on the GPU. In the first simulation the particles do not The Verlet velocity integration defines the new position x(t + ∆t)
interact with each other. However, we will allow particles to be re- and velocity v(t + ∆t) of a particle as:
moved and added to the system, hence we have to perform a stream
compaction that will remove the empty fields in the array. The sec- 1
x(t + ∆t) = x(t) + v(t)∆t + a(t)∆t2
ond assignment will then introduce local constraints that will affect 2
the movement of particles to simulate materials such as cloth. 1
v(t + ∆t) = v(t) + (a(t) + a(t + ∆t))∆t
2
The deadline for the assignment is on 22th of June.
The ∆t refers to a discrete time step that we take between two (cur-
1 Scientific Simulation rent and the new) simulation steps. a(t) refers to the acceleration
(i.e. first derivation of the velocity, second derivation of the posi-
Particle systems are widely used to simulate natural phenomena like tion) of the particle. As a thorough explanation of the underlying
smoke, fire, or water. In numerical simulation, we distinguish two motivation for using this integration scheme is beyond the scope of
approaches that conceptually differ in how they handle the simula- this text, we refer you to Wikipedia for a more detail description of
tion domain. the velocity Verlet integration.

Lagrangian Approach The Lagrangian approach is based on the 1.2 Collision Detection
discretization of the domain into a set of finite mass elements.
These particles are then allowed to move freely, or with respect to Since the particles can interact with the environment, e.g. collide
defined rules, through the environment. Each particle is usually de- with a wall, we need to simulate these interactions. Consider for
scribed by its position and velocity. According to our needs, it can instance a particle that is on one side of a wall at time t and on the
also have any other quantity such as mass, temperature, or radius. opposite side of the same wall at t + ∆t. Obviously, the veloc-
Methods based on tracking particles over time and altering their ity and acceleration of the particle have moved it through the wall.
quantities according to the prescribed laws are called particle Sys- We should try to detect such interactions and correct the position
tems. Since they directly track chunks of matter through the space, of the particle accordingly. The most robust approach is to advance
particle systems can easily guarantee the conservation of mass. On only by ∆t that is collision free: there is no collision within the
the other hand, handling of incompressibility (e.g. in simulation of whole system until ∆t. This is called continuous collision detec-
fluids) might be more difficult to achieve. tion, which requires computing the time to the next collision for
each particle, and setting ∆t to the minimum of these values. In
systems with frequent collisions, this can make the time step very
Eulerian Approach Instead of simulating the phenomena using small considerably slowing down the overall simulation.
chunks of matter, we can partition the domain into a set of small
symmetric cells: voxels. This approach is called Eulerian dis- We will take a simplified approach: we always advance by a con-
cretization and the main idea resides in allowing the matter to freely stant time step and correct for all previously occurring collisions in
move through the fixed grid, while tracing the simulated quantities. a post process. Furthermore, we will only account for one collision
In other words, we do not track the matter itself, but we simulate the during a single time step. Therefore, given the old and new particle
quantities at fixed positions in the domain. The discretization is in position, x(t) and x(t + ∆t), we will search for an intersection of
the simplest case a regular Cartesian grid, however, more advanced the line between these two points and all the objects in the scene. If
techniques often employ hierarchical or multi-resolution structures there is such an intersection, we should correct the position of the
to devote more computation and resolution to areas, where the particle using one of the approaches shown in Figure 1.
quantities change with higher frequency.
In this assignment we will use the Lagrangian approach simulating normal x(t) x(t)
a number of particles that do not interact with each other (first task), ∆t ∆t
v(t) v(t+∆t) v(t)
and a system where the particles represent a cloth and interact with v(t+∆t)
their neighbors under some constraints (e.g. springs). x(t+∆t) surface x(t+∆t) surface

1.1 Verlet Integration corrected corrected

v(t+∆t) x(t+∆t) new position v(t+∆t) x(t+∆t) new position

Having a system with the matter represented as particles, we want

original new position original new position
to find out the position of each particle at an arbitrary time. As
long as the system is just moderately complicated, we cannot use
the closed form solutions and have to solve the equations numer- Figure 1: Left: We correct the position of the particle by pulling it
ically. There are various approaches to numerical integration that back where the collision occurred. Right: We take the new position
of the particle and project it onto the surface using the surface nor-
∗ e-mail: [Link]@[Link] mal. The corrected velocity is in both cases computed by reflecting
† e-mail: [Link]@[Link] the original velocity about the normal of the surface.
‡ e-mail: dachsbacher@[Link]
2 OpenGL Interoperation • Kill old particles

During the simulation the particles will move along complex paths • Possibly create new particles
in 3-dimensional space and collide with solid surfaces. To evaluate • Store new particle data
the results, we need to provide real-time visual feedback about the
status of the simulation. OpenCL is suitable for general-purpose 3.1.1 Integration
programming, but we should employ one of the standard rendering
APIs to efficiently display 3D content on the screen. In order to compute the acceleration of the particle, use the
This assignment demonstrates the basics of interoperability. We gravitational acceleration and the acceleration defined by the
will use OpenCL to update the state of the simulated world (particle mass and the force at the current position of the particle.
positions, velocities, etc.) and OpenGL to display them. In this For fetching the force from the 3D texture use the func-
case, interoperability means that the same resources will be used tion read imagef(gForceField, sampler, lookUp),
in multiple contexts. For example, triplets of floats in an OpenCL where the gForceField and the sampler are parameters that
buffer memory object can be reinterpreted as vertex positions by the kernel obtains and the lookUp is a float4 with the first
OpenGL and displayed as a set of points in space. xyz componets specifying the position and the w defines the mip
level (in our case it should be set to 0). These are already initial-
OpenGL 2.0 defines buffer objects to hold rendering data. Based ized to perform the trilinear interpolation. Since you will have to
on their usage, we can talk about vertex (VBO), pixel (PBO) and perform some arithmetic operations on vectors, we advise you to
texture (TBO) buffer objects. If we want an OpenCL kernel to be use float4 data type for the 3-component position and velocity
able to modify the rendered geometry, we can create a buffer mem- with the fourth component set to 0. You can also define a custom
ory object from a VBO. As OpenCL and OpenGL coexist on the float3 type and write a few basic functions for them, however,
same device in parallel, there can be conflicts when accessing the you will not be able to define custom operators for them, as OpenCL
shared resources. Therefore, before executing OpenCL kernels that does not support C++ constructs.
use the shared buffers, the OpenCL must place a lock these buffers.
3.1.2 Collision Detection
As graphics programming using OpenGL is not covered by our
course, the implementation of the rendering part will be already In order to compute the collision of particles with the scene ob-
provided. To complete this assignment, you only need a basic un- jects, the kernel is also given a pointer to a global array with all
derstanding about the OpenGL context sharing. triangles in the scene. The triangles are stored as a triangle soup:
the triangles are stored in a consecutive chunk of memory, each tri-
3 Task 1: Simple Particle System angle is represented by three vertices, and each vertex is defined
as a float4 variable. Given the old and the new position of the
In this task we will implement a simple particle system that is driven particle, you will construct a ray segment and try to intersect this
by a force field. As there are no interactions or collisions between segment with each of these triangles. The straightforward solution
particles, the algorithm can be trivially parallelized using one thread is to iterate over all vertices in the global memory, construct a trian-
for each particle. We will use a 3D vector field to define the force gle, and call a function that will determine whether the ray segment
field within the simulation domain. This force field is loaded from intersects the triangle.
a file and uploaded to the GPU as a 3D texture. Using a 3D tex- Since all threads are iterating over the same values, there is a great
ture instead of a regular linear array has two advantages: first, we chance to use the local memory to cache the triangles. One of the
can use a 3D vector to conveniently address the 3D texture, second, possible approaches is to read one vertex from the triangle soup by
the hardware can automatically perform a trilinear interpolation of each thread and store it in the local memory. Then all threads can it-
the eight nearest neighbors. In other words, if we address a point erate over the cached triangles in the local memory and perform the
within the texture that is not exactly one of the discrete positions collision test. If the number of triangles is higher then the number
at which the texture samples the signal, we would have to load the of cached triangles, we have to repeat the process of loading and
eight nearest neighbors and perform a trilinear interpolation manu- testing the triangles multiple times. Notice, that if the number of
ally. As this is a frequently used operation in rendering, GPUs have threads is not a multiple of 3 (the number of vertices for one trian-
a hardware implementation that significantly speeds up the trilinear gle), it can happen that we will not read the whole last triangle: the
filtering. threads might load only one or two vertices of the triangle. Unless
we want to take a special care of this case, we have to make sure
3.1 Integration and Collision Detection that the number of threads within a work-group is divisible by 3.

Each particle in our system is defined using a position, velocity, We provide you with a LineTriangleIntersection func-
mass, and age. The first two characteristics are 3-component vec- tion, which computes an intersection of a line with a triangle. Your
tors, whereas the mass and age are scalars, therefore, we can pack task is to call this function for each triangle and find the closest in-
all the particle data into two float4 arrays. This will enable coa- tersection, as the particle should collide with the nearest triangle.
lesced accesses and minimize the fillrate of the application. Once you have it, you can adjust the new position and velocity of
the particle using one of the approaches illustrated in Figure 1.
In order to implement the Verlet integration and collision de-
tection, add you implementation in the Integrate kernel in 3.1.3 Removing and Adding Particles
[Link]. This kernel should perform the follow-
ing steps: So far the problem could be very easily parallelized. In order to
introduce a little bit of complexity, we will allow the particles to
• Load the particle data die and be reborn. After you account for the collisions, you should
• Perform the Verlet integration decrease the age of the particle and check if it is less or equal to
zero. If yes, the particle should be removed. In order to mark the
• Check for collisions particle as dead, use the gAlive buffer. Each record in this array
maps to exactly one particle (as in the case of the gPosAge and Particle data
gVelMass that store the position and the age, and the velocity and
the mass).
Flag array
We will also allow new particles to be born. However, we do not 1 0 1 1 0 1 1 1 1 1 0 1 1 1 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
want to just regenerate the dead particles. We want the new parti-
cles to be born whenever some criterion is met. For instance, we Parallel Prefix Sum
want to generate a new particle whenever another particle exceeds
some velocity. In such case, the fast particle will be virtually split Rank array
into two slower particles. Another example of generating newborns 0 1 1 2 3 3 4 5 6 7 8 8 9 10 11 1213131313131313131414141414141414
can be when a particle strongly bounces into an obstacle. The cri-
terion for splitting the particles is left up to you. The important fact
Compaction
is that each current particle can possibly generate a new particle.
Therefore, the arrays for holding the particle data are initialized to
have double size. If a particle k generates a new particle, the new- 1 0 1 1 0 1 1 1 1 1 0 1 1 1 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
born will be stored at position N + k, where N is the number of
0 1 1 2 3 3 4 5 6 7 8 8 9 10 11 1213131313131313131414141414141414
particles. Figure 2 shows an example, where three particles die and
one particle is split into two.

N 2N Figure 3: The upper half of the figure shows an example of a par-

ticle data array, flag buffer signaling the status of the particle, and
an outcome of the exclusive parallel prefix sum, the rank array. The
Dead particles Particle generating new particle New particle lower part of the figure demonstrates how the rank and flag arrays
are used to compact the particles.
1 0 1 1 0 1 1 1 1 1 0 1 1 1 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
Buffer of flags
the particle data with some random positions, mass, and age. It
Figure 2: An example of dead and newly created particles. The bot-
also prepares most of the device arrays that you will need. The
tom array shows the flags that should be written into the gAlive
provided structure of kernels is not mandatory, so you can adjust
array.
it to suit your implementation better. There are two scene files:
[Link] and [Link] that contain a box with
3.2 Stream Compaction some additional geometry. You can select the actual scene that is
loaded as the beginning of the initialization function.
Up to now all the computation could have been concentrated in a
single kernel. This kernel can possibly produce new particles that 3.4 Visualization
are in the second half of the array. In order to enable rendering and
further splitting of these particles, we need to compact the whole A particle simulated by the OpenCL kernel is mapped to a single
array and get rid of the empty records. This step can be performed vertex during OpenGL rendering. We render the particle system
using a parallel prefix sum. We will execute the PPS on the array by generating a sprite at the location of each particle. A sprite is
of flags gAlive that mark alive particles with 1 and the other with a quad, which always faces towards the camera. These sprites are
zero. For each particle, the exclusive PPS will give us the total num- then rendered with additive blending to the screen, giving the im-
ber of preceding living particles, which equals to the new position pression of small, point light sources (Figure 4).
of the particle in the compacted array. The process is demonstrated
in Figure 3.
In order to perform the stream compaction, you can use the PPS
that you implemented in the second assignment. Once you com-
pute the exclusive prefix sum you can compact the particles in an-
other kernel that loads the particle data from input array and us-
ing the gAlive and gRank write the data from gPosAgeIn
and gVelMassIn to the right positions of gPosAgeOut and
gVelMassOut. Since we do not write into all positions in the
output array, you should also execute a tiny kernel before the reor-
ganization to set the particle data to zero. Otherwise there might be
Figure 4: Left: 192K particles are simulated in a force field, col-
some old particle data that would be rendered.
liding with the surrounding geometry. The color coding shows the
dynamically born particles, generated using speed and height cri-
3.3 Implementation Details teria, at the top left corner of the image. Right: The force field is
visualized using line primitives as indicators.
Notice that if the number of particles is greater than N (because the
number of added particles is higher than the number of removed)
we will only use the first N particles. This is a the usual way of To better understand the behavior of the particles, we can modulate
handling overflows within the context of GPUs: we have a given the color of the rendered sprites as a function of their weight, age
budget of memory and we cannot extend it. or speed. This color coding can not only create interesting, colorful
visualizations, but can also help for debugging the code. For ex-
The skeleton of the assignment already performs most of the initial- ample, coloring the particles based on their age will clearly show
ization that you will need. It creates the force 3D texture, initializes newly generated particles in order to check their behavior. We do
not get into the details of visualization, but we explain the color- 4 Task 2: Cloth Simulation
coding mechanism of the vertex shader, [Link]:
Real-time cloth simulation is a popular physically motivated sim-
ulation even in computer games. In this assignment we will learn
uniform samplerBuffer tboSampler ;
how to implement a basic cloth that uses a spring model to main-
v o i d main ( ) { tain the structure of the material and achieve the cloth-like behavior.
The springs (or constraints) are used to mimic the thread-like struc-
/ / g e t p o s i t i o n and a g e ( bound a s v e r t e x b u f f e r ) ture of the material. The model is easy to parallelize, however, the
vec4 v = vec4 ( gl_Vertex ) ;
resulting cloth behaves somewhat more elastic than most of the ma-
vec4 position = vec4 ( v . xyz , 1 ) ; terials we know from the real life. Despite this drawback, it gained
f l o a t age = v . w ; popularity in game industry, since a basic cloth can be simulated
only with a few particles and springs.
/ / g e t s p e e d and mass ( bound a s a t e x t u r e b u f f e r )
v = texelFetchBuffer ( tboSampler , gl_VertexID ) ;

vec3 speed = v . xyz ;

f l o a t mass = v . w ;

/ / render the p a r t i c l e s using t h e i r speed

gl_FrontColor = colorCode ( length ( speed ) ) ;

[...]

gl_Position = gl_ModelViewProjectionMatrix ∗ position ;

}

Figure 5: Cloth simulation based on an explicit model defining the

This GLSL code (the standard shading language of OpenGL) gets structure via virtual springs between the particles.
executed for each particle vertex. Each vertex will be replaced by
a point sprite later on, but the vertex shader will define the posi- The implementation of the cloth simulation requires the following
tion and color of the sprite. The gPosAge buffer is mapped as kernels:
a vertex buffer, so the shader can get the position and age of the
particles using the gl Vertex built-in variable. The gVelMass • Integration kernel
buffer is used as a texture buffer (TBO), and fetched using the • Constraints kernel
texelFetchBuffer() command. This particular sample then • Collision detection kernel
color codes the particle using the length of its speed vector. The last
line of the shader is a standard OpenGL transformation from world • Kernel for recomputing normals
space (where the particles are simulated) to another coordinate sys- Each of the kernels is detailed in the following sections. In our
tem which OpenGL uses to project the particles to the screen. Feel implementation, we first executed the integration kernel and then
free to modify this shader to develop different meaningful visual- entered a loop that launched the constraints and collision detection
izations of particle behavior. kernel several times (e.g. two times the resolution of the cloth). At
the end we recomputed the normals to achieve more correct shad-
The application also displays the force field as a set of colored lines. ing of the cloth. The suggested order and structure of kernels is not
These lines are scattered inside the volume, showing the direction mandatory. You can adjust the kernels (for instance by concatenat-
and magnitude of the forces at sample locations. You can toggle the ing some of them into a single one), number of iterations, how often
rendering of the force field during simulation by hitting ’f’. you check for collisions, or the data that you exchange in between
Note that the application loads the collision geometry from kernels, if you can justify your decision by meaningful arguments.
a Wavefront OBJ file. That means, you can experiment
with different collision objects by replacing the input file 4.1 Integration
(CParticleSystem::InitResources) to a custom one. As
you increase the number of triangles in the collision object, you will The first executed kernel is responsible for moving the particles ac-
notice that the performance of your simulation will drop dramati- cording to all external forces (e.g. gravity or wind). In contrast to
cally. This is because each particle is tested against each triangle the previous task, we suggest using the Verlet position integration,
in the scene. Using a regular grid to spatially index the triangles which allows handling the correction of the position due to con-
based on their position would make the search for collisions a lot straints more conveniently. Since we want the cloth to behave like
more efficient. You can think of such extensions for the free-style it was attached on top to a bar, we should not trigger the computa-
assignment. tion for some of the particles in the very first row. This condition is
already provided in the code. Notice that the same condition should
also appear in the kernel for satisfying the constraints.
3.5 Evaluation
4.2 Satisfaction of Constraints
The total amount of points reserved for this task is 10:
In order to simulate the cloth, we will use a set of constraints that
• Verlet velocity integration (2 points). will act as springs between the particles. These springs will try
• Collision detection of particles with triangles (1.5 points) to preserve the initial spacing between particles, pushing them to-
stored in the local memory (1.5 points). gether if they are too far, and pulling them away if they are too
close. Figure 6 shows three different types of springs that we will
• Aging of particles (1 point). use. The structural constraints preserve the original adjacency,
shear constraints support the grid-like structure, and the bend con-
• Adding particles and correct stream compaction (4 points). straints affect the stiffness of the cloth.
 The only step left to finish the simulation part of the task is to make
the cloth interact with the solid sphere placed under the bar. To
resolve particle collisions with a sphere, you can use the approach
illustrated in Figure 1(right).
The collision detection with the sphere is simple: we have to de-
termine if any particle is closer to the sphere center than the sphere
radius. In case of a collision, we can push the particle to the surface
of the sphere using a vector in the radial direction. Note, that this
simple implementation of collision response will always keep the
particles on the surface of the sphere. Compared to the collision
detection in the particle system, this approach is not very robust: it
can happen that the sphere travels through the cloth if you push it
Figure 6: Different types of constraints (springs) that we use to hard enough. You can improve the detection akin to the previous
preserve the structure and cloth-like behavior. task, but as it should be straightforward, it is not mandatory.

4.4 Recomputation of Normals

The kernel for satisfying the constraints should be executed mul-
tiple times. During each iteration, we shift the particles a little in The last kernel in the execution flow does not add to the simulation
order to relax the constraints to get closer to the state with minimum itself, but to the rendering of the cloth. In order to correctly shade
energy. Since we cannot satisfy the constraints in one pass, we have the cloth, we need to know the normals of the surface. Vast ma-
to take an iterative approach and trigger the satisfy constraints ker- jority of shading models (e.g. the Phong shading model) requires a
nel. Since all the constraints are evaluated in parallel, we must not surface normal for each point to be shaded. The normal is used to
move one particle by more than d/2 during one iteration, where d find out the orientation of the surface towards the light source (and
is the distance between two particles in the rest state. If we did, two the viewer), which we must know in order to compute the amount
neighboring particles could possibly get attracted or repulsed too of reflected light. Figure 8 shows the difference between using the
much and the computation would ”blow up”. Therefore, we have to original normals and recomputing the normals at each time step.
weight the contribution of each appropriately, or clamp the change
in position to d/2.
Figure 7 depicts particles that act on a single particle in a horizontal
direction only. Since we do not want to offset the particle by more
than d/2, we need to conservatively add the forces; otherwise the
simulation might blow up. We included some default weights in
the code(WEIGHT ORTHO, WEIGHT DIAG, etc.), which you can
use to scale the contribution of each spring. Feel free to adjust the
values or use different approach. You can also adjust the bending
and shear constraints to change the stiffness of the cloth.
Figure 8: Left: Incorrect shading due to incorrect normals. Right:
Shading with correctly recomputed normals.

The normals should be computed in the following manner. For each

particle (cloth vertex), you have to read the eight nearest neigh-
bors and iterate over them in a loop. In each iteration, read two
subsequent vertices, construct two edges from the center vertex to
these two vertices, and compute the normalized cross product of
the edges. This will give you the normal of the artificial triangle
defined by the center and the two neighbor vertices. Accumulate
all the normals of the eight surrounding triangles and normalize the
sum to get the recomputed surface normal at the center vertex. Fig-
ure 9 illustrates the construction of the first two triangles and their
normals.
Figure 7: Neighboring particles acting on the center particle in the
horizontal direction only.
normal 1 normal 2

In order to efficiently to implement the constraint satisfaction effi- edge 2 edge 2

ciently, you should use the local memory in the same spirit as in
edge 1 edge 3
the non-separable image convolution assignment. The halo region
will in this case correspond to a two particles wide ring around the
processed tile.

4.3 Collision Detection

If you run the simulation at this point, the cloth should perform a Figure 9: Construction of triangles and normals during the recom-
realistic waving motion according to the applied forces, but keeping putation of per-vertex normals.
the distances among particles by satisfying the cloth constraints.
4.5 Evaluation

The total amount of points reserved for this task is 10:

• Verlet integration accounting for gravity and wind (2 points).
• Satisfaction of constraints (4 points).
• Collision detection with a sphere (2 points).
• Recomputation of normals (2 points).
There are several options to enhance the simulation. One of them
is for instance ripping of the cloth: if the springs between particles
stretch too much, they break and they are not considered anymore.
We are willing the give you 2 extra points for implementing this
behavior (or anything else that has similar complexity as the cloth
ripping).