Prototype port to 3D

What is it with OpenCL @ GPU and arrays of 3-dimensional vectors?!
Luckily I vaguely remembered encountering this problem before when
implementing my OpenCL LBM code but why on earth is the compiler not
warning about this? I guess the underlying reason has to do with the
data alignment requirements on the GPU but still...

5 files changed, 187 insertions, 52 deletions

diff --git a/boltzgas/initial_condition.py b/boltzgas/initial_condition.py
index 4d66263..2526566 100644
--- a/boltzgas/initial_condition.py
+++ b/boltzgas/initial_condition.py
@@ -1,15 +1,20 @@
 import numpy as np
 
 def grid_of_random_velocity_particles(width, radius, u_scale):
-    np_position = np.ndarray((width**2, 2))
-    np_velocity = np.ndarray((width**2, 2))
+    np_position = np.ndarray((width**3, 4))
+    np_velocity = np.ndarray((width**3, 4))
 
     grid = np.meshgrid(np.linspace(2*radius, 1-2*radius, width),
+                       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_position[:,2] = grid[2].flatten()
+    np_position[:,3] = 0
 
-    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,)))
+    np_velocity[:,0] = u_scale*(-0.5 + np.random.random_sample((width**3,)))
+    np_velocity[:,1] = u_scale*(-0.5 + np.random.random_sample((width**3,)))
+    np_velocity[:,2] = u_scale*(-0.5 + np.random.random_sample((width**3,)))
+    np_velocity[:,3] = 0
 
     return (np_position, np_velocity)
diff --git a/boltzgas/kernel.cl b/boltzgas/kernel.cl
index ad984bb..aabd1af 100644
--- a/boltzgas/kernel.cl
+++ b/boltzgas/kernel.cl
@@ -2,22 +2,44 @@
 #define N_PARTICLES $n_particles
 #define DELTA_T     $delta_t
 
+#define ENABLE_3D
+
 typedef float  scalar_t;
+
+#ifdef ENABLE_3D
+typedef float3 vec_t;
+typedef float4 data_vec_t;
+#define N_WALLS 6
+#else
 typedef float2 vec_t;
+typedef float2 data_vec_t;
+#define N_WALLS 4
+#endif
 
 #define SCALAR_MAX FLT_MAX
 
