Extract GL moments, particle buffers and add texture buffer

11 files changed, 426 insertions, 231 deletions

diff --git a/channel_2d_gl_interop.py b/channel_2d_gl_interop.py
index f85c29b..cae2f5b 100644
--- a/channel_2d_gl_interop.py
+++ b/channel_2d_gl_interop.py
@@ -2,6 +2,7 @@ import numpy
 from string import Template
 
 from simulation         import Lattice, Geometry
+from utility.opengl     import MomentsVertexBuffer
 from utility.particles  import Particles
 from symbolic.generator import LBM
 
@@ -14,6 +15,7 @@ from OpenGL.GL import shaders
 
 from pyrr import matrix44
 
+
 lattice_x = 480
 lattice_y = 300
 
@@ -174,10 +176,11 @@ lattice.apply_material_map(
     get_channel_material_map(lattice.geometry))
 lattice.sync_material()
 
+moments_vbo = MomentsVertexBuffer(lattice)
+
 particles = Particles(
-    lattice.context,
-    lattice.queue,
-    lattice.memory.float_type,
+    lattice,
+    moments_vbo,
     numpy.mgrid[
         lattice.geometry.size_x//20:2*lattice.geometry.size_x//20:100j,
         4*lattice.geometry.size_y//9:5*lattice.geometry.size_y//9:100000/100j
@@ -187,39 +190,27 @@ def on_display():
     for i in range(0,updates_per_frame):
         lattice.evolve()
 
-    lattice.collect_gl_moments()
+    moments_vbo.collect()
 
     for i in range(0,updates_per_frame):
-        lattice.update_gl_particles(particles, aging = False)
+        particles.update(aging = False)
 
     lattice.sync()
 
     glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
 
-    lattice.memory.gl_moments.bind()
-    glEnableClientState(GL_VERTEX_ARRAY)
-
     shaders.glUseProgram(shader_program)
     glUniformMatrix4fv(projection_id, 1, False, numpy.ascontiguousarray(projection))
-
-    glVertexPointer(4, GL_FLOAT, 0, lattice.memory.gl_moments)
-
+    moments_vbo.bind()
     glPointSize(point_size)
     glDrawArrays(GL_POINTS, 0, lattice.geometry.volume)
 
-    particles.gl_particles.bind()
-    glEnableClientState(GL_VERTEX_ARRAY)
-
     shaders.glUseProgram(particle_program)
     glUniformMatrix4fv(projection_id, 1, False, numpy.asfortranarray(projection))
-
-    glVertexPointer(4, GL_FLOAT, 0, particles.gl_particles)
-
+    particles.bind()
     glPointSize(point_size)
     glDrawArrays(GL_POINTS, 0, particles.count)
 
-    glDisableClientState(GL_VERTEX_ARRAY)
-
     glutSwapBuffers()
 
 def on_reshape(width, height):
diff --git a/channel_3d_gl_interop.py b/channel_3d_gl_interop.py
index a9633c4..e509e45 100644
--- a/channel_3d_gl_interop.py
+++ b/channel_3d_gl_interop.py
@@ -3,6 +3,7 @@ from string import Template
 
 from simulation         import Lattice, Geometry
 from utility.particles  import Particles
+from utility.opengl     import MomentsVertexBuffer
 from symbolic.generator import LBM
 
 import symbolic.D3Q27 as D3Q27
@@ -238,10 +239,11 @@ primitives   = list(map(lambda material: material[0], filter(lambda material: no
 lattice.apply_material_map(material_map)
 lattice.sync_material()
 
+moments_vbo = MomentsVertexBuffer(lattice)
+
 particles = Particles(
-    lattice.context,
-    lattice.queue,
-    lattice.memory.float_type,
+    lattice,
+    moments_vbo,
     numpy.mgrid[
         2*lattice.geometry.size_x//100:4*lattice.geometry.size_x//100:particle_count/10000j,
         lattice.geometry.size_y//16:15*lattice.geometry.size_y//16:100j,
@@ -256,10 +258,10 @@ def on_display():
     for i in range(0,updates_per_frame):
         lattice.evolve()
 
-    lattice.collect_gl_moments()
+    moments_vbo.collect()
 
     for i in range(0,updates_per_frame):
-        lattice.update_gl_particles(particles, aging = True)
+        particles.update(aging = True)
 
     lattice.sync()
 
@@ -268,13 +270,10 @@ def on_display():
     glEnable(GL_DEPTH_TEST)
     glDepthFunc(GL_LESS)
 
-    glEnableClientState(GL_VERTEX_ARRAY)
-    particles.gl_particles.bind()
-
     shaders.glUseProgram(particle_program)
     glUniformMatrix4fv(projection_id, 1, False, numpy.ascontiguousarray(projection))
     glUniformMatrix4fv(rotation_id, 1, False, numpy.ascontiguousarray(rotation.get()))
-    glVertexPointer(4, GL_FLOAT, 0, particles.gl_particles)
+    particles.bind()
     glEnable(GL_POINT_SMOOTH)
     glPointSize(point_size)
     glDrawArrays(GL_POINTS, 0, particles.count)
diff --git a/ldc_2d_gl_interop.py b/ldc_2d_gl_interop.py
index 631b5c8..b7223e0 100644
--- a/ldc_2d_gl_interop.py
+++ b/ldc_2d_gl_interop.py
@@ -13,6 +13,8 @@ from OpenGL.GL import shaders
 
 from pyrr import matrix44
 
+from utility.opengl import MomentsTexture
+
 lattice_x = 480
 lattice_y = 300
 
@@ -70,15 +72,32 @@ lbm = LBM(D2Q9)
 
 window = glut_window(fullscreen = False)
 
-vertex_shader = shaders.compileShader(Template("""
+vertex_shader = shaders.compileShader("""
 #version 430
 
-layout (location=0) in vec4 CellMoments;
-
-out vec3 color;
+layout (location=0) in vec4 vertex;
+                   out vec2 frag_pos;
 
 uniform mat4 projection;
 
+void main() {
+    gl_Position = projection * vertex;
+    frag_pos    = vertex.xy;
+}""", GL_VERTEX_SHADER)
+
+fragment_shader = shaders.compileShader(Template("""
+#version 430
+
+in vec2 frag_pos;
+
+uniform sampler2D moments;
+
+out vec4 result;
+
+vec2 unit(vec2 v) {
+    return vec2(v[0] / $size_x, v[1] / $size_y);
+}
+
 vec3 blueRedPalette(float x) {
     return mix(
         vec3(0.0, 0.0, 1.0),
@@ -87,39 +106,17 @@ vec3 blueRedPalette(float x) {
     );
 }
 
-vec2 fluidVertexAtIndex(uint i) {
-    const float y = floor(float(i) / $size_x);
-    return vec2(
-        i - $size_x*y,
-        y
-    );
-}
-
-void main() {
-    const vec2 idx = fluidVertexAtIndex(gl_VertexID);
-
-    gl_Position = projection * vec4(
-        idx.x,
-        idx.y,
-        0.,
-        1.
-    );
-
-    color = blueRedPalette(CellMoments[3] / $lid_speed);
-}""").substitute({
-    'size_x'   : lattice_x,
-    'lid_speed': lid_speed
-}), GL_VERTEX_SHADER)
-
-fragment_shader = shaders.compileShader("""
-#version 430
-
-in vec3 color;
-
 void main(){
-	gl_FragColor = vec4(color.xyz, 0.0);
-}""", GL_FRAGMENT_SHADER)
-
+    const vec2 sample_pos = unit(frag_pos);
+    const vec4 data = texture(moments, sample_pos);
+    result.a = 1.0;
+    result.rgb = blueRedPalette(data[3] / $lid_speed);
+}
+""").substitute({
+    "size_x": lattice_x,
+    "size_y": lattice_y,
+    "lid_speed": lid_speed
+}), GL_FRAGMENT_SHADER)
 
 shader_program = shaders.compileProgram(vertex_shader, fragment_shader)
 projection_id = shaders.glGetUniformLocation(shader_program, 'projection')
@@ -137,26 +134,28 @@ lattice.apply_material_map(
     get_cavity_material_map(lattice.geometry))
 lattice.sync_material()
 
+moments_texture = MomentsTexture(lattice)
+
+cube_vertices, cube_edges = lattice.geometry.wireframe()
+
 def on_display():
     for i in range(0,updates_per_frame):
         lattice.evolve()
 
-    lattice.collect_gl_moments()
+    moments_texture.collect()
 
     glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
 
-    lattice.memory.gl_moments.bind()
-    glEnableClientState(GL_VERTEX_ARRAY)
-
     shaders.glUseProgram(shader_program)
     glUniformMatrix4fv(projection_id, 1, False, numpy.asfortranarray(projection))
-
-    glVertexPointer(4, GL_FLOAT, 0, lattice.memory.gl_moments)
-
-    glPointSize(point_size)
-    glDrawArrays(GL_POINTS, 0, lattice.geometry.volume)
-
-    glDisableClientState(GL_VERTEX_ARRAY)
+    moments_texture.bind()
+
+    glBegin(GL_POLYGON)
+    glVertex(0,0,0)
+    glVertex(lattice.geometry.size_x,0,0)
+    glVertex(lattice.geometry.size_x,lattice.geometry.size_y,0)
+    glVertex(0,lattice.geometry.size_y,0)
+    glEnd()
 
     glutSwapBuffers()
 
diff --git a/ldc_3d_gl_interop.py b/ldc_3d_gl_interop.py
index 4065881..a8a10b6 100644
--- a/ldc_3d_gl_interop.py
+++ b/ldc_3d_gl_interop.py
@@ -2,6 +2,7 @@ import numpy
 from string import Template
 
 from simulation         import Lattice, Geometry
+from utility.opengl     import MomentsVertexBuffer
 from utility.particles  import Particles
 from symbolic.generator import LBM
 
@@ -171,10 +172,11 @@ lattice.apply_material_map(
     get_cavity_material_map(lattice.geometry))
 lattice.sync_material()
 
+moments_vbo = MomentsVertexBuffer(lattice)
+
 particles = Particles(
-    lattice.context,
-    lattice.queue,
-    lattice.memory.float_type,
+    lattice,
+    moments_vbo,
     numpy.mgrid[
         8*lattice.geometry.size_x//10:9*lattice.geometry.size_x//10:10j,
         lattice.geometry.size_y//10:9*lattice.geometry.size_y//10:particle_count/100j,
@@ -189,22 +191,19 @@ def on_display():
     for i in range(0,updates_per_frame):
         lattice.evolve()
 
-    lattice.collect_gl_moments()
+    moments_vbo.collect()
 
     for i in range(0,updates_per_frame):
-        lattice.update_gl_particles(particles, aging = True)
+        particles.update(aging = True)
 
     lattice.sync()
 
     glClear(GL_COLOR_BUFFER_BIT)
-    glEnableClientState(GL_VERTEX_ARRAY)
-
-    particles.gl_particles.bind()
 
     shaders.glUseProgram(particle_program)
     glUniformMatrix4fv(projection_id, 1, False, numpy.ascontiguousarray(projection))
     glUniformMatrix4fv(rotation_id, 1, False, numpy.ascontiguousarray(rotation.get()))
-    glVertexPointer(4, GL_FLOAT, 0, particles.gl_particles)
+    particles.bind()
     glEnable(GL_POINT_SMOOTH)
     glPointSize(point_size)
     glDrawArrays(GL_POINTS, 0, particles.count)
diff --git a/simulation.py b/simulation.py
index a274be4..de30e64 100644
--- a/simulation.py
+++ b/simulation.py
@@ -10,8 +10,6 @@ from mako.template import Template
 from pathlib import Path
 
 from pyopencl.tools import get_gl_sharing_context_properties
-import OpenGL.GL as gl
-from OpenGL.arrays import vbo
 
 class Geometry:
     def __init__(self, size_x, size_y, size_z = 1):
@@ -37,21 +35,32 @@ class Geometry:
             return (self.size_x-2, self.size_y-2, self.size_z-2)
 
     def wireframe(self):
-        return ([
-            [0          , 0          , 0          ],
-            [self.size_x, 0          , 0          ],
-            [self.size_x, self.size_y, 0          ],
-            [0          , self.size_y, 0          ],
-            [0          , 0          , self.size_z],
-            [self.size_x, 0          , self.size_z],
-            [self.size_x, self.size_y, self.size_z],
-            [0          , self.size_y, self.size_z]
-        ],
-        [
-            (0,1), (1,2), (2,3), (3,0),
-            (4,5), (5,6), (6,7), (7,4),
-            (0,4), (1,5), (2,6), (3,7)
-        ])
+        if self.size_z > 1:
+            return ([
+                [0          , 0          , 0          ],
+                [self.size_x, 0          , 0          ],
+                [self.size_x, self.size_y, 0          ],
+                [0          , self.size_y, 0          ],
+                [0          , 0          , self.size_z],
+                [self.size_x, 0          , self.size_z],
+                [self.size_x, self.size_y, self.size_z],
+                [0          , self.size_y, self.size_z]
+            ],
+            [
+                (0,1), (1,2), (2,3), (3,0),
+                (4,5), (5,6), (6,7), (7,4),
+                (0,4), (1,5), (2,6), (3,7)
+            ])
+        else:
+            return ([
+                [0          , 0          ],
+                [self.size_x, 0          ],
+                [self.size_x, self.size_y],
+                [0          , self.size_y],
+            ],
+            [
+                (0,1), (1,2), (2,3), (3,0)
+            ])
 
 def pad(n, m):
     return (n // m + min(1,n % m)) * m
@@ -103,14 +112,7 @@ class Memory:
         self.cl_pop_a = cl.Buffer(self.context, mf.READ_WRITE, size=self.pop_size)
         self.cl_pop_b = cl.Buffer(self.context, mf.READ_WRITE, size=self.pop_size)
 
-        if opengl:
-            self.np_moments = numpy.ndarray(shape=(self.volume, 4), dtype=self.float_type)
-            self.gl_moments = vbo.VBO(data=self.np_moments, usage=gl.GL_DYNAMIC_DRAW, target=gl.GL_ARRAY_BUFFER)
-            self.gl_moments.bind()
-            self.cl_gl_moments  = cl.GLBuffer(self.context, mf.READ_WRITE, int(self.gl_moments))
-        else:
-            self.cl_moments  = cl.Buffer(self.context, mf.WRITE_ONLY, size=self.moments_size)
-
+        self.cl_moments  = cl.Buffer(self.context, mf.WRITE_ONLY, size=self.moments_size)
         self.cl_material = cl.Buffer(self.context, mf.READ_ONLY, size=self.volume * numpy.int32(0).nbytes)
 
     def gid(self, x, y, z = 0):
@@ -244,28 +246,3 @@ class Lattice:
         cl.enqueue_copy(self.queue, moments, self.memory.cl_moments).wait();
 
         return moments
-
-    def collect_gl_moments(self):
-        cl.enqueue_acquire_gl_objects(self.queue, [self.memory.cl_gl_moments])
-
-        if self.tick:
-            self.program.collect_gl_moments(
-                self.queue, self.grid.size(), self.layout, self.memory.cl_pop_b, self.memory.cl_material, self.memory.cl_gl_moments)
-        else:
-            self.program.collect_gl_moments(
-                self.queue, self.grid.size(), self.layout, self.memory.cl_pop_a, self.memory.cl_material, self.memory.cl_gl_moments)
-
-    def update_gl_particles(self, particles, aging = False):
-        cl.enqueue_acquire_gl_objects(self.queue, [particles.cl_gl_particles])
-
-        if aging:
-            age = numpy.float32(0.000006)
-        else:
-            age = numpy.float32(0.0)
-
-        self.program.update_particles(
-            self.queue, (particles.count,1), None,
-            self.memory.cl_gl_moments,
-            self.memory.cl_material,
-            particles.cl_gl_particles, particles.cl_init_particles,
-            age)
diff --git a/template/kernel.mako b/template/kernel.mako
index a57c2ec..025baf5 100644
--- a/template/kernel.mako
+++ b/template/kernel.mako
@@ -98,82 +98,3 @@ __kernel void collect_moments(__global __read_only  ${float_type}* f,
     moments[${pop_offset(i)} + gid] = ${ccode(expr.rhs)};
 % endfor
 }
-
-__kernel void collect_gl_moments(__global __read_only  ${float_type}* f,
-                                 __global __read_only  int* material,
-                                 __global __write_only float4* moments)
-{
-    const unsigned int gid = ${gid()};
-
-    __global __read_only ${float_type}* preshifted_f = f + gid;
-
-% for i in range(0,descriptor.q):
-    const ${float_type} f_curr_${i} = preshifted_f[${pop_offset(i)}];
-% endfor
-
-% for i, expr in enumerate(moments_subexpr):
-    const ${float_type} ${expr[0]} = ${ccode(expr[1])};
-% endfor
-
-    float4 data;
-
-    if (material[gid] == 1) {
-% if descriptor.d == 2:
-      data.x = ${ccode(moments_assignment[0].rhs)};
-      data.y = ${ccode(moments_assignment[1].rhs)};
-      data.z = ${ccode(moments_assignment[2].rhs)};
-      data.w = sqrt(data.y*data.y + data.z*data.z);
-% elif descriptor.d == 3:
-      data.x = ${ccode(moments_assignment[0].rhs)};
-      data.y = ${ccode(moments_assignment[1].rhs)};
-      data.z = ${ccode(moments_assignment[2].rhs)};
-      data.w = ${ccode(moments_assignment[3].rhs)};
-% endif
-    } else {
-      data.x = 0.0;
-      data.y = 0.0;
-      data.z = 0.0;
-      data.w = -material[gid];
-    }
-
-    moments[gid] = data;
-}
-
-__kernel void update_particles(__global __read_only  float4* moments,
-                               __global __read_only  int*    material,
-                               __global __write_only float4* particles,
-                               __global __read_only  float4* init_particles,
-                               float aging)
-{
-  const unsigned int pid = get_global_id(0);
-
-  float4 particle = particles[pid];
-
-% if descriptor.d == 2:
-  const unsigned int gid = floor(particle.y)*${memory.size_x} + floor(particle.x);
-% elif descriptor.d == 3:
-  const unsigned int gid = floor(particle.z)*${memory.size_x*memory.size_y} + floor(particle.y)*${memory.size_x} + floor(particle.x);
-% endif
-
-  const float4 moment = moments[gid];
-
-  if (material[gid] == 1) {
-    particle.x += moment.y;
-    particle.y += moment.z;
-% if descriptor.d == 2:
-    particle.w += min(particle.x, particle.y) * aging;
-% elif descriptor.d == 3:
-    particle.z += moment.w;
-    float dy = (particle.y-${geometry.size_y/2.0});
-    float dz = (particle.z-${geometry.size_z/2.0});
-    dy *= dy;
-    dz *= dz;
-    particle.w = 10.0*sqrt(moment.y*moment.y+moment.z*moment.z+moment.w*moment.w);
-% endif
-  } else {
-    particle.xyz = init_particles[pid].xyz;
-    particle.w   = 0.0;
-  }
-
-  particles[pid] = particle;
-}
diff --git a/template/opengl.mako b/template/opengl.mako
new file mode 100644
index 0000000..6b09f69
--- /dev/null
+++ b/template/opengl.mako
@@ -0,0 +1,100 @@
+<%
+def gid():
+    return {
+        2: 'get_global_id(1)*%d + get_global_id(0)' % memory.size_x,
+        3: 'get_global_id(2)*%d + get_global_id(1)*%d + get_global_id(0)' % (memory.size_x*memory.size_y, memory.size_x)
+    }.get(descriptor.d)
+
+def pop_offset(i):
+    return i * memory.volume
+
+def moments_cell():
+    return {
+        2: '(int2)(get_global_id(0), get_global_id(1))',
+        3: '(int4)(get_global_id(0), get_global_id(1), get_global_id(2), 0)'
+    }.get(descriptor.d)
+%>
+
+__kernel void collect_gl_moments(__global __read_only  ${float_type}* f,
+                                 __global __read_only  int* material,
+                                 __global __write_only float4* moments)
+{
+    const unsigned int gid = ${gid()};
+
+    __global __read_only ${float_type}* preshifted_f = f + gid;
+
+% for i in range(0,descriptor.q):
+    const ${float_type} f_curr_${i} = preshifted_f[${pop_offset(i)}];
+% endfor
+
+% for i, expr in enumerate(moments_subexpr):
+    const ${float_type} ${expr[0]} = ${ccode(expr[1])};
+% endfor
+
+    float4 data;
+
+    if (material[gid] == 1) {
+% if descriptor.d == 2:
+      data.x = ${ccode(moments_assignment[0].rhs)};
+      data.y = ${ccode(moments_assignment[1].rhs)};
+      data.z = ${ccode(moments_assignment[2].rhs)};
+      data.w = sqrt(data.y*data.y + data.z*data.z);
+% elif descriptor.d == 3:
+      data.x = ${ccode(moments_assignment[0].rhs)};
+      data.y = ${ccode(moments_assignment[1].rhs)};
+      data.z = ${ccode(moments_assignment[2].rhs)};
+      data.w = ${ccode(moments_assignment[3].rhs)};
+% endif
+    } else {
+      data.x = 0.0;
+      data.y = 0.0;
+      data.z = 0.0;
+      data.w = -material[gid];
+    }
+
+    moments[gid] = data;
+}
+
+__kernel void collect_gl_moments_to_texture(__global __read_only  ${float_type}* f,
+                                            __global __read_only  int* material,
+% if descriptor.d == 2:
+                                            __write_only image2d_t moments)
+% elif descriptor.d == 3:
+                                            __write_only image3d_t moments)
+% endif
+{
+    const unsigned int gid = ${gid()};
+
+    __global __read_only ${float_type}* preshifted_f = f + gid;
+
+% for i in range(0,descriptor.q):
+    const ${float_type} f_curr_${i} = preshifted_f[${pop_offset(i)}];
+% endfor
+
+% for i, expr in enumerate(moments_subexpr):
+    const ${float_type} ${expr[0]} = ${ccode(expr[1])};
+% endfor
+
+    float4 data;
+
+    if (material[gid] != 5) {
+% if descriptor.d == 2:
+      data.x = ${ccode(moments_assignment[0].rhs)};
+      data.y = ${ccode(moments_assignment[1].rhs)};
+      data.z = ${ccode(moments_assignment[2].rhs)};
+      data.w = sqrt(data.y*data.y + data.z*data.z);
+% elif descriptor.d == 3:
+      data.x = ${ccode(moments_assignment[0].rhs)};
+      data.y = ${ccode(moments_assignment[1].rhs)};
+      data.z = ${ccode(moments_assignment[2].rhs)};
+      data.w = ${ccode(moments_assignment[3].rhs)};
+% endif
+    } else {
+      data.x = 0.0;
+      data.y = 0.0;
+      data.z = 0.0;
+      data.w = 1.0;
+    }
+
+    write_imagef(moments, ${moments_cell()}, data);
+}
diff --git a/template/particles.mako b/template/particles.mako
new file mode 100644
index 0000000..6235307
--- /dev/null
+++ b/template/particles.mako
@@ -0,0 +1,38 @@
+__kernel void update_particles(__global __read_only  float4* moments,
+                               __global __read_only  int*    material,
+                               __global __write_only float4* particles,
+                               __global __read_only  float4* init_particles,
+                               float aging)
+{
+  const unsigned int pid = get_global_id(0);
+
+  float4 particle = particles[pid];
+
+% if descriptor.d == 2:
+  const unsigned int gid = floor(particle.y)*${memory.size_x} + floor(particle.x);
+% elif descriptor.d == 3:
+  const unsigned int gid = floor(particle.z)*${memory.size_x*memory.size_y} + floor(particle.y)*${memory.size_x} + floor(particle.x);
+% endif
+
+  const float4 moment = moments[gid];
+
+  if (material[gid] == 1 && particle.w < 1.0) {
+    particle.x += moment.y;
+    particle.y += moment.z;
+% if descriptor.d == 2:
+    particle.w += min(particle.x, particle.y) * aging;
+% elif descriptor.d == 3:
+    particle.z += moment.w;
+    float dy = (particle.y-${geometry.size_y/2.0});
+    float dz = (particle.z-${geometry.size_z/2.0});
+    dy *= dy;
+    dz *= dz;
+    particle.w = 10.0*sqrt(moment.y*moment.y+moment.z*moment.z+moment.w*moment.w);
+% endif
+  } else {
+    particle.xyz = init_particles[pid].xyz;
+    particle.w   = particle.w-1.0;
+  }
+
+  particles[pid] = particle;
+}
diff --git a/trugfeuer_2d_gl_interop.py b/trugfeuer_2d_gl_interop.py
index 6d66efa..969fe4f 100644
--- a/trugfeuer_2d_gl_interop.py
+++ b/trugfeuer_2d_gl_interop.py
@@ -2,6 +2,7 @@ import numpy
 from string import Template
 
 from simulation         import Lattice, Geometry
+from utility.opengl     import MomentsVertexBuffer
 from utility.particles  import Particles
 from symbolic.generator import LBM
 
@@ -174,10 +175,11 @@ lattice.apply_material_map(
     get_channel_material_map(lattice.geometry))
 lattice.sync_material()
 
+moments_vbo = MomentsVertexBuffer(lattice)
+
 particles = Particles(
-    lattice.context,
-    lattice.queue,
-    lattice.memory.float_type,
+    lattice,
+    moments_vbo,
     numpy.mgrid[
         4*lattice.geometry.size_x//9:5*lattice.geometry.size_x//9:particle_count/100j,
         lattice.geometry.size_y//20:2*lattice.geometry.size_y//20:100j,
@@ -187,30 +189,25 @@ def on_display():
     for i in range(0,updates_per_frame):
         lattice.evolve()
 
-    lattice.collect_gl_moments()
+    moments_vbo.collect()
 
     for i in range(0,updates_per_frame):
-        lattice.update_gl_particles(particles, aging = True)
+        particles.update(aging = True)
 
     lattice.sync()
 
     glClear(GL_COLOR_BUFFER_BIT)
-    glEnableClientState(GL_VERTEX_ARRAY)
-
-    lattice.memory.gl_moments.bind()
 
     shaders.glUseProgram(moment_program)
     glUniformMatrix4fv(projection_id, 1, False, numpy.ascontiguousarray(projection))
-    glVertexPointer(4, GL_FLOAT, 0, lattice.memory.gl_moments)
+    moments_vbo.bind()
     glPointSize(point_size)
     glDisable(GL_POINT_SMOOTH)
     glDrawArrays(GL_POINTS, 0, lattice.geometry.volume)
 
-    particles.gl_particles.bind()
-
     shaders.glUseProgram(particle_program)
     glUniformMatrix4fv(projection_id, 1, False, numpy.ascontiguousarray(projection))
-    glVertexPointer(4, GL_FLOAT, 0, particles.gl_particles)
+    particles.bind()
     glPointSize(point_size)
     glEnable(GL_POINT_SMOOTH)
     glDrawArrays(GL_POINTS, 0, particles.count)
diff --git a/utility/opengl.py b/utility/opengl.py
new file mode 100644
index 0000000..afbfa00
--- /dev/null
+++ b/utility/opengl.py
@@ -0,0 +1,136 @@
+import pyopencl as cl
+mf = cl.mem_flags
+
+import numpy
+import sympy
+
+from mako.template import Template
+from pathlib import Path
+
+from OpenGL.GL import *
+from OpenGL.arrays import vbo
+
+class MomentsVertexBuffer:
+    def __init__(self, lattice):
+        self.lattice = lattice
+
+        self.np_moments = numpy.ndarray(shape=(self.lattice.memory.volume, 4), dtype=self.lattice.memory.float_type)
+        self.gl_moments = vbo.VBO(data=self.np_moments, usage=GL_DYNAMIC_DRAW, target=GL_ARRAY_BUFFER)
+        self.gl_moments.bind()
+        self.cl_gl_moments  = cl.GLBuffer(self.lattice.context, mf.READ_WRITE, int(self.gl_moments))
+
+        self.build_kernel()
+
+    def build_kernel(self):
+        program_src = Template(filename = str(Path(__file__).parent/'../template/opengl.mako')).render(
+            descriptor = self.lattice.descriptor,
+            geometry   = self.lattice.geometry,
+            memory     = self.lattice.memory,
+
+            moments_subexpr    = self.lattice.moments[0],
+            moments_assignment = self.lattice.moments[1],
+            collide_subexpr    = self.lattice.collide[0],
+            collide_assignment = self.lattice.collide[1],
+
+            float_type = self.lattice.float_type[1],
+
+            ccode = sympy.ccode
+        )
+        self.program = cl.Program(self.lattice.context, program_src).build(self.lattice.compiler_args)
+
+    def bind(self):
+        self.gl_moments.bind()
+        glEnableClientState(GL_VERTEX_ARRAY)
+        glVertexPointer(4, GL_FLOAT, 0, self.gl_moments)
+
+    def collect(self):
+        cl.enqueue_acquire_gl_objects(self.lattice.queue, [self.cl_gl_moments])
+
+        if self.lattice.tick:
+            self.program.collect_gl_moments(
+                self.lattice.queue,
+                self.lattice.grid.size(),
+                self.lattice.layout,
+                self.lattice.memory.cl_pop_b,
+                self.lattice.memory.cl_material,
+                self.cl_gl_moments)
+        else:
+            self.program.collect_gl_moments(
+                self.lattice.queue,
+                self.lattice.grid.size(),
+                self.lattice.layout,
+                self.lattice.memory.cl_pop_a,
+                self.lattice.memory.cl_material,
+                self.cl_gl_moments)
+
+
+class MomentsTexture:
+    def __init__(self, lattice):
+        self.lattice = lattice
+        self.gl_texture_buffer = numpy.ndarray(shape=(self.lattice.memory.volume, 4), dtype=self.lattice.memory.float_type)
+        self.gl_texture_buffer[:,:] = 0.0
+
+        self.gl_moments = glGenTextures(1)
+        self.gl_texture_type = {2: GL_TEXTURE_2D, 3: GL_TEXTURE_3D}.get(self.lattice.descriptor.d)
+        glBindTexture(self.gl_texture_type, self.gl_moments)
+
+        if self.gl_texture_type == GL_TEXTURE_3D:
+            glTexImage3D(self.gl_texture_type, 0, GL_RGBA, self.lattice.memory.size_x, self.lattice.memory.size_y, self.lattice.memory.size_z, 0, GL_RGBA, GL_FLOAT, self.gl_texture_buffer)
+            glTexParameteri(self.gl_texture_type, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
+            glTexParameteri(self.gl_texture_type, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
+            glTexParameteri(self.gl_texture_type, GL_TEXTURE_WRAP_T,     GL_CLAMP_TO_EDGE)
+            glTexParameteri(self.gl_texture_type, GL_TEXTURE_WRAP_S,     GL_CLAMP_TO_EDGE)
+            glTexParameteri(self.gl_texture_type, GL_TEXTURE_WRAP_R,     GL_CLAMP_TO_EDGE)
+            self.cl_gl_moments  = cl.GLTexture(self.lattice.context, mf.READ_WRITE, self.gl_texture_type, 0, self.gl_moments, 3)
+        elif self.gl_texture_type == GL_TEXTURE_2D:
+            glTexImage2D(self.gl_texture_type, 0, GL_RGBA, self.lattice.memory.size_x, self.lattice.memory.size_y, 0, GL_RGBA, GL_FLOAT, self.gl_texture_buffer)
+            glTexParameteri(self.gl_texture_type, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
+            glTexParameteri(self.gl_texture_type, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
+            glTexParameteri(self.gl_texture_type, GL_TEXTURE_WRAP_T,     GL_CLAMP_TO_EDGE)
+            glTexParameteri(self.gl_texture_type, GL_TEXTURE_WRAP_S,     GL_CLAMP_TO_EDGE)
+            self.cl_gl_moments  = cl.GLTexture(self.lattice.context, mf.READ_WRITE, self.gl_texture_type, 0, self.gl_moments, 2)
+
+        self.build_kernel()
+
+    def build_kernel(self):
+        program_src = Template(filename = str(Path(__file__).parent/'../template/opengl.mako')).render(
+            descriptor = self.lattice.descriptor,
+            geometry   = self.lattice.geometry,
+            memory     = self.lattice.memory,
+
+            moments_subexpr    = self.lattice.moments[0],
+            moments_assignment = self.lattice.moments[1],
+            collide_subexpr    = self.lattice.collide[0],
+            collide_assignment = self.lattice.collide[1],
+
+            float_type = self.lattice.float_type[1],
+
+            ccode = sympy.ccode
+        )
+        self.program = cl.Program(self.lattice.context, program_src).build(self.lattice.compiler_args)
+
+    def bind(self, location = GL_TEXTURE0):
+        glEnable(self.gl_texture_type)
+        glActiveTexture(location);
+        glBindTexture(self.gl_texture_type, self.gl_moments)
+
+    def collect(self):
+        cl.enqueue_acquire_gl_objects(self.lattice.queue, [self.cl_gl_moments])
+
+        if self.lattice.tick:
+            self.program.collect_gl_moments_to_texture(
+                self.lattice.queue,
+                self.lattice.grid.size(),
+                self.lattice.layout,
+                self.lattice.memory.cl_pop_b,
+                self.lattice.memory.cl_material,
+                self.cl_gl_moments)
+        else:
+            self.program.collect_gl_moments_to_texture(
+                self.lattice.queue,
+                self.lattice.grid.size(),
+                self.lattice.layout,
+                self.lattice.memory.cl_pop_a,
+                self.lattice.memory.cl_material,
+                self.cl_gl_moments)
+
diff --git a/utility/particles.py b/utility/particles.py
index a7ef777..4b770bd 100644
--- a/utility/particles.py
+++ b/utility/particles.py
@@ -3,16 +3,23 @@ mf = cl.mem_flags
 
 import numpy
 
+from mako.template import Template
+from pathlib import Path
+
 import OpenGL.GL as gl
 from OpenGL.arrays import vbo
 
 class Particles:
-    def __init__(self, context, queue, float_type, grid):
-        self.context = context
+    def __init__(self, lattice, moments, grid):
+        self.lattice = lattice
+        self.context = self.lattice.context
+        self.queue   = self.lattice.queue
+        self.float_type = self.lattice.memory.float_type
+        self.moments = moments
         self.count = len(grid)
 
-        self.np_particles = numpy.ndarray(shape=(self.count, 4), dtype=float_type)
-        self.np_init_particles = numpy.ndarray(shape=(self.count, 4), dtype=float_type)
+        self.np_particles = numpy.ndarray(shape=(self.count, 4), dtype=self.float_type)
+        self.np_init_particles = numpy.ndarray(shape=(self.count, 4), dtype=self.float_type)
 
         if len(grid[0,:]) == 2:
             self.np_particles[:,0:2] = grid
@@ -28,5 +35,36 @@ class Particles:
         self.gl_particles.bind()
         self.cl_gl_particles = cl.GLBuffer(self.context, mf.READ_WRITE, int(self.gl_particles))
 
-        self.cl_init_particles = cl.Buffer(context, mf.READ_ONLY, size=self.count * 4*numpy.float32(0).nbytes)
-        cl.enqueue_copy(queue, self.cl_init_particles, self.np_init_particles).wait();
+        self.cl_init_particles = cl.Buffer(self.context, mf.READ_ONLY, size=self.count * 4*numpy.float32(0).nbytes)
+        cl.enqueue_copy(self.queue, self.cl_init_particles, self.np_init_particles).wait();
+
+        self.build_kernel()
+
+    def build_kernel(self):
+        program_src = Template(filename = str(Path(__file__).parent/'../template/particles.mako')).render(
+            descriptor = self.lattice.descriptor,
+            geometry   = self.lattice.geometry,
+            memory     = self.lattice.memory,
+            float_type = self.float_type,
+        )
+        self.program = cl.Program(self.lattice.context, program_src).build(self.lattice.compiler_args)
+
+    def bind(self):
+        gl.glEnableClientState(gl.GL_VERTEX_ARRAY)
+        self.gl_particles.bind()
+        gl.glVertexPointer(4, gl.GL_FLOAT, 0, self.gl_particles)
+
+    def update(self, aging = False):
+        cl.enqueue_acquire_gl_objects(self.queue, [self.cl_gl_particles])
+
+        if aging:
+            age = numpy.float32(0.000006)
+        else:
+            age = numpy.float32(0.0)
+
+        self.program.update_particles(
+            self.queue, (self.count,1), None,
+            self.moments.cl_gl_moments,
+            self.lattice.memor