From 7a79d49dd27b3ec923cc70daadeae6f0e48a7590 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Tue, 8 Oct 2019 20:17:46 +0200 Subject: Play around with 2d streamlines --- channel_2d_streamlines_gl_interop.py | 180 +++++++++++++++++++++++++++++++++++ template/streamline.mako | 59 ++++++++++++ utility/streamline.py | 72 ++++++++++++++ 3 files changed, 311 insertions(+) create mode 100644 channel_2d_streamlines_gl_interop.py create mode 100644 template/streamline.mako create mode 100644 utility/streamline.py diff --git a/channel_2d_streamlines_gl_interop.py b/channel_2d_streamlines_gl_interop.py new file mode 100644 index 0000000..c1aa1e8 --- /dev/null +++ b/channel_2d_streamlines_gl_interop.py @@ -0,0 +1,180 @@ +import numpy +from string import Template + +from simulation import Lattice, Geometry +from utility.opengl import MomentsVertexBuffer +from utility.streamline import Streamlines +from symbolic.generator import LBM + +import symbolic.D2Q9 as D2Q9 + +from OpenGL.GL import * +from OpenGL.GLUT import * + +from OpenGL.GL import shaders + +from pyrr import matrix44 + +lattice_x = 1024 +lattice_y = 256 + +updates_per_frame = 10 + +inflow = 0.01 +relaxation_time = 0.51 + +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: x == 1, 3), # inflow + (lambda x, y: x == geometry.size_x-2, 4), # outflow + (lambda x, y: y == 1, 2), # bottom + (lambda x, y: y == geometry.size_y-2, 2), # top + + (circle(1.0*geometry.size_x//6, 1*geometry.size_y//3, geometry.size_y//5), 2), + (circle(1.5*geometry.size_x//6, 2*geometry.size_y//3, geometry.size_y//6), 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 = min(time/10000.0 * $inflow, $inflow); + u_1 = 0.0; + } + if ( m == 4 ) { + rho = 1.0; + } +""").substitute({ + 'inflow': inflow +}) + +def get_projection(width, height): + world_width = lattice_x + world_height = world_width / width * height + + projection = matrix44.create_orthogonal_projection(-world_width/2, world_width/2, -world_height/2, world_height/2, -1, 1) + translation = matrix44.create_from_translation([-lattice_x/2, -lattice_y/2, 0]) + + point_size = width / world_width + + return numpy.matmul(translation, projection), point_size + +def glut_window(fullscreen = False): + glutInit(sys.argv) + glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH) + + if fullscreen: + window = glutEnterGameMode() + else: + glutInitWindowSize(800, 600) + glutInitWindowPosition(0, 0) + window = glutCreateWindow("LBM") + + return window + +lbm = LBM(D2Q9) + +window = glut_window(fullscreen = False) + +vertex_shader = shaders.compileShader(""" +#version 430 + +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); +} + +void main(){ + const vec2 sample_pos = unit(frag_pos); + result = texture(moments, sample_pos); +} +""").substitute({ + "size_x": lattice_x, + "size_y": lattice_y, + "inflow": inflow +}), GL_FRAGMENT_SHADER) + +shader_program = shaders.compileProgram(vertex_shader, fragment_shader) +projection_id = shaders.glGetUniformLocation(shader_program, 'projection') + +lattice = Lattice( + descriptor = D2Q9, + geometry = Geometry(lattice_x, lattice_y), + 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() + +moments_vbo = MomentsVertexBuffer(lattice) +streamline_texture = Streamlines( + lattice, moments_vbo, + list(map(lambda y: [2, y*lattice.geometry.size_y//48], range(1,48)))) + +def on_display(): + for i in range(0,updates_per_frame): + lattice.evolve() + + moments_vbo.collect() + streamline_texture.update() + + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) + + shaders.glUseProgram(shader_program) + glUniformMatrix4fv(projection_id, 1, False, numpy.asfortranarray(projection)) + streamline_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() + +def on_reshape(width, height): + global projection, point_size + glViewport(0,0,width,height) + projection, point_size = get_projection(width, height) + +def on_timer(t): + glutTimerFunc(t, on_timer, t) + glutPostRedisplay() + +glutDisplayFunc(on_display) +glutReshapeFunc(on_reshape) +glutTimerFunc(10, on_timer, 10) + +glutMainLoop() diff --git a/template/streamline.mako b/template/streamline.mako new file mode 100644 index 0000000..73c80e2 --- /dev/null +++ b/template/streamline.mako @@ -0,0 +1,59 @@ +<% +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) +%> + +__kernel void dillute(__global int* material, + __read_write image2d_t streamlines) +{ + const unsigned int gid = ${gid()}; + const int2 pos = (int2)(get_global_id(0), get_global_id(1)); + + float4 color = read_imagef(streamlines, pos); + + if (material[gid] == 1) { + color.xyz *= 0.975; + } else { + color.xyz = 0.2; + } + + write_imagef(streamlines, pos, color); +} + +float3 blueRedPalette(float x) { + return mix( + (float3)(0.0, 0.0, 1.0), + (float3)(1.0, 0.0, 0.0), + x + ); +} + +__kernel void draw_streamline(__global float4* moments, + __global int* material, + __global float2* origins, + __read_write image2d_t streamlines) +{ + float2 particle = origins[get_global_id(0)]; + + for (int i = 0; i < ${2*memory.size_x}; ++i) { + const unsigned int gid = round(particle.y)*${memory.size_x} + round(particle.x); + const float4 moment = moments[gid]; + + if (material[gid] != 1) { + break; + } + + particle.x += 0.5 * moment.y / 0.01; + particle.y += 0.5 * moment.z / 0.01; + + const int2 pos = (int2)(round(particle.x), round(particle.y)); + + float4 color = read_imagef(streamlines, pos); + color.xyz += 0.05; + + write_imagef(streamlines, pos, color); + } +} diff --git a/utility/streamline.py b/utility/streamline.py new file mode 100644 index 0000000..5111fc2 --- /dev/null +++ b/utility/streamline.py @@ -0,0 +1,72 @@ +import pyopencl as cl +mf = cl.mem_flags + +import numpy + +from mako.template import Template +from pathlib import Path + +from OpenGL.GL import * +from OpenGL.arrays import vbo + +class Streamlines: + def __init__(self, lattice, moments, origins): + 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(origins) + self.np_origins = numpy.ndarray(shape=(self.count, 2), dtype=self.float_type) + + for i, pos in enumerate(origins): + self.np_origins[i,:] = pos + + self.cl_origins = cl.Buffer(self.context, mf.READ_ONLY, size=self.count * 2*numpy.float32(0).nbytes) + cl.enqueue_copy(self.queue, self.cl_origins, self.np_origins).wait(); + + 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_streamlines = glGenTextures(1) + self.gl_texture_type = GL_TEXTURE_2D + glBindTexture(self.gl_texture_type, self.gl_streamlines) + + glTexImage2D(self.gl_texture_type, 0, GL_RGBA32F, 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_LINEAR); + glTexParameteri(self.gl_texture_type, GL_TEXTURE_MAG_FILTER, GL_LINEAR) + 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_streamlines = cl.GLTexture(self.lattice.context, mf.READ_WRITE, self.gl_texture_type, 0, self.gl_streamlines, 2) + + self.build_kernel() + + def build_kernel(self): + program_src = Template(filename = str(Path(__file__).parent/'../template/streamline.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, location = GL_TEXTURE0): + glEnable(self.gl_texture_type) + glActiveTexture(location); + glBindTexture(self.gl_texture_type, self.gl_streamlines) + + def update(self): + cl.enqueue_acquire_gl_objects(self.queue, [self.cl_gl_streamlines]) + + self.program.dillute( + self.queue, (self.lattice.memory.size_x,self.lattice.memory.size_y), None, + self.lattice.memory.cl_material, + self.cl_gl_streamlines) + + self.program.draw_streamline( + self.queue, (self.count,1), None, + self.moments.cl_gl_moments, + self.lattice.memory.cl_material, + self.cl_origins, + self.cl_gl_streamlines) -- cgit v1.2.3