From 908c1dfbd41a1bfe2210e3cfba0e41fda0940371 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Mon, 23 Mar 2020 16:01:56 +0100 Subject: Add some structure to view implementation --- gas.py | 232 +++++++++++++++++++++++++++++++++---------------------------- gasplot.py | 11 ++- kernel.cl | 20 +++--- 3 files changed, 142 insertions(+), 121 deletions(-) diff --git a/gas.py b/gas.py index 78612df..de047a1 100644 --- a/gas.py +++ b/gas.py @@ -8,157 +8,177 @@ import numpy as np from particles import GasFlow, HardSphereSetup, grid_of_random_velocity_particles -def on_timer(t): - glutTimerFunc(t, on_timer, t) - glutPostRedisplay() +glutInit() +glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH) +glutInitWindowPosition(0, 0) +window = glutCreateWindow("gas") -def get_projection(width, height): - world_height = 1.0 - world_width = world_height / height * width - pixels_per_unit = height / world_height +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) - 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.0/2, -1.0/2, 0]) + def use(self): + shaders.glUseProgram(self.program) - return (np.matmul(translation, projection), pixels_per_unit) +particle_shader = Shader( + fragment_src = """ + #version 430 -def on_reshape(width, height): - global projection, pixels_per_unit - glViewport(0,0,width,height) - projection, pixels_per_unit = get_projection(width, height) + in vec3 color; -glutInit() -glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH) -glutInitWindowPosition(0, 0) -window = glutCreateWindow("gas") -glutTimerFunc(10, on_timer, 10) + void main(){ + if (length(gl_PointCoord - vec2(0.5)) > 0.5) { + discard; + } -fragment_shader = shaders.compileShader(""" -#version 430 + gl_FragColor = vec4(color.xyz, 0.0); + }""", + vertex_src = """ + #version 430 -in vec3 color; + layout (location=0) in vec2 particles; -void main(){ - if (length(gl_PointCoord - vec2(0.5)) > 0.5) { - discard; - } + out vec3 color; - gl_FragColor = vec4(color.xyz, 0.0); -}""", GL_FRAGMENT_SHADER) + uniform mat4 projection; -particle_shader = shaders.compileShader(""" -#version 430 + void main() { + gl_Position = projection * vec4(particles, 0., 1.); + color = vec3(0.0); + }""", + uniform = ['projection'] +) -layout (location=0) in vec2 particles; +decoration_shader = Shader( + fragment_src = """ + #version 430 -out vec3 color; + in vec3 color; -uniform mat4 projection; + void main(){ + gl_FragColor = vec4(color.xyz, 0.0); + }""", + vertex_src = """ + #version 430 -void main() { - gl_Position = projection * vec4(particles, 0., 1.); - color = vec3(0.0); -}""", GL_VERTEX_SHADER) + in vec3 vertex; -particle_program = shaders.compileProgram(particle_shader, fragment_shader) + out vec3 color; -background_fragment_shader = shaders.compileShader(""" -#version 430 + uniform mat4 projection; + uniform vec3 face_color; -in vec3 color; + void main() { + gl_Position = projection * vec4(vertex, 1.); + color = face_color; + }""", + uniform = ['projection', 'face_color'] +) -void main(){ - gl_FragColor = vec4(color.xyz, 0.0); -}""", GL_FRAGMENT_SHADER) +class View: + def __init__(self, gas, decorations): + self.gas = gas + self.decorations = decorations -background_vertex_shader = shaders.compileShader(""" -#version 430 + def reshape(self, width, height): + glViewport(0,0,width,height) -in vec3 vertex; + world_height = 1.0 + world_width = world_height / height * width -out vec3 color; + 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.0/2, -1.0/2, 0]) -uniform mat4 projection; -uniform vec3 face_color; + self.pixels_per_unit = height / world_height + self.projection = np.matmul(translation, projection) -void main() { - gl_Position = projection * vec4(vertex, 1.); - color = face_color; -}""", GL_VERTEX_SHADER) -background_program = shaders.compileProgram(background_vertex_shader, background_fragment_shader) + def display(self): + glClearColor(0.4,0.4,0.4,1.) + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) -particle_projection_id = shaders.glGetUniformLocation(particle_program, 'projection') -background_projection_id = shaders.glGetUniformLocation(background_program, 'projection') -background_color_id = shaders.glGetUniformLocation(background_program, 'face_color') + decoration_shader.use() + glUniformMatrix4fv(decoration_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection)) + for decoration in self.decorations: + decoration.display(decoration_shader.uniform) -class View: - def __init__(self, gas): + particle_shader.use() + glUniformMatrix4fv(particle_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection)) + 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.energy = 0 - self.tracer = [ ] + self.iParticle = iParticle + self.trace = [ ] - def update_trace(self): - positions = self.gas.get_positions() - self.tracer.append((positions[5][0],positions[5][1])) - if len(self.tracer) > 100: - self.tracer.pop(0) + def update(self): + position = self.gas.get_positions()[self.iParticle] + self.trace.append((position[0], position[1])) - def draw_trace(self): - glUniform3f(background_color_id, 1., 0., 0.) + def display(self, uniform): + glUniform3f(uniform['face_color'], 1., 0., 0.) glLineWidth(2) glBegin(GL_LINE_STRIP) - for v in self.tracer: + for v in self.trace: glVertex(*v, 0.) glEnd() - def on_display(self): - for i in range(0,10): - self.gas.evolve() - - glClearColor(0.4,0.4,0.4,1.) - glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) +class ColoredBox: + def __init__(self, origin, extend, color): + self.origin = origin + self.extend = extend + self.color = color - shaders.glUseProgram(background_program) - glUniformMatrix4fv(background_projection_id, 1, False, np.asfortranarray(projection)) - glUniform3f(background_color_id, 1., 1., 1.) + def display(self, uniform): + glUniform3f(uniform['face_color'], *self.color) glBegin(GL_TRIANGLE_STRIP) - glVertex(0.,0.,0.) - glVertex(1.,0.,0.) - glVertex(0.,1.,0.) - glVertex(1.,1.,0.) + 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() - if trace: - self.draw_trace() +grid_width = 30 +radius = 0.002 +char_u = 1120 - shaders.glUseProgram(particle_program) - glUniformMatrix4fv(particle_projection_id, 1, False, np.asfortranarray(projection)) - glEnable(GL_POINT_SPRITE) - glPointSize(2*radius*pixels_per_unit) - gas.gl_draw_particles() +position, velocity = grid_of_random_velocity_particles(grid_width, radius, char_u) +velocity[:,:] = 0 +velocity[0,0] = 10*char_u +velocity[0,1] = 1*char_u - glutSwapBuffers() +config = HardSphereSetup(radius, char_u, position, velocity) +gas = GasFlow(config, opengl = True, t_scale = 1.0) - if trace: - self.update_trace() +tracer = Tracer(gas, 5) +view = View(gas, [ColoredBox([0,0], [1,1], (1,1,1)), tracer]) - velocities = gas.get_velocities() - energy = np.sum(np.array([np.linalg.norm(v)**2 for v in velocities])) - if abs(energy - self.energy) > 1e-4: - print("energy = %.05f" % energy) - self.energy = energy +def on_display(): + for i in range(0,10): + gas.evolve() -grid_width = 10 -radius = 0.005 -char_u = 1 -trace = True + tracer.update() -config = HardSphereSetup(radius, char_u, *grid_of_random_velocity_particles(grid_width, radius, char_u)) -gas = GasFlow(config, opengl = True) + view.display() + +def on_reshape(width, height): + view.reshape(width, height) -view = View(gas) +def on_timer(t): + glutTimerFunc(t, on_timer, t) + glutPostRedisplay() -glutDisplayFunc(view.on_display) +glutDisplayFunc(on_display) glutReshapeFunc(on_reshape) +glutTimerFunc(10, on_timer, 10) glutMainLoop() diff --git a/gasplot.py b/gasplot.py index 8adfe4d..44b0a53 100644 --- a/gasplot.py +++ b/gasplot.py @@ -12,7 +12,12 @@ grid_width = 30 radius = 0.002 char_u = 1120 -config = HardSphereSetup(radius, char_u, *grid_of_random_velocity_particles(grid_width, radius, char_u)) +position, velocity = grid_of_random_velocity_particles(grid_width, radius, char_u) +velocity[:,:] = 0 +velocity[0,0] = 10.75*char_u +velocity[0,1] = -.25*char_u +config = HardSphereSetup(radius, char_u, position, velocity) +#config = HardSphereSetup(radius, char_u, *grid_of_random_velocity_particles(grid_width, radius, char_u)) gas = GasFlow(config) m_nitrogen = 0.028 / const.N_A @@ -21,7 +26,7 @@ def plot(step, velocities): velocities = np.array([np.linalg.norm(v) for v in velocities]) maxwellian = stats.maxwell.fit(velocities) - print("T = %.0f K; u_mean = %.0f [m/s]; energy = %.05f" % ((maxwellian[1]**2 / const.k * m_nitrogen, stats.maxwell.mean(*maxwellian), np.sum([x*2 for x in velocities])))) + print("T = %.0f K; u_mean = %.0f [m/s]; energy = %.05f" % ((maxwellian[1]**2 / const.k * m_nitrogen, stats.maxwell.mean(*maxwellian), np.sum([x**2 for x in velocities])))) plt.figure() @@ -52,4 +57,4 @@ def simulate(n_steps, section): plot(i, velocities) -simulate(10000, 500) +simulate(100000, 1000) diff --git a/kernel.cl b/kernel.cl index 8bbe28a..0a8056b 100644 --- a/kernel.cl +++ b/kernel.cl @@ -108,22 +108,18 @@ __kernel void evolve(__global vec_t* pos_a, pos_b[i] = p; vel_b[i] = v; } else { - vec_t p_ = pos_a[jParticle]; - vec_t v_ = vel_a[jParticle]; - - p += min2intersect * v; - p_ += min2intersect * v_; + if (i < jParticle) { + vec_t p_ = pos_a[jParticle]; + vec_t v_ = vel_a[jParticle]; - vec_t omega = normalize(p - p_); + p += min2intersect * v; + p_ += min2intersect * v_; - while (dot(omega,omega) < 1.0) { - omega *= 1.001; - } + vec_t omega = normalize(p - p_); - v -= dot(vel_a[i] - vel_a[jParticle], omega) * omega; - v_ -= dot(vel_a[jParticle] - vel_a[i], omega) * omega; + v -= dot(vel_a[i] - vel_a[jParticle], omega) * omega; + v_ -= dot(vel_a[jParticle] - vel_a[i], omega) * omega; - if (i < jParticle) { p += (DELTA_T - min2intersect) * v; p_ += (DELTA_T - min2intersect) * v_; -- cgit v1.2.3