From 235123e65b5edb3d5f0c8f14b05e2d2e56f9e40a Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Tue, 24 Mar 2020 20:51:06 +0100 Subject: Separate classes into modules --- boltzgas/__init__.py | 4 + boltzgas/initial_condition.py | 15 ++ boltzgas/kernel.cl | 156 ++++++++++++++++++++ boltzgas/simulation.py | 131 +++++++++++++++++ boltzgas/visual/__init__.py | 4 + boltzgas/visual/box.py | 16 +++ boltzgas/visual/histogram.py | 125 ++++++++++++++++ boltzgas/visual/shader.py | 15 ++ boltzgas/visual/tracer.py | 20 +++ boltzgas/visual/view.py | 154 ++++++++++++++++++++ gas.py | 325 ++---------------------------------------- kernel.cl | 156 -------------------- particles.py | 142 ------------------ shell.nix | 3 +- 14 files changed, 653 insertions(+), 613 deletions(-) create mode 100644 boltzgas/__init__.py create mode 100644 boltzgas/initial_condition.py create mode 100644 boltzgas/kernel.cl create mode 100644 boltzgas/simulation.py create mode 100644 boltzgas/visual/__init__.py create mode 100644 boltzgas/visual/box.py create mode 100644 boltzgas/visual/histogram.py create mode 100644 boltzgas/visual/shader.py create mode 100644 boltzgas/visual/tracer.py create mode 100644 boltzgas/visual/view.py delete mode 100644 kernel.cl delete mode 100644 particles.py diff --git a/boltzgas/__init__.py b/boltzgas/__init__.py new file mode 100644 index 0000000..3b623dd --- /dev/null +++ b/boltzgas/__init__.py @@ -0,0 +1,4 @@ +import boltzgas.initial_condition +import boltzgas.visual + +from boltzgas.simulation import HardSphereSetup, HardSphereSimulation diff --git a/boltzgas/initial_condition.py b/boltzgas/initial_condition.py new file mode 100644 index 0000000..4d66263 --- /dev/null +++ b/boltzgas/initial_condition.py @@ -0,0 +1,15 @@ +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)) + + 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) diff --git a/boltzgas/kernel.cl b/boltzgas/kernel.cl new file mode 100644 index 0000000..ad984bb --- /dev/null +++ b/boltzgas/kernel.cl @@ -0,0 +1,156 @@ +#define RADIUS $radius +#define N_PARTICLES $n_particles +#define DELTA_T $delta_t + +typedef float scalar_t; +typedef float2 vec_t; + +#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 scalar_t wall_loc[4] = { + -RADIUS, 1.-RADIUS, 1.-RADIUS, -RADIUS +}; + +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; + return (loc - dot(p,n)) / dot(wall_v,n); + } else { + return SCALAR_MAX; + } +} + +scalar_t pos_min_or_infty(scalar_t t0, scalar_t t1) { + if (t0 >= 0.) { + if (t1 >= 0.) { + return min(min(t0, t1), SCALAR_MAX); + } else { + return min(t0, SCALAR_MAX); + } + } else { + if (t1 >= 0.) { + return min(t1, SCALAR_MAX); + } else { + return SCALAR_MAX; + } + } +} + +scalar_t solve_collision(vec_t p, vec_t v, vec_t p_, vec_t v_) { + scalar_t a = dot(v-v_, v-v_); + scalar_t b = 2.*dot(p-p_, v-v_); + scalar_t c = dot(p-p_, p-p_) - 4.*RADIUS*RADIUS; + scalar_t d = b*b - 4.*a*c; + + if (d >= 0.) { + scalar_t t0 = (-b + sqrt(d))/(2.*a); + scalar_t t1 = (-b - sqrt(d))/(2.*a); + + return pos_min_or_infty(t0, t1); + } else { + return SCALAR_MAX; + } +} + +__kernel void evolve(__global vec_t* pos_a, + __global vec_t* vel_a, + __global vec_t* pos_b, + __global 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]; + + 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]; + + scalar_t time2intersect = solve_collision(p, v, p_, v_); + + if (time2intersect < min2intersect) { + min2intersect = time2intersect; + jParticle = iParticle; + } + } + } + + 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); + if (time2wall < min2wall) { + min2wall = time2wall; + jWall = iWall; + } + } + + if (min2intersect < DELTA_T) { + if (min2wall < min2intersect) { + p += min2wall * v; + v -= 2*dot(v,wall_normals[jWall])*wall_normals[jWall]; + p += (DELTA_T - min2wall) * v; + last_collide[i] = N_PARTICLES; + + pos_b[i] = p; + vel_b[i] = v; + } else { + if (i < jParticle) { + vec_t p_ = pos_a[jParticle]; + vec_t v_ = 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; + + 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_; + + last_collide[i] = jParticle; + last_collide[jParticle] = i; + } + } + } else { + if (min2wall < DELTA_T) { + p += min2wall * v; + v -= 2*dot(v,wall_normals[jWall])*wall_normals[jWall]; + p += (DELTA_T - min2wall) * v; + last_collide[i] = N_PARTICLES; + } else { + p += DELTA_T * v; + } + + pos_b[i] = p; + vel_b[i] = v; + } +} + +__kernel void get_velocity_norms(__global vec_t* velocities, __global scalar_t* norms) +{ + unsigned int i = get_global_id(0); + + norms[i] = length(velocities[i]); +} diff --git a/boltzgas/simulation.py b/boltzgas/simulation.py new file mode 100644 index 0000000..a4ae4da --- /dev/null +++ b/boltzgas/simulation.py @@ -0,0 +1,131 @@ +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('boltzgas/kernel.cl', 'r') as kernel_src: + return Template(kernel_src.read()).substitute( + delta_t = delta_t, + n_particles = n_particles, + radius = radius) + +class HardSphereSimulation: + 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) + self.cl_particle_velocity_norms = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_particle_velocity_norms) + + 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.radius = setup.radius + self.char_u = setup.char_u + + 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.np_particle_velocity_norms = np.ndarray((self.n_particles, 1), dtype=np.float32) + + self.kernel_src = build_kernel(self.t_scale*self.radius/self.char_u, self.n_particles, self.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_velocity_norms(self): + if self.tick: + self.program.get_velocity_norms(self.queue, (self.n_particles,), None, self.cl_particle_velocity_b, self.cl_particle_velocity_norms) + cl.enqueue_copy(self.queue, self.np_particle_velocity_norms, self.cl_particle_velocity_norms).wait() + return self.np_particle_velocity_norms + else: + self.program.get_velocity_norms(self.queue, (self.n_particles,), None, self.cl_particle_velocity_a, self.cl_particle_velocity_norms) + cl.enqueue_copy(self.queue, self.np_particle_velocity_norms, self.cl_particle_velocity_norms).wait() + return self.np_particle_velocity_norms + + 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 diff --git a/boltzgas/visual/__init__.py b/boltzgas/visual/__init__.py new file mode 100644 index 0000000..3462d96 --- /dev/null +++ b/boltzgas/visual/__init__.py @@ -0,0 +1,4 @@ +from .box import ColoredBox +from .histogram import VelocityHistogram +from .tracer import Tracer +from .view import View diff --git a/boltzgas/visual/box.py b/boltzgas/visual/box.py new file mode 100644 index 0000000..0a6ab5a --- /dev/null +++ b/boltzgas/visual/box.py @@ -0,0 +1,16 @@ +from OpenGL.GL import * + +class ColoredBox: + def __init__(self, origin, extend, color): + self.origin = origin + self.extend = extend + self.color = color + + def display(self, uniform): + glUniform3f(uniform['face_color'], *self.color) + glBegin(GL_TRIANGLE_STRIP) + glVertex(self.origin[0], self.origin[1] , 0.) + glVertex(self.origin[0] + self.extend[0], self.origin[1] , 0.) + glVertex(self.origin[0] , self.origin[1] + self.extend[1], 0.) + glVertex(self.origin[0] + self.extend[1], self.origin[1] + self.extend[1], 0.) + glEnd() diff --git a/boltzgas/visual/histogram.py b/boltzgas/visual/histogram.py new file mode 100644 index 0000000..2f90b8b --- /dev/null +++ b/boltzgas/visual/histogram.py @@ -0,0 +1,125 @@ +from OpenGL.GL import * + +import numpy as np + +import scipy.stats as stats +import scipy.constants as const +from scipy.optimize import minimize + +import matplotlib +import matplotlib.pyplot as plt + +from concurrent.futures import Future, ProcessPoolExecutor + +def get_histogram(velocities, char_u): + maxwellian = stats.maxwell.fit(velocities) + + plt.style.use('dark_background') + + fig = plt.figure(figsize=(8,8)) + ax = fig.add_axes([0.1, 0.06, 0.88, 0.92]) + + plt.ylim(0, 0.003) + plt.ylabel('Probability') + + plt.xlim(0, 1.2*char_u) + plt.xlabel('Velocity magnitude [m/s]') + + plt.hist(velocities, bins=50, density=True, alpha=0.5, label='Simulated velocities') + + xs = np.linspace(0, 1.2*char_u, 100) + plt.plot(xs, stats.maxwell.pdf(xs, *maxwellian), label='Maxwell-Boltzmann distribution') + + plt.legend(loc='upper right') + + fig.canvas.draw() + buf = fig.canvas.tostring_rgb() + width, height = fig.canvas.get_width_height() + texture = np.frombuffer(buf, dtype=np.uint8).reshape(width, height, 3) + + plt.close(fig) + + return (texture, height, width) + + +class VelocityHistogram: + def __init__(self, gas, origin, extend): + self.gas = gas + self.origin = origin + self.extend = extend + self.steps = -1 + + self.pool = ProcessPoolExecutor(max_workers=1) + self.plotter = None + + self.tick = False + self.mixing = 0.0 + + def setup(self): + self.vertices = np.array([ + self.origin[0] , self.origin[1] , 0., 1., + self.origin[0] + self.extend[0], self.origin[1] , 1., 1., + self.origin[0] , self.origin[1] + self.extend[1], 0., 0., + self.origin[0] + self.extend[0], self.origin[1] + self.extend[1], 1., 0. + ], dtype=np.float32) + + self.vao = glGenVertexArrays(1) + self.vbo = glGenBuffers(1) + + glBindVertexArray(self.vao) + glBindBuffer(GL_ARRAY_BUFFER, self.vbo) + glBufferData(GL_ARRAY_BUFFER, self.vertices.tostring(), GL_STATIC_DRAW) + + glEnableVertexAttribArray(0) + glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 4*np.dtype(np.float32).itemsize, ctypes.c_void_p(0)) + glEnableVertexAttribArray(1) + glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 4*np.dtype(np.float32).itemsize, ctypes.c_void_p(2*np.dtype(np.float32).itemsize)) + + self.texture_id = glGenTextures(2) + + glBindTexture(GL_TEXTURE_2D, self.texture_id[0]) + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + + glBindTexture(GL_TEXTURE_2D, self.texture_id[1]) + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + + def update(self): + self.steps = self.steps + 1 + + if self.steps % 50 == 0 and self.plotter == None: + self.plotter = self.pool.submit(get_histogram, self.gas.get_velocity_norms(), self.gas.char_u) + + else: + if self.plotter != None and self.plotter.done(): + texture, width, height = self.plotter.result() + if self.tick: + glBindTexture(GL_TEXTURE_2D, self.texture_id[0]) + glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, texture) + self.tick = False + self.mixing = 1.0 + + else: + glBindTexture(GL_TEXTURE_2D, self.texture_id[1]) + glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, texture) + self.tick = True + self.mixing = 0.0 + + self.plotter = None + + + def display(self, uniform): + if self.tick: + self.mixing = min(self.mixing+0.1, 1.0); + else: + self.mixing = max(self.mixing-0.1, 0.0); + + glBindTextures(self.texture_id[0], 2, self.texture_id) + glUniform1iv(uniform['picture'], len(self.texture_id), self.texture_id) + glUniform1f(uniform['mixing'], self.mixing) + + glBindVertexArray(self.vao); + glDrawArrays(GL_TRIANGLE_STRIP, 0, 4) + glBindVertexArray(0) + diff --git a/boltzgas/visual/shader.py b/boltzgas/visual/shader.py new file mode 100644 index 0000000..b4295a8 --- /dev/null +++ b/boltzgas/visual/shader.py @@ -0,0 +1,15 @@ +from OpenGL.GL import GL_VERTEX_SHADER, GL_FRAGMENT_SHADER +from OpenGL.GL import shaders + +class Shader: + def __init__(self, fragment_src, vertex_src, uniform): + self.program = shaders.compileProgram( + shaders.compileShader(vertex_src, GL_VERTEX_SHADER), + shaders.compileShader(fragment_src, GL_FRAGMENT_SHADER)) + self.uniform = { } + for name in uniform: + self.uniform[name] = shaders.glGetUniformLocation(self.program, name) + + def use(self): + shaders.glUseProgram(self.program) + diff --git a/boltzgas/visual/tracer.py b/boltzgas/visual/tracer.py new file mode 100644 index 0000000..a5b2203 --- /dev/null +++ b/boltzgas/visual/tracer.py @@ -0,0 +1,20 @@ +from OpenGL.GL import * + +class Tracer: + def __init__(self, gas, iParticle): + self.gas = gas + self.iParticle = iParticle + self.trace = [ ] + + def update(self): + position = self.gas.get_positions()[self.iParticle] + self.trace.append((position[0], position[1])) + + def display(self, uniform): + glUniform3f(uniform['face_color'], 1., 0., 0.) + glLineWidth(2) + glBegin(GL_LINE_STRIP) + for v in self.trace: + glVertex(*v, 0.) + glEnd() + diff --git a/boltzgas/visual/view.py b/boltzgas/visual/view.py new file mode 100644 index 0000000..d6e188b --- /dev/null +++ b/boltzgas/visual/view.py @@ -0,0 +1,154 @@ +from OpenGL.GL import * +from OpenGL.GLUT import * + +from pyrr import matrix44 + +import numpy as np + +from .shader import Shader + +particle_shader = ( + """ + #version 430 + + in vec3 color; + + void main(){ + if (length(gl_PointCoord - vec2(0.5)) > 0.5) { + discard; + } + + gl_FragColor = vec4(color.xyz, 0.0); + } + """, + """ + #version 430 + + layout (location=0) in vec2 particles; + + out vec3 color; + + uniform mat4 projection; + uniform vec3 face_color; + uniform vec3 trace_color; + uniform uint trace_id; + + void main() { + gl_Position = projection * vec4(particles, 0., 1.); + + if (gl_VertexID == trace_id) { + color = trace_color; + } else { + color = face_color; + } + } + """, + ['projection', 'face_color', 'trace_color', 'trace_id'] +) + +decoration_shader = ( + """ + #version 430 + + in vec3 color; + + void main(){ + gl_FragColor = vec4(color.xyz, 0.0); + } + """, + """ + #version 430 + + in vec3 vertex; + + out vec3 color; + + uniform mat4 projection; + uniform vec3 face_color; + + void main() { + gl_Position = projection * vec4(vertex, 1.); + color = face_color; + } + """, + ['projection', 'face_color'] +) + +texture_shader = ( + """ + #version 430 + + in vec2 tex_coord; + + uniform sampler2D picture[2]; + + uniform float mixing; + + void main() { + gl_FragColor = mix(texture(picture[0], tex_coord), texture(picture[1], tex_coord), mixing); + } + """, + """ + #version 430 + + layout (location=0) in vec2 screen_vertex; + layout (location=1) in vec2 texture_vertex; + + out vec2 tex_coord; + + uniform mat4 projection; + + void main() { + gl_Position = projection * vec4(screen_vertex, 0.0, 1.0); + tex_coord = texture_vertex; + } + """, + ['picture','projection','mixing'] +) + +class View: + def __init__(self, gas, decorations, windows): + self.gas = gas + self.decorations = decorations + self.windows = windows + + self.texture_shader = Shader(*texture_shader) + self.particle_shader = Shader(*particle_shader) + self.decoration_shader = Shader(*decoration_shader) + + def reshape(self, width, height): + glViewport(0,0,width,height) + + world_height = 1.4 + world_width = world_height / height * width + + 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]) + + self.pixels_per_unit = height / world_height + self.projection = np.matmul(translation, projection) + + def display(self): + glClearColor(0.,0.,0.,1.) + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) + + self.decoration_shader.use() + glUniformMatrix4fv(self.decoration_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection)) + for decoration in self.decorations: + decoration.display(self.decoration_shader.uniform) + + self.texture_shader.use() + glUniformMatrix4fv(self.texture_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection)) + for window in self.windows: + window.display(self.texture_shader.uniform) + + self.particle_shader.use() + glUniformMatrix4fv(self.particle_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection)) + 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'], 4) + glEnable(GL_POINT_SPRITE) + glPointSize(2*self.gas.radius*self.pixels_per_unit) + self.gas.gl_draw_particles() + + glutSwapBuffers() diff --git a/gas.py b/gas.py index c1ccbce..9261d87 100644 --- a/gas.py +++ b/gas.py @@ -1,319 +1,16 @@ from OpenGL.GL import * from OpenGL.GLUT import * -from OpenGL.GL import shaders - -from pyrr import matrix44 import numpy as np -import scipy.stats as stats -import scipy.constants as const -from scipy.optimize import minimize - -import matplotlib -import matplotlib.pyplot as plt - -from particles import GasFlow, HardSphereSetup, grid_of_random_velocity_particles - -from concurrent.futures import Future, ProcessPoolExecutor +from boltzgas import HardSphereSetup, HardSphereSimulation +from boltzgas.initial_condition import grid_of_random_velocity_particles +from boltzgas.visual import View, VelocityHistogram, Tracer, ColoredBox glutInit() glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH) glutInitWindowPosition(0, 0) -window = glutCreateWindow("gas") - -class Shader: - def __init__(self, vertex_src, fragment_src, uniform): - self.program = shaders.compileProgram( - shaders.compileShader(vertex_src, GL_VERTEX_SHADER), - shaders.compileShader(fragment_src, GL_FRAGMENT_SHADER)) - self.uniform = { } - for name in uniform: - self.uniform[name] = shaders.glGetUniformLocation(self.program, name) - - def use(self): - shaders.glUseProgram(self.program) - -particle_shader = Shader( - fragment_src = """ - #version 430 - - in vec3 color; - - void main(){ - if (length(gl_PointCoord - vec2(0.5)) > 0.5) { - discard; - } - - gl_FragColor = vec4(color.xyz, 0.0); - }""", - vertex_src = """ - #version 430 - - layout (location=0) in vec2 particles; - - out vec3 color; - - uniform mat4 projection; - uniform vec3 face_color; - uniform vec3 trace_color; - uniform uint trace_id; - - void main() { - gl_Position = projection * vec4(particles, 0., 1.); - - if (gl_VertexID == trace_id) { - color = trace_color; - } else { - color = face_color; - } - }""", - uniform = ['projection', 'face_color', 'trace_color', 'trace_id'] -) - -decoration_shader = Shader( - fragment_src = """ - #version 430 - - in vec3 color; - - void main(){ - gl_FragColor = vec4(color.xyz, 0.0); - }""", - vertex_src = """ - #version 430 - - in vec3 vertex; - - out vec3 color; - - uniform mat4 projection; - uniform vec3 face_color; - - void main() { - gl_Position = projection * vec4(vertex, 1.); - color = face_color; - }""", - uniform = ['projection', 'face_color'] -) - -texture_shader = Shader( - fragment_src = """ - #version 430 - - in vec2 tex_coord; - - uniform sampler2D picture[2]; - - uniform float mixing; - - void main() { - gl_FragColor = mix(texture(picture[0], tex_coord), texture(picture[1], tex_coord), mixing); - }""", - vertex_src = """ - #version 430 - - layout (location=0) in vec2 screen_vertex; - layout (location=1) in vec2 texture_vertex; - - out vec2 tex_coord; - - uniform mat4 projection; - - void main() { - gl_Position = projection * vec4(screen_vertex, 0.0, 1.0); - tex_coord = texture_vertex; - }""", - uniform = ['picture','projection','mixing'] -) - -class View: - def __init__(self, gas, decorations, windows): - self.gas = gas - self.decorations = decorations - self.windows = windows - - def reshape(self, width, height): - glViewport(0,0,width,height) - - world_height = 1.4 - world_width = world_height / height * width - - 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]) - - self.pixels_per_unit = height / world_height - self.projection = np.matmul(translation, projection) - - def display(self): - glClearColor(0.,0.,0.,1.) - glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) - - decoration_shader.use() - glUniformMatrix4fv(decoration_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection)) - for decoration in self.decorations: - decoration.display(decoration_shader.uniform) - - texture_shader.use() - glUniformMatrix4fv(texture_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection)) - for window in self.windows: - window.display(texture_shader.uniform) - - particle_shader.use() - glUniformMatrix4fv(particle_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection)) - glUniform3f(particle_shader.uniform['face_color'], 1., 1., 1.) - glUniform3f(particle_shader.uniform['trace_color'], 1., 0., 0.) - glUniform1ui(particle_shader.uniform['trace_id'], 4) - glEnable(GL_POINT_SPRITE) - glPointSize(2*radius*self.pixels_per_unit) - self.gas.gl_draw_particles() - - glutSwapBuffers() - -class Tracer: - def __init__(self, gas, iParticle): - self.gas = gas - self.iParticle = iParticle - self.trace = [ ] - - def update(self): - position = self.gas.get_positions()[self.iParticle] - self.trace.append((position[0], position[1])) - - def display(self, uniform): - glUniform3f(uniform['face_color'], 1., 0., 0.) - glLineWidth(2) - glBegin(GL_LINE_STRIP) - for v in self.trace: - glVertex(*v, 0.) - glEnd() - -class ColoredBox: - def __init__(self, origin, extend, color): - self.origin = origin - self.extend = extend - self.color = color - - def display(self, uniform): - glUniform3f(uniform['face_color'], *self.color) - glBegin(GL_TRIANGLE_STRIP) - glVertex(self.origin[0], self.origin[1] , 0.) - glVertex(self.origin[0] + self.extend[0], self.origin[1] , 0.) - glVertex(self.origin[0] , self.origin[1] + self.extend[1], 0.) - glVertex(self.origin[0] + self.extend[1], self.origin[1] + self.extend[1], 0.) - glEnd() - -def get_histogram(velocities): - maxwellian = stats.maxwell.fit(velocities) - - plt.style.use('dark_background') - - fig = plt.figure(figsize=(10,10)) - ax = fig.add_axes([0.1, 0.05, 0.88, 0.93]) - - plt.ylim(0, 0.003) - plt.ylabel('Probability') - - plt.xlim(0, 1.2*char_u) - plt.xlabel('Velocity magnitude [m/s]') - - plt.hist(velocities, bins=50, density=True, alpha=0.5, label='Simulated velocities') - - xs = np.linspace(0, 1.2*char_u, 100) - plt.plot(xs, stats.maxwell.pdf(xs, *maxwellian), label='Maxwell-Boltzmann distribution') - - plt.legend(loc='upper right') - - fig.canvas.draw() - buf = fig.canvas.tostring_rgb() - width, height = fig.canvas.get_width_height() - texture = np.frombuffer(buf, dtype=np.uint8).reshape(width, height, 3) - - plt.close(fig) - - return (texture, height, width) - - -class VelocityHistogram: - def __init__(self, gas, origin, extend): - self.gas = gas - self.origin = origin - self.extend = extend - self.steps = -1 - - self.pool = ProcessPoolExecutor(max_workers=1) - self.plotter = None - - self.tick = False - self.mixing = 0.0 - - def setup(self): - self.vertices = np.array([ - self.origin[0] , self.origin[1] , 0., 1., - self.origin[0] + self.extend[0], self.origin[1] , 1., 1., - self.origin[0] , self.origin[1] + self.extend[1], 0., 0., - self.origin[0] + self.extend[0], self.origin[1] + self.extend[1], 1., 0. - ], dtype=np.float32) - - self.vao = glGenVertexArrays(1) - self.vbo = glGenBuffers(1) - - glBindVertexArray(self.vao) - glBindBuffer(GL_ARRAY_BUFFER, self.vbo) - glBufferData(GL_ARRAY_BUFFER, self.vertices.tostring(), GL_STATIC_DRAW) - - glEnableVertexAttribArray(0) - glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 4*np.dtype(np.float32).itemsize, ctypes.c_void_p(0)) - glEnableVertexAttribArray(1) - glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 4*np.dtype(np.float32).itemsize, ctypes.c_void_p(2*np.dtype(np.float32).itemsize)) - - self.texture_id = glGenTextures(2) - - glBindTexture(GL_TEXTURE_2D, self.texture_id[0]) - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); - - glBindTexture(GL_TEXTURE_2D, self.texture_id[1]) - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); - - def update(self): - self.steps = self.steps + 1 - - if self.steps % 50 == 0 and self.plotter == None: - self.plotter = self.pool.submit(get_histogram, self.gas.get_velocity_norms()) - - else: - if self.plotter != None and self.plotter.done(): - texture, width, height = self.plotter.result() - if self.tick: - glBindTexture(GL_TEXTURE_2D, self.texture_id[0]) - glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, texture) - self.tick = False - self.mixing = 1.0 - - else: - glBindTexture(GL_TEXTURE_2D, self.texture_id[1]) - glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, texture) - self.tick = True - self.mixing = 0.0 - - self.plotter = None - - - def display(self, uniform): - if self.tick: - self.mixing = min(self.mixing+0.1, 1.0); - else: - self.mixing = max(self.mixing-0.1, 0.0); - - glBindTextures(self.texture_id[0], 2, self.texture_id) - glUniform1iv(uniform['picture'], len(self.texture_id), self.texture_id) - glUniform1f(uniform['mixing'], self.mixing) - - glBindVertexArray(self.vao); - glDrawArrays(GL_TRIANGLE_STRIP, 0, 4) - glBindVertexArray(0) +glutCreateWindow("BoltzGas") grid_width = 32 radius = 0.004 @@ -325,7 +22,7 @@ velocity[0,0] = 10*char_u velocity[0,1] = 4*char_u config = HardSphereSetup(radius, char_u, position, velocity) -gas = GasFlow(config, opengl = True, t_scale = 1.0) +gas = HardSphereSimulation(config, opengl = True, t_scale = 0.1) tracer = Tracer(gas, 4) histogram = VelocityHistogram(gas, [1.1,0], [1,1]) @@ -358,11 +55,11 @@ def on_keyboard(key, x, y): def on_close(): histogram.pool.shutdown() -if __name__ == "__main__": - glutDisplayFunc(on_display) - glutReshapeFunc(on_reshape) - glutTimerFunc(10, on_timer, 10) - glutKeyboardFunc(on_keyboard) - glutCloseFunc(on_close) +glutDisplayFunc(on_display) +glutReshapeFunc(on_reshape) +glutTimerFunc(10, on_timer, 10) +glutKeyboardFunc(on_keyboard) +glutCloseFunc(on_close) +if __name__ == "__main__": glutMainLoop() diff --git a/kernel.cl b/kernel.cl deleted file mode 100644 index ad984bb..0000000 --- a/kernel.cl +++ /dev/null @@ -1,156 +0,0 @@ -#define RADIUS $radius -#define N_PARTICLES $n_particles -#define DELTA_T $delta_t - -typedef float scalar_t; -typedef float2 vec_t; - -#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 scalar_t wall_loc[4] = { - -RADIUS, 1.-RADIUS, 1.-RADIUS, -RADIUS -}; - -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; - return (loc - dot(p,n)) / dot(wall_v,n); - } else { - return SCALAR_MAX; - } -} - -scalar_t pos_min_or_infty(scalar_t t0, scalar_t t1) { - if (t0 >= 0.) { - if (t1 >= 0.) { - return min(min(t0, t1), SCALAR_MAX); - } else { - return min(t0, SCALAR_MAX); - } - } else { - if (t1 >= 0.) { - return min(t1, SCALAR_MAX); - } else { - return SCALAR_MAX; - } - } -} - -scalar_t solve_collision(vec_t p, vec_t v, vec_t p_, vec_t v_) { - scalar_t a = dot(v-v_, v-v_); - scalar_t b = 2.*dot(p-p_, v-v_); - scalar_t c = dot(p-p_, p-p_) - 4.*RADIUS*RADIUS; - scalar_t d = b*b - 4.*a*c; - - if (d >= 0.) { - scalar_t t0 = (-b + sqrt(d))/(2.*a); - scalar_t t1 = (-b - sqrt(d))/(2.*a); - - return pos_min_or_infty(t0, t1); - } else { - return SCALAR_MAX; - } -} - -__kernel void evolve(__global vec_t* pos_a, - __global vec_t* vel_a, - __global vec_t* pos_b, - __global 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]; - - 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]; - - scalar_t time2intersect = solve_collision(p, v, p_, v_); - - if (time2intersect < min2intersect) { - min2intersect = time2intersect; - jParticle = iParticle; - } - } - } - - 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); - if (time2wall < min2wall) { - min2wall = time2wall; - jWall = iWall; - } - } - - if (min2intersect < DELTA_T) { - if (min2wall < min2intersect) { - p += min2wall * v; - v -= 2*dot(v,wall_normals[jWall])*wall_normals[jWall]; - p += (DELTA_T - min2wall) * v; - last_collide[i] = N_PARTICLES; - - pos_b[i] = p; - vel_b[i] = v; - } else { - if (i < jParticle) { - vec_t p_ = pos_a[jParticle]; - vec_t v_ = 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; - - 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_; - - last_collide[i] = jParticle; - last_collide[jParticle] = i; - } - } - } else { - if (min2wall < DELTA_T) { - p += min2wall * v; - v -= 2*dot(v,wall_normals[jWall])*wall_normals[jWall]; - p += (DELTA_T - min2wall) * v; - last_collide[i] = N_PARTICLES; - } else { - p += DELTA_T * v; - } - - pos_b[i] = p; - vel_b[i] = v; - } -} - -__kernel void get_velocity_norms(__global vec_t* velocities, __global scalar_t* norms) -{ - unsigned int i = get_global_id(0); - - norms[i] = length(velocities[i]); -} diff --git a/particles.py b/particles.py deleted file mode 100644 index 048e52c..0000000 --- a/particles.py +++ /dev/null @@ -1,142 +0,0 @@ -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) - self.cl_particle_velocity_norms = cl.Buffer(self.context, mf.COPY_HOST_PTR, hostbuf=self.np_particle_velocity_norms) - - 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.np_particle_velocity_norms = np.ndarray((self.n_particles, 1), dtype=np.float32) - - 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_velocity_norms(self): - if self.tick: - self.program.get_velocity_norms(self.queue, (self.n_particles,), None, self.cl_particle_velocity_b, self.cl_particle_velocity_norms) - cl.enqueue_copy(self.queue, self.np_particle_velocity_norms, self.cl_particle_velocity_norms).wait() - return self.np_particle_velocity_norms - else: - self.program.get_velocity_norms(self.queue, (self.n_particles,), None, self.cl_particle_velocity_a, self.cl_particle_velocity_norms) - cl.enqueue_copy(self.queue, self.np_particle_velocity_norms, self.cl_particle_velocity_norms).wait() - return self.np_particle_velocity_norms - - 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 diff --git a/shell.nix b/shell.nix index d82cbe3..d2ea65e 100644 --- a/shell.nix +++ b/shell.nix @@ -37,7 +37,8 @@ pkgs.stdenvNoCC.mkDerivation rec { shellHook = '' export NIX_SHELL_NAME="${name}" + export PYTHONPATH="$PWD/boltzgas:$PYTHONPATH" export PYOPENCL_COMPILER_OUTPUT=1 - export PYTHONPATH="$PWD:$PYTHONPATH" + export PYTHONDONTWRITEBYTECODE=1 ''; } -- cgit v1.2.3