1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
import numpy as np
import pyopencl as cl
mf = cl.mem_flags
from pyopencl.tools import get_gl_sharing_context_properties
import OpenGL.GL as gl
from OpenGL.arrays import vbo
from string import Template
class HardSphereSetup:
def __init__(self, radius, char_u, position, velocity):
self.radius = radius
self.char_u = char_u
self.position = position
self.velocity = velocity
self.n_particles = len(position[:,0])
def build_kernel(delta_t, n_particles, radius):
with open('kernel.cl', 'r') as kernel_src:
return Template(kernel_src.read()).substitute(
delta_t = delta_t,
n_particles = n_particles,
radius = radius)
def grid_of_random_velocity_particles(width, radius, u_scale):
np_position = np.ndarray((width**2, 2))
np_velocity = np.ndarray((width**2, 2))
grid = np.meshgrid(np.linspace(2*radius, 1-2*radius, width),
np.linspace(2*radius, 1-2*radius, width))
np_position[:,0] = grid[0].flatten()
np_position[:,1] = grid[1].flatten()
np_velocity[:,0] = u_scale*(-0.5 + np.random.random_sample((width**2,)))
np_velocity[:,1] = u_scale*(-0.5 + np.random.random_sample((width**2,)))
return (np_position, np_velocity)
class GasFlow:
def setup_cl(self):
self.platform = cl.get_platforms()[0]
if self.opengl:
self.context = cl.Context(
properties=[(cl.context_properties.PLATFORM, self.platform)] + get_gl_sharing_context_properties())
else:
self.context = cl.Context(
properties=[(cl.context_properties.PLATFORM, self.platform)])
self.queue = cl.CommandQueue(self.context)
self.program = cl.Program(self.context, self.kernel_src).build(
'-cl-single-precision-constant -cl-opt-disable')
self.cl_particle_velocity_a = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_particle_velocity)
self.cl_particle_velocity_b = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_particle_velocity)
if self.opengl:
self.gl_particle_position_a = vbo.VBO(data=self.np_particle_position, usage=gl.GL_DYNAMIC_DRAW, target=gl.GL_ARRAY_BUFFER)
self.gl_particle_position_a.bind()
self.gl_particle_position_b = vbo.VBO(data=self.np_particle_position, usage=gl.GL_DYNAMIC_DRAW, target=gl.GL_ARRAY_BUFFER)
self.gl_particle_position_b.bind()
self.cl_particle_position_a = cl.GLBuffer(self.context, mf.READ_WRITE, int(self.gl_particle_position_a))
self.cl_particle_position_b = cl.GLBuffer(self.context, mf.READ_WRITE, int(self.gl_particle_position_b))
else:
self.cl_particle_position_a = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_particle_position)
self.cl_particle_position_b = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_particle_position)
self.cl_last_collide = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_last_collide)
def __init__(self, setup, opengl = False, t_scale = 1.0):
self.np_particle_position = setup.position.astype(np.float32)
self.np_particle_velocity = setup.velocity.astype(np.float32)
self.n_particles = setup.n_particles
self.opengl = opengl
self.t_scale = t_scale
self.np_last_collide = np.ndarray((self.n_particles, 1), dtype=np.uint32)
self.np_last_collide[:,0] = self.n_particles
self.kernel_src = build_kernel(self.t_scale*setup.radius/setup.char_u, self.n_particles, setup.radius)
self.setup_cl()
self.tick = True
def evolve(self):
if self.opengl:
cl.enqueue_acquire_gl_objects(self.queue, [self.cl_particle_position_a, self.cl_particle_position_b])
if self.tick:
self.tick = False
kernelargs = (self.cl_particle_position_a, self.cl_particle_velocity_a, self.cl_particle_position_b, self.cl_particle_velocity_b, self.cl_last_collide)
else:
self.tick = True
kernelargs = (self.cl_particle_position_b, self.cl_particle_velocity_b, self.cl_particle_position_a, self.cl_particle_velocity_a, self.cl_last_collide)
self.program.evolve(self.queue, (self.n_particles,), None, *(kernelargs))
def gl_draw_particles(self):
gl.glEnableClientState(gl.GL_VERTEX_ARRAY)
if self.tick:
self.gl_particle_position_b.bind()
gl.glVertexPointer(2, gl.GL_FLOAT, 0, self.gl_particle_position_b)
else:
self.gl_particle_position_a.bind()
gl.glVertexPointer(2, gl.GL_FLOAT, 0, self.gl_particle_position_a)
gl.glDrawArrays(gl.GL_POINTS, 0, self.n_particles)
gl.glDisableClientState(gl.GL_VERTEX_ARRAY)
def get_velocities(self):
if self.tick:
cl.enqueue_copy(self.queue, self.np_particle_velocity, self.cl_particle_velocity_b).wait()
return self.np_particle_velocity
else:
cl.enqueue_copy(self.queue, self.np_particle_velocity, self.cl_particle_velocity_a).wait()
return self.np_particle_velocity
def get_positions(self):
if self.tick:
cl.enqueue_copy(self.queue, self.np_particle_position, self.cl_particle_position_b).wait()
return self.np_particle_position
else:
cl.enqueue_copy(self.queue, self.np_particle_position, self.cl_particle_position_a).wait()
return self.np_particle_position
|