1 //---------------------------------------------------------------------------//
2 // Copyright (c) 2014 Fabian Köhler <fabian2804@googlemail.com>
3 //
4 // Distributed under the Boost Software License, Version 1.0
5 // See accompanying file LICENSE_1_0.txt or copy at
6 // http://www.boost.org/LICENSE_1_0.txt
7 //
8 // See http://boostorg.github.com/compute for more information.
9 //---------------------------------------------------------------------------//
10
11 #include <iostream>
12
13 #define GL_GLEXT_PROTOTYPES
14 #ifdef __APPLE__
15 #include <OpenGL/gl.h>
16 #include <OpenGL/glext.h>
17 #else
18 #include <GL/gl.h>
19 #include <GL/glext.h>
20 #endif
21
22 #include <QtGlobal>
23 #if QT_VERSION >= 0x050000
24 #include <QtWidgets>
25 #else
26 #include <QtGui>
27 #endif
28 #include <QtOpenGL>
29 #include <QTimer>
30
31 #include <boost/program_options.hpp>
32 #include <boost/random/uniform_real_distribution.hpp>
33 #include <boost/random/mersenne_twister.hpp>
34
35 #ifndef Q_MOC_RUN
36 #include <boost/compute/system.hpp>
37 #include <boost/compute/algorithm.hpp>
38 #include <boost/compute/container/vector.hpp>
39 #include <boost/compute/interop/opengl.hpp>
40 #include <boost/compute/utility/source.hpp>
41 #endif // Q_MOC_RUN
42
43 namespace compute = boost::compute;
44 namespace po = boost::program_options;
45
46 using compute::uint_;
47 using compute::float4_;
48
49 const char source[] = BOOST_COMPUTE_STRINGIZE_SOURCE(
50 __kernel void updateVelocity(__global const float4* position, __global float4* velocity, float dt, uint N)
51 {
52 uint gid = get_global_id(0);
53
54 float4 r = { 0.0f, 0.0f, 0.0f, 0.0f };
55 float f = 0.0f;
56 for(uint i = 0; i != gid; i++) {
57 if(i != gid) {
58 r = position[i]-position[gid];
59 f = length(r)+0.001f;
60 f *= f*f;
61 f = dt/f;
62 velocity[gid] += f*r;
63 }
64 }
65 }
66 __kernel void updatePosition(__global float4* position, __global const float4* velocity, float dt)
67 {
68 uint gid = get_global_id(0);
69
70 position[gid].xyz += dt*velocity[gid].xyz;
71 }
72 );
73
74 class NBodyWidget : public QGLWidget
75 {
76 Q_OBJECT
77
78 public:
79 NBodyWidget(std::size_t particles, float dt, QWidget* parent = 0);
80 ~NBodyWidget();
81
82 void initializeGL();
83 void resizeGL(int width, int height);
84 void paintGL();
85 void updateParticles();
86 void keyPressEvent(QKeyEvent* event);
87
88 private:
89 QTimer* timer;
90
91 compute::context m_context;
92 compute::command_queue m_queue;
93 compute::program m_program;
94 compute::opengl_buffer m_position;
95 compute::vector<float4_>* m_velocity;
96 compute::kernel m_velocity_kernel;
97 compute::kernel m_position_kernel;
98
99 bool m_initial_draw;
100
101 const uint_ m_particles;
102 const float m_dt;
103 };
104
NBodyWidget(std::size_t particles,float dt,QWidget * parent)105 NBodyWidget::NBodyWidget(std::size_t particles, float dt, QWidget* parent)
106 : QGLWidget(parent), m_initial_draw(true), m_particles(particles), m_dt(dt)
107 {
108 // create a timer to redraw as fast as possible
109 timer = new QTimer(this);
110 connect(timer, SIGNAL(timeout()), this, SLOT(updateGL()));
111 timer->start(1);
112 }
113
~NBodyWidget()114 NBodyWidget::~NBodyWidget()
115 {
116 delete m_velocity;
117
118 // delete the opengl buffer
119 GLuint vbo = m_position.get_opengl_object();
120 glDeleteBuffers(1, &vbo);
121 }
122
initializeGL()123 void NBodyWidget::initializeGL()
124 {
125 // create context, command queue and program
126 m_context = compute::opengl_create_shared_context();
127 m_queue = compute::command_queue(m_context, m_context.get_device());
128 m_program = compute::program::create_with_source(source, m_context);
129 m_program.build();
130
131 // prepare random particle positions that will be transferred to the vbo
132 float4_* temp = new float4_[m_particles];
133 boost::random::uniform_real_distribution<float> dist(-0.5f, 0.5f);
134 boost::random::mt19937_64 gen;
135 for(size_t i = 0; i < m_particles; i++) {
136 temp[i][0] = dist(gen);
137 temp[i][1] = dist(gen);
138 temp[i][2] = dist(gen);
139 temp[i][3] = 1.0f;
140 }
141
142 // create an OpenGL vbo
143 GLuint vbo = 0;
144 glGenBuffers(1, &vbo);
145 glBindBuffer(GL_ARRAY_BUFFER, vbo);
146 glBufferData(GL_ARRAY_BUFFER, m_particles*sizeof(float4_), temp, GL_DYNAMIC_DRAW);
147
148 // create a OpenCL buffer from the vbo
149 m_position = compute::opengl_buffer(m_context, vbo);
150 delete[] temp;
151
152 // create buffer for velocities
153 m_velocity = new compute::vector<float4_>(m_particles, m_context);
154 compute::fill(m_velocity->begin(), m_velocity->end(), float4_(0.0f, 0.0f, 0.0f, 0.0f), m_queue);
155
156 // create compute kernels
157 m_velocity_kernel = m_program.create_kernel("updateVelocity");
158 m_velocity_kernel.set_arg(0, m_position);
159 m_velocity_kernel.set_arg(1, m_velocity->get_buffer());
160 m_velocity_kernel.set_arg(2, m_dt);
161 m_velocity_kernel.set_arg(3, m_particles);
162 m_position_kernel = m_program.create_kernel("updatePosition");
163 m_position_kernel.set_arg(0, m_position);
164 m_position_kernel.set_arg(1, m_velocity->get_buffer());
165 m_position_kernel.set_arg(2, m_dt);
166 }
resizeGL(int width,int height)167 void NBodyWidget::resizeGL(int width, int height)
168 {
169 // update viewport
170 glViewport(0, 0, width, height);
171 }
paintGL()172 void NBodyWidget::paintGL()
173 {
174 // clear buffer
175 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
176
177 // check if this is the first draw
178 if(m_initial_draw) {
179 // do not update particles
180 m_initial_draw = false;
181 } else {
182 // update particles
183 updateParticles();
184 }
185
186 // draw
187 glVertexPointer(4, GL_FLOAT, 0, 0);
188 glEnableClientState(GL_VERTEX_ARRAY);
189 glDrawArrays(GL_POINTS, 0, m_particles);
190 glFinish();
191 }
updateParticles()192 void NBodyWidget::updateParticles()
193 {
194 // enqueue kernels to update particles and make sure that the command queue is finished
195 compute::opengl_enqueue_acquire_buffer(m_position, m_queue);
196 m_queue.enqueue_1d_range_kernel(m_velocity_kernel, 0, m_particles, 0).wait();
197 m_queue.enqueue_1d_range_kernel(m_position_kernel, 0, m_particles, 0).wait();
198 m_queue.finish();
199 compute::opengl_enqueue_release_buffer(m_position, m_queue);
200 }
keyPressEvent(QKeyEvent * event)201 void NBodyWidget::keyPressEvent(QKeyEvent* event)
202 {
203 if(event->key() == Qt::Key_Escape) {
204 this->close();
205 }
206 }
207
main(int argc,char ** argv)208 int main(int argc, char** argv)
209 {
210 // parse command line arguments
211 po::options_description options("options");
212 options.add_options()
213 ("help", "show usage")
214 ("particles", po::value<uint_>()->default_value(1000), "number of particles")
215 ("dt", po::value<float>()->default_value(0.00001f), "width of each integration step");
216 po::variables_map vm;
217 po::store(po::parse_command_line(argc, argv, options), vm);
218 po::notify(vm);
219
220 if(vm.count("help") > 0) {
221 std::cout << options << std::endl;
222 return 0;
223 }
224
225 const uint_ particles = vm["particles"].as<uint_>();
226 const float dt = vm["dt"].as<float>();
227
228 QApplication app(argc, argv);
229 NBodyWidget nbody(particles, dt);
230
231 nbody.show();
232
233 return app.exec();
234 }
235
236 #include "nbody.moc"
237