-__constant vec_t wall_normals[4] = {
-	(vec_t)(-1.,0.),
-	(vec_t)( 1.,0.),
-	(vec_t)(0., 1.),
-	(vec_t)(0.,-1.),
+__constant float4 wall_normals[N_WALLS] = {
+	(float4)(-1., 0., 0., 0.),
+	(float4)( 1., 0., 0., 0.),
+	(float4)( 0., 1., 0., 0.),
+	(float4)( 0.,-1., 0., 0.),
+	(float4)( 0., 0.,-1., 0.),
+	(float4)( 0., 0., 1., 0.),
 };
 
-__constant scalar_t wall_loc[4] = {
-	-RADIUS, 1.-RADIUS, 1.-RADIUS, -RADIUS
+
+__constant scalar_t wall_loc[N_WALLS] = {
+	-RADIUS, 1.-RADIUS, 1.-RADIUS, -RADIUS, -RADIUS, 1.-RADIUS
 };
 
+vec_t wall_normal(unsigned int iWall) {
+#ifdef ENABLE_3D
+	return wall_normals[iWall].xyz;
+#else
+	return wall_normals[iWall].xy;
+#endif
+}
+
 scalar_t solve_wall_collision(vec_t n, scalar_t loc, vec_t p, vec_t v) {
 	if (dot(n,v) > 0.) {
 		vec_t wall_v = dot(n,v) * n;
@@ -59,24 +81,59 @@ scalar_t solve_collision(vec_t p, vec_t v, vec_t p_, vec_t v_) {
 	}
 }
 
-__kernel void evolve(__global vec_t* pos_a,
-                     __global vec_t* vel_a,
-                     __global vec_t* pos_b,
-                     __global vec_t* vel_b,
+vec_t getPosition(__global data_vec_t* pos,
+                  unsigned int iParticle)
+{
+#ifdef ENABLE_3D
+	return pos[iParticle].xyz;
+#else
+	return pos[iParticle];
+#endif
+}
+
+vec_t getVelocity(__global data_vec_t* vel,
+                  unsigned int iParticle)
+{
+#ifdef ENABLE_3D
+	return vel[iParticle].xyz;
+#else
+	return vel[iParticle];
+#endif
+}
+
+void set(__global data_vec_t* pos,
+         __global data_vec_t* vel,
+         unsigned int iParticle,
+         vec_t p,
+         vec_t v)
+{
+#ifdef ENABLE_3D
+	pos[iParticle].xyz = p;
+	vel[iParticle].xyz = v;
+#else
+	pos[iParticle] = p;
+	vel[iParticle] = v;
+#endif
+}
+
+__kernel void evolve(__global data_vec_t* pos_a,
+                     __global data_vec_t* vel_a,
+                     __global data_vec_t* pos_b,
+                     __global data_vec_t* vel_b,
                      __global volatile unsigned int* last_collide)
 {
 	unsigned int i = get_global_id(0);
 
-	vec_t p = pos_a[i];
-	vec_t v = vel_a[i];
+	vec_t p = getPosition(pos_a, i);
+	vec_t v = getVelocity(vel_a, i);
 
 	unsigned int jParticle = N_PARTICLES;
 	scalar_t min2intersect = SCALAR_MAX;
 
 	for (unsigned int iParticle=0; iParticle < N_PARTICLES; ++iParticle) {
 		if (iParticle != i && !(last_collide[i] == iParticle && last_collide[iParticle] == i)) {
-			vec_t p_ = pos_a[iParticle];
-			vec_t v_ = vel_a[iParticle];
+			vec_t p_ = getPosition(pos_a, iParticle);
+			vec_t v_ = getVelocity(vel_a, iParticle);
 
 			scalar_t time2intersect = solve_collision(p, v, p_, v_);
 
@@ -90,8 +147,8 @@ __kernel void evolve(__global vec_t* pos_a,
 	unsigned int jWall = N_PARTICLES;
 	scalar_t min2wall = SCALAR_MAX;
 
-	for (unsigned int iWall=0; iWall < 4; ++iWall) {
-		scalar_t time2wall = solve_wall_collision(wall_normals[iWall], wall_loc[iWall], p, v);
+	for (unsigned int iWall=0; iWall < N_WALLS; ++iWall) {
+		scalar_t time2wall = solve_wall_collision(wall_normal(iWall), wall_loc[iWall], p, v);
 		if (time2wall < min2wall) {
 			min2wall = time2wall;
 			jWall = iWall;
@@ -101,33 +158,29 @@ __kernel void evolve(__global vec_t* pos_a,
 	if (min2intersect < DELTA_T) {
 		if (min2wall < min2intersect) {
 			p += min2wall * v;
-			v -= 2*dot(v,wall_normals[jWall])*wall_normals[jWall];
+			v -= 2*dot(v,wall_normal(jWall)) * wall_normal(jWall);
 			p += (DELTA_T - min2wall) * v;
 			last_collide[i] = N_PARTICLES;
 
-			pos_b[i] = p;
-			vel_b[i] = v;
+			set(pos_b, vel_b, i, p, v);
 		} else {
 			if (i < jParticle) {
-				vec_t p_ = pos_a[jParticle];
-				vec_t v_ = vel_a[jParticle];
+				vec_t p_ = getPosition(pos_a, jParticle);
+				vec_t v_ = getVelocity(vel_a, jParticle);
 
 				p  += min2intersect * v;
 				p_ += min2intersect * v_;
 
 				vec_t omega = normalize(p - p_);
 
-				v  -= dot(vel_a[i] - vel_a[jParticle], omega) * omega;
-				v_ -= dot(vel_a[jParticle] - vel_a[i], omega) * omega;
+				v  -= dot(getVelocity(vel_a, i) - getVelocity(vel_a, jParticle), omega) * omega;
+				v_ -= dot(getVelocity(vel_a, jParticle) - getVelocity(vel_a, i), omega) * omega;
 
 				p  += (DELTA_T - min2intersect) * v;
 				p_ += (DELTA_T - min2intersect) * v_;
 
-				pos_b[i] = p;
-				vel_b[i] = v;
-
-				pos_b[jParticle] = p_;
-				vel_b[jParticle] = v_;
+				set(pos_b, vel_b, i,         p,  v);
+				set(pos_b, vel_b, jParticle, p_, v_);
 
 				last_collide[i] = jParticle;
 				last_collide[jParticle] = i;
@@ -136,21 +189,20 @@ __kernel void evolve(__global vec_t* pos_a,
 	} else {
 		if (min2wall < DELTA_T) {
 			p += min2wall * v;
-			v -= 2*dot(v,wall_normals[jWall])*wall_normals[jWall];
+			v -= 2*dot(v,wall_normal(jWall)) * wall_normal(jWall);
 			p += (DELTA_T - min2wall) * v;
 			last_collide[i] = N_PARTICLES;
 		} else {
 			p += DELTA_T * v;
 		}
 
-		pos_b[i] = p;
-		vel_b[i] = v;
+		set(pos_b, vel_b, i, p, v);
 	}
 }
 
-__kernel void get_velocity_norms(__global vec_t* velocities, __global scalar_t* norms)
+__kernel void get_velocity_norms(__global data_vec_t* velocities, __global scalar_t* norms)
 {
 	unsigned int i = get_global_id(0);
 
-	norms[i] = length(velocities[i]);
+	norms[i] = length(getVelocity(velocities, i));
 }
diff --git a/boltzgas/simulation.py b/boltzgas/simulation.py
index d818bab..905a04f 100644
--- a/boltzgas/simulation.py
+++ b/boltzgas/simulation.py
@@ -95,10 +95,10 @@ class HardSphereSimulation:
 
         if self.tick:
             self.gl_particle_position_b.bind()
-            gl.glVertexPointer(2, gl.GL_FLOAT, 0, self.gl_particle_position_b)
+            gl.glVertexPointer(4, 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.glVertexPointer(4, gl.GL_FLOAT, 0, self.gl_particle_position_a)
 
         gl.glDrawArrays(gl.GL_POINTS, 0, self.n_particles)
         gl.glDisableClientState(gl.GL_VERTEX_ARRAY)
diff --git a/boltzgas/visual/view.py b/boltzgas/visual/view.py
index e2746e2..327214a 100644
--- a/boltzgas/visual/view.py
+++ b/boltzgas/visual/view.py
@@ -1,7 +1,7 @@
 from OpenGL.GL   import *
 from OpenGL.GLUT import *
 
-from pyrr import matrix44
+from pyrr import matrix44, quaternion
 
 import numpy as np
 
@@ -24,26 +24,31 @@ particle_shader = (
     """
         #version 430
 
-        layout (location=0) in vec2 particles;
+        layout (location=0) in vec4 particles;
 
         out vec3 color;
 
         uniform mat4 projection;
+        uniform mat4 rotation;
+
         uniform vec3 face_color;
         uniform vec3 trace_color;
         uniform uint trace_id;
 
         void main() {
-            gl_Position = projection * vec4(particles, 0., 1.);
+            gl_Position = projection * rotation * vec4(particles.xyz, 1.);
 
-            if (gl_VertexID == trace_id) {
+            if (particles.x < 0.0 || particles.x > 1.0 ||
+                particles.y < 0.0 || particles.y > 1.0 ||
+                particles.z < 0.0 || particles.z > 1.0
+            ) {
                 color = trace_color;
             } else {
                 color = face_color;
             }
         }
     """,
-    ['projection', 'face_color', 'trace_color', 'trace_id']
+    ['projection', 'rotation', 'face_color', 'trace_color', 'trace_id']
 )
 
 decoration_shader = (
@@ -106,6 +111,69 @@ texture_shader = (
     ['picture','projection','mixing']
 )
 
+
+class Projection:
+    def __init__(self, distance):
+        self.distance = distance
+        self.ratio    = 4./3.
+        self.update()
+
+    def update(self):
+        projection = matrix44.create_perspective_projection(20.0, self.ratio, 0.1, 5000.0)
+        look = matrix44.create_look_at(
+            eye    = [0, 0, -self.distance],
+            target = [0, 0, 0],
+            up     = [0,-1, 0])
+
+        self.matrix = np.matmul(look, projection)
+
+    def update_ratio(self, width, height, update_viewport = True):
+        if update_viewport:
+            glViewport(0,0,width,height)
+
+        self.ratio = width/height
+        self.update()
+
+    def update_distance(self, change):
+        self.distance += change
+        self.update()
+
+    def get(self):
+        return self.matrix
+
+class Rotation:
+    def __init__(self, shift, x = np.pi, z = np.pi):
+        self.matrix = matrix44.create_from_translation(shift),
+        self.rotation_x = quaternion.Quaternion()
+        self.update(x,z)
+
+    def shift(self, x, z):
+        self.matrix = np.matmul(
+            self.matrix,
+            matrix44.create_from_translation([x,0,z])
+        )
+        self.inverse_matrix = np.linalg.inv(self.matrix)
+
+    def update(self, x, z):
+        rotation_x = quaternion.Quaternion(quaternion.create_from_eulers([x,0,0]))
+        rotation_z = self.rotation_x.conjugate.cross(
+                quaternion.Quaternion(quaternion.create_from_eulers([0,0,z])))
+        self.rotation_x = self.rotation_x.cross(rotation_x)
+
+        self.matrix = np.matmul(
+            self.matrix,
+            matrix44.create_from_quaternion(rotation_z.cross(self.rotation_x))
+        )
+        self.inverse_matrix = np.linalg.inv(self.matrix)
+
+    def get(self):
+        return self.matrix
+
+    def get_inverse(self):
+        return self.inverse_matrix
+
+
+
 class View:
     def __init__(self, gas, decorations, windows):
         self.gas = gas
@@ -116,12 +184,17 @@ class View:
         self.particle_shader   = Shader(*particle_shader)
         self.decoration_shader = Shader(*decoration_shader)
 
+        self.projection3d = Projection(distance = 7)
+        self.rotation3d = Rotation([-1/2, -1/2, -1/2], 5*np.pi/4, np.pi/4)
+
     def reshape(self, width, height):
         glViewport(0,0,width,height)
 
         world_height = 1.4
         world_width = world_height / height * width
 
+        projection = Projection(10)
+
         projection  = matrix44.create_orthogonal_projection(-world_width/2, world_width/2, -world_height/2, world_height/2, -1, 1)
         translation = matrix44.create_from_translation([-1.05, -1.0/2, 0])
 
@@ -143,10 +216,15 @@ class View:
             window.display(self.texture_shader.uniform)
 
         self.particle_shader.use()
-        glUniformMatrix4fv(self.particle_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection))
+        #glUniformMatrix4fv(self.particle_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection))
+
+        glUniformMatrix4fv(self.particle_shader.uniform['projection'], 1, False, np.ascontiguousarray(self.projection3d.get()))
+        glUniformMatrix4fv(self.particle_shader.uniform['rotation'],   1, False, np.ascontiguousarray(self.rotation3d.get()))
+
         glUniform3f(self.particle_shader.uniform['face_color'],  1., 1., 1.)
         glUniform3f(self.particle_shader.uniform['trace_color'], 1., 0., 0.)
         glUniform1ui(self.particle_shader.uniform['trace_id'], -1)
+
         glEnable(GL_POINT_SPRITE)
         glPointSize(2*self.gas.radius*self.pixels_per_unit)
         self.gas.gl_draw_particles()
diff --git a/random_velocities.py b/random_velocities.py
index cf7a503..2d60f3b 100644
--- a/random_velocities.py
+++ b/random_velocities.py
@@ -4,20 +4,20 @@ from boltzgas import HardSphereSetup, HardSphereSimulation
 from boltzgas.initial_condition import grid_of_random_velocity_particles
 from boltzgas.visual import VelocityHistogram, ColoredBox, Tracer
 
-grid_width = 50
-radius = 0.002
+grid_width = 10
+radius = 0.005
 char_u = 1120
 
 position, velocity = grid_of_random_velocity_particles(grid_width, radius, char_u)
 
 config = HardSphereSetup(radius, char_u, position, velocity)
-gas = HardSphereSimulation(config, opengl = True)
+gas = HardSphereSimulation(config, opengl = True, t_scale=1)
 
-tracer = Tracer(gas, int((grid_width**2)/2+grid_width/2))
+#tracer = Tracer(gas, int((grid_width**2)/2+grid_width/2))
 histogram = VelocityHistogram(gas, [1.1,0], [1,1])
 
-decorations = [ ColoredBox([0,0], [1,1], (0.2,0.2,0.2)), tracer ]
-instruments = [ histogram, tracer ]
+decorations = [ ]
+instruments = [ histogram ]
 windows = [ histogram ]
 
 boltzgas.visualizer.simulate(config, gas, instruments, decorations, windows)