Add a fun little fake bonfire _simulation_

…using appropriately colored aging particles

5 files changed, 261 insertions, 8 deletions

diff --git a/channel_2d_gl_interop.py b/channel_2d_gl_interop.py
index 3ca623c..377bbd0 100644
--- a/channel_2d_gl_interop.py
+++ b/channel_2d_gl_interop.py
@@ -171,6 +171,7 @@ projection = get_projection()
 
 particles = Particles(
     lattice.context,
+    lattice.queue,
     lattice.memory.float_type,
     numpy.mgrid[
         lattice.geometry.size_x//20:2*lattice.geometry.size_x//20:100j,
@@ -184,7 +185,7 @@ def on_display():
     lattice.collect_gl_moments()
 
     for i in range(0,updates_per_frame):
-        lattice.update_gl_particles(particles)
+        lattice.update_gl_particles(particles, aging = False)
 
     lattice.sync()
 
diff --git a/simulation.py b/simulation.py
index f8ce4eb..13fa962 100644
--- a/simulation.py
+++ b/simulation.py
@@ -234,8 +234,17 @@ class Lattice:
             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):
+    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)
+            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 1ee39cd..d0d6d30 100644
--- a/template/kernel.mako
+++ b/template/kernel.mako
@@ -141,7 +141,9 @@ __kernel void collect_gl_moments(__global __read_only  ${float_type}* f,
 
 __kernel void update_particles(__global __read_only  float4* moments,
                                __global __read_only  int*    material,
-                               __global __write_only float4* particles)
+                               __global __write_only float4* particles,
+                               __global __read_only  float2* init_particles,
+                               float aging)
 {
   const unsigned int pid = get_global_id(0);
 
@@ -151,12 +153,15 @@ __kernel void update_particles(__global __read_only  float4* moments,
 
   float4 moment = moments[gid];
 
-  if (material[gid] == 1) {
+  if (material[gid] == 1 && particle.z < 1.0) {
     particle.x += moment.y;
     particle.y += moment.z;
+    particle.z += min(particle.x, particle.y) * aging;
   } else {
-    particle.x = particle.z;
-    particle.y = particle.w;
+    float2 orig = init_particles[pid];
+    particle.x = orig.x;
+    particle.y = orig.y;
+    particle.z = particle.z-1.0;
   }
 
   particles[pid] = particle;
diff --git a/trugfeuer_2d_gl_interop.py b/trugfeuer_2d_gl_interop.py
new file mode 100644
index 0000000..c23a355
--- /dev/null
+++ b/trugfeuer_2d_gl_interop.py
@@ -0,0 +1,233 @@
+import numpy
+from string import Template
+
+from simulation         import Lattice, Geometry
+from utility.particles  import Particles
+from symbolic.generator import LBM
+
+import symbolic.D2Q9 as D2Q9
+
+from OpenGL.GL   import *
+from OpenGL.GLUT import *
+
+from OpenGL.GL import shaders
+
+screen_x = 1920
+screen_y = 1200
+pixels_per_cell   = 4
+updates_per_frame = 40
+particle_count    = 100000
+
+inflow = 0.005
+relaxation_time = 0.515
+
+def circle(cx, cy, r):
+    return lambda x, y: (x - cx)**2 + (y - cy)**2 < r*r
+
+def get_channel_material_map(geometry):
+    return [
+        (lambda x, y: x > 0 and x < geometry.size_x-1 and y > 0 and y < geometry.size_y-1, 1), # bulk fluid
+        (lambda x, y: y == 1,                 3), # inflow
+        (lambda x, y: y == geometry.size_y-2, 4), # outflow
+        (lambda x, y: x == 1,                 2), # bottom
+        (lambda x, y: x == geometry.size_x-2, 2), # top
+        (lambda x, y: y > geometry.size_y//20 and y < 2*geometry.size_y//20 and x < 4*geometry.size_x//9, 2),
+        (lambda x, y: y > geometry.size_y//20 and y < 2*geometry.size_y//20 and x > 5*geometry.size_x//9, 2),
+        (circle(geometry.size_x//2   , geometry.size_y//8, 3), 2),
+        (circle(geometry.size_x//2-10, geometry.size_y//8, 3), 2),
+        (circle(geometry.size_x//2+10, geometry.size_y//8, 3), 2),
+        (lambda x, y: x == 0 or x == geometry.size_x-1 or y == 0 or y == geometry.size_y-1, 0) # ghost cells
+    ]
+
+boundary = Template("""
+    if ( m == 2 ) {
+        u_0 = 0.0;
+        u_1 = 0.0;
+    }
+    if ( m == 3 ) {
+        u_0 = 0.0;
+        u_1 = min(time/10000.0 * $inflow, $inflow);
+    }
+    if ( m == 4 ) {
+        rho = 1.0;
+    }
+""").substitute({
+    'inflow': inflow
+})
+
+def get_projection():
+    scale = numpy.diag([(2.0*pixels_per_cell)/screen_x, (2.0*pixels_per_cell)/screen_y, 1.0, 1.0])
+    translation        = numpy.matrix(numpy.identity(4))
+    translation[3,0:3] = [-1.0, -1.0, 0.0]
+    return scale * translation;
+
+def glut_window(fullscreen = False):
+    glutInit(sys.argv)
+    glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH)
+
+    if fullscreen:
+        window = glutEnterGameMode()
+    else:
+        glutInitWindowSize(screen_x, screen_y)
+        glutInitWindowPosition(0, 0)
+        window = glutCreateWindow("LBM")
+
+    return window
+
+lbm = LBM(D2Q9)
+
+window = glut_window(fullscreen = False)
+
+vertex_shader = shaders.compileShader(Template("""
+#version 430
+
+layout (location=0) in vec4 CellMoments;
+
+out vec3 color;
+
+uniform mat4 projection;
+
+vec3 blueRedPalette(float x) {
+    return mix(
+        vec3(0.0, 0.0, 1.0),
+        vec3(1.0, 0.0, 0.0),
+        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.
+    );
+
+    if (CellMoments[3] > 0.0) {
+        //color = blueRedPalette(CellMoments[3] / 0.2);
+        color = vec3(0.0,0.0,0.0);
+    } else {
+        color = vec3(0.4,0.4,0.4);
+    }
+}""").substitute({
+    'size_x': screen_x//pixels_per_cell,
+    'inflow': inflow
+}), GL_VERTEX_SHADER)
+
+fragment_shader = shaders.compileShader("""
+#version 430
+
+in vec3 color;
+
+void main(){
+    gl_FragColor = vec4(color.xyz, 0.0);
+}""", GL_FRAGMENT_SHADER)
+
+particle_shader = shaders.compileShader(Template("""
+#version 430
+
+layout (location=0) in vec4 Particles;
+
+out vec3 color;
+
+uniform mat4 projection;
+
+vec3 fire(float x) {
+    return mix(
+        vec3(1.0, 1.0, 0.0),
+        vec3(1.0, 0.0, 0.0),
+        x
+    );
+}
+
+void main() {
+    gl_Position = projection * vec4(
+        Particles[0],
+        Particles[1],
+        0.,
+        1.
+    );
+
+    color = fire(1.0-Particles[2]);
+}""").substitute({}), GL_VERTEX_SHADER)
+
+moment_program = shaders.compileProgram(vertex_shader, fragment_shader)
+particle_program = shaders.compileProgram(particle_shader, fragment_shader)
+projection_id = shaders.glGetUniformLocation(moment_program, 'projection')
+
+lattice = Lattice(
+    descriptor   = D2Q9,
+    geometry     = Geometry(screen_x//pixels_per_cell, screen_y//pixels_per_cell),
+    moments      = lbm.moments(optimize = False),
+    collide      = lbm.bgk(f_eq = lbm.equilibrium(), tau = relaxation_time),
+    boundary_src = boundary,
+    opengl       = True
+)
+
+lattice.apply_material_map(
+    get_channel_material_map(lattice.geometry))
+lattice.sync_material()
+
+particles = Particles(
+    lattice.context,
+    lattice.queue,
+    lattice.memory.float_type,
+    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,
+    ].reshape(2,-1).T)
+
+projection = get_projection()
+
+def on_display():
+    for i in range(0,updates_per_frame):
+        lattice.evolve()
+
+    lattice.collect_gl_moments()
+
+    for i in range(0,updates_per_frame):
+        lattice.update_gl_particles(particles, 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.asfortranarray(projection))
+    glVertexPointer(4, GL_FLOAT, 0, lattice.memory.gl_moments)
+    glPointSize(pixels_per_cell)
+    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.asfortranarray(projection))
+    glVertexPointer(4, GL_FLOAT, 0, particles.gl_particles)
+    glPointSize(1)
+    glEnable(GL_POINT_SMOOTH)
+    glDrawArrays(GL_POINTS, 0, particles.count)
+
+    glutSwapBuffers()
+
+def on_timer(t):
+    glutTimerFunc(t, on_timer, t)
+    glutPostRedisplay()
+
+glutDisplayFunc(on_display)
+glutTimerFunc(10, on_timer, 10)
+
+glutMainLoop()
diff --git a/utility/particles.py b/utility/particles.py
index b838790..43f32a3 100644
--- a/utility/particles.py
+++ b/utility/particles.py
@@ -7,7 +7,7 @@ import OpenGL.GL as gl
 from OpenGL.arrays import vbo
 
 class Particles:
-    def __init__(self, context, float_type, grid):
+    def __init__(self, context, queue, float_type, grid):
         self.context = context
         self.count = len(grid)
 
@@ -16,6 +16,11 @@ class Particles:
         self.np_particles[:,0:2] = grid
         self.np_particles[:,2:4] = self.np_particles[:,0:2]
 
+        self.np_init_particles = numpy.ndarray(shape=(self.count, 2), dtype=float_type)
+        self.np_init_particles[:,0:2] = grid
+        self.cl_init_particles = cl.Buffer(context, mf.READ_ONLY, size=self.count * 2*numpy.float32(0).nbytes)
+        cl.enqueue_copy(queue, self.cl_init_particles, self.np_init_particles).wait();
+
         self.gl_particles = vbo.VBO(data=self.np_particles, usage=gl.GL_DYNAMIC_DRAW, target=gl.GL_ARRAY_BUFFER)
         self.gl_particles.bind()
         self.cl_gl_particles = cl.GLBuffer(self.context, mf.READ_WRITE, int(self.gl_particles))