diff options
| -rw-r--r-- | boltzgas/initial_condition.py | 13 | ||||
| -rw-r--r-- | boltzgas/kernel.cl | 120 | ||||
| -rw-r--r-- | boltzgas/simulation.py | 4 | ||||
| -rw-r--r-- | boltzgas/visual/view.py | 90 | ||||
| -rw-r--r-- | random_velocities.py | 12 |
5 files changed, 187 insertions, 52 deletions
diff --git a/boltzgas/initial_condition.py b/boltzgas/initial_condition.py index 4d66263..2526566 100644 --- a/boltzgas/initial_condition.py +++ b/boltzgas/initial_condition.py @@ -1,15 +1,20 @@ 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)) + np_position = np.ndarray((width**3, 4)) + np_velocity = np.ndarray((width**3, 4)) grid = np.meshgrid(np.linspace(2*radius, 1-2*radius, width), + 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_position[:,2] = grid[2].flatten() + np_position[:,3] = 0 - 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,))) + np_velocity[:,0] = u_scale*(-0.5 + np.random.random_sample((width**3,))) + np_velocity[:,1] = u_scale*(-0.5 + np.random.random_sample((width**3,))) + np_velocity[:,2] = u_scale*(-0.5 + np.random.random_sample((width**3,))) + np_velocity[:,3] = 0 return (np_position, np_velocity) diff --git a/boltzgas/kernel.cl b/boltzgas/kernel.cl index ad984bb..aabd1af 100644 --- a/boltzgas/kernel.cl +++ b/boltzgas/kernel.cl @@ -2,22 +2,44 @@ #define N_PARTICLES $n_particles #define DELTA_T $delta_t +#define ENABLE_3D + typedef float scalar_t; + +#ifdef ENABLE_3D +typedef float3 vec_t; +typedef float4 data_vec_t; +#define N_WALLS 6 +#else typedef float2 vec_t; +typedef float2 data_vec_t; +#define N_WALLS 4 +#endif #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 float4 wall_normals[N_WALLS] = { + (float4)(-1., 0., 0., 0.), + (float4)( 1., 0., 0., 0.), + (float4)( 0., 1., 0., 0.), + (float4)( 0.,-1., 0., 0.), + (float4)( 0., 0.,-1., 0.), + (float4)( 0., 0., 1., 0.), }; -__constant scalar_t wall_loc[4] = { - -RADIUS, 1.-RADIUS, 1.-RADIUS, -RADIUS + +__constant scalar_t wall_loc[N_WALLS] = { + -RADIUS, 1.-RADIUS, 1.-RADIUS, -RADIUS, -RADIUS, 1.-RADIUS }; +vec_t wall_normal(unsigned int iWall) { +#ifdef ENABLE_3D + return wall_normals[iWall].xyz; +#else + return wall_normals[iWall].xy; +#endif +} + 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; @@ -59,24 +81,59 @@ scalar_t solve_collision(vec_t p, vec_t v, vec_t p_, vec_t v_) { } } -__kernel void evolve(__global vec_t* pos_a, - __global vec_t* vel_a, - __global vec_t* pos_b, - __global vec_t* vel_b, +vec_t getPosition(__global data_vec_t* pos, + unsigned int iParticle) +{ +#ifdef ENABLE_3D + return pos[iParticle].xyz; +#else + return pos[iParticle]; +#endif +} + +vec_t getVelocity(__global data_vec_t* vel, + unsigned int iParticle) +{ +#ifdef ENABLE_3D + return vel[iParticle].xyz; +#else + return vel[iParticle]; +#endif +} + +void set(__global data_vec_t* pos, + __global data_vec_t* vel, + unsigned int iParticle, + vec_t p, + vec_t v) +{ +#ifdef ENABLE_3D + pos[iParticle].xyz = p; + vel[iParticle].xyz = v; +#else + pos[iParticle] = p; + vel[iParticle] = v; +#endif +} + +__kernel void evolve(__global data_vec_t* pos_a, + __global data_vec_t* vel_a, + __global data_vec_t* pos_b, + __global data_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]; + vec_t p = getPosition(pos_a, i); + vec_t v = getVelocity(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]; + vec_t p_ = getPosition(pos_a, iParticle); + vec_t v_ = getVelocity(vel_a, iParticle); scalar_t time2intersect = solve_collision(p, v, p_, v_); @@ -90,8 +147,8 @@ __kernel void evolve(__global vec_t* pos_a, 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); + for (unsigned int iWall=0; iWall < N_WALLS; ++iWall) { + scalar_t time2wall = solve_wall_collision(wall_normal(iWall), wall_loc[iWall], p, v); if (time2wall < min2wall) { min2wall = time2wall; jWall = iWall; @@ -101,33 +158,29 @@ __kernel void evolve(__global vec_t* pos_a, if (min2intersect < DELTA_T) { if (min2wall < min2intersect) { p += min2wall * v; - v -= 2*dot(v,wall_normals[jWall])*wall_normals[jWall]; + v -= 2*dot(v,wall_normal(jWall)) * wall_normal(jWall); p += (DELTA_T - min2wall) * v; last_collide[i] = N_PARTICLES; - pos_b[i] = p; - vel_b[i] = v; + set(pos_b, vel_b, i, p, v); } else { if (i < jParticle) { - vec_t p_ = pos_a[jParticle]; - vec_t v_ = vel_a[jParticle]; + vec_t p_ = getPosition(pos_a, jParticle); + vec_t v_ = getVelocity(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; + v -= dot(getVelocity(vel_a, i) - getVelocity(vel_a, jParticle), omega) * omega; + v_ -= dot(getVelocity(vel_a, jParticle) - getVelocity(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_; + set(pos_b, vel_b, i, p, v); + set(pos_b, vel_b, jParticle, p_, v_); last_collide[i] = jParticle; last_collide[jParticle] = i; @@ -136,21 +189,20 @@ __kernel void evolve(__global vec_t* pos_a, } else { if (min2wall < DELTA_T) { p += min2wall * v; - v -= 2*dot(v,wall_normals[jWall])*wall_normals[jWall]; + v -= 2*dot(v,wall_normal(jWall)) * wall_normal(jWall); p += (DELTA_T - min2wall) * v; last_collide[i] = N_PARTICLES; } else { p += DELTA_T * v; } - pos_b[i] = p; - vel_b[i] = v; + set(pos_b, vel_b, i, p, v); } } -__kernel void get_velocity_norms(__global vec_t* velocities, __global scalar_t* norms) +__kernel void get_velocity_norms(__global data_vec_t* velocities, __global scalar_t* norms) { unsigned int i = get_global_id(0); - norms[i] = length(velocities[i]); + norms[i] = length(getVelocity(velocities, i)); } diff --git a/boltzgas/simulation.py b/boltzgas/simulation.py index d818bab..905a04f 100644 --- a/boltzgas/simulation.py +++ b/boltzgas/simulation.py @@ -95,10 +95,10 @@ class HardSphereSimulation: if self.tick: self.gl_particle_position_b.bind() - gl.glVertexPointer(2, gl.GL_FLOAT, 0, self.gl_particle_position_b) + gl.glVertexPointer(4, 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.glVertexPointer(4, gl.GL_FLOAT, 0, self.gl_particle_position_a) gl.glDrawArrays(gl.GL_POINTS, 0, self.n_particles) gl.glDisableClientState(gl.GL_VERTEX_ARRAY) diff --git a/boltzgas/visual/view.py b/boltzgas/visual/view.py index e2746e2..327214a 100644 --- a/boltzgas/visual/view.py +++ b/boltzgas/visual/view.py @@ -1,7 +1,7 @@ from OpenGL.GL import * from OpenGL.GLUT import * -from pyrr import matrix44 +from pyrr import matrix44, quaternion import numpy as np @@ -24,26 +24,31 @@ particle_shader = ( """ #version 430 - layout (location=0) in vec2 particles; + layout (location=0) in vec4 particles; out vec3 color; uniform mat4 projection; + uniform mat4 rotation; + uniform vec3 face_color; uniform vec3 trace_color; uniform uint trace_id; void main() { - gl_Position = projection * vec4(particles, 0., 1.); + gl_Position = projection * rotation * vec4(particles.xyz, 1.); - if (gl_VertexID == trace_id) { + if (particles.x < 0.0 || particles.x > 1.0 || + particles.y < 0.0 || particles.y > 1.0 || + particles.z < 0.0 || particles.z > 1.0 + ) { color = trace_color; } else { color = face_color; } } """, - ['projection', 'face_color', 'trace_color', 'trace_id'] + ['projection', 'rotation', 'face_color', 'trace_color', 'trace_id'] ) decoration_shader = ( @@ -106,6 +111,69 @@ texture_shader = ( ['picture','projection','mixing'] ) + +class Projection: + def __init__(self, distance): + self.distance = distance + self.ratio = 4./3. + self.update() + + def update(self): + projection = matrix44.create_perspective_projection(20.0, self.ratio, 0.1, 5000.0) + look = matrix44.create_look_at( + eye = [0, 0, -self.distance], + target = [0, 0, 0], + up = [0,-1, 0]) + + self.matrix = np.matmul(look, projection) + + def update_ratio(self, width, height, update_viewport = True): + if update_viewport: + glViewport(0,0,width,height) + + self.ratio = width/height + self.update() + + def update_distance(self, change): + self.distance += change + self.update() + + def get(self): + return self.matrix + +class Rotation: + def __init__(self, shift, x = np.pi, z = np.pi): + self.matrix = matrix44.create_from_translation(shift), + self.rotation_x = quaternion.Quaternion() + self.update(x,z) + + def shift(self, x, z): + self.matrix = np.matmul( + self.matrix, + matrix44.create_from_translation([x,0,z]) + ) + self.inverse_matrix = np.linalg.inv(self.matrix) + + def update(self, x, z): + rotation_x = quaternion.Quaternion(quaternion.create_from_eulers([x,0,0])) + rotation_z = self.rotation_x.conjugate.cross( + quaternion.Quaternion(quaternion.create_from_eulers([0,0,z]))) + self.rotation_x = self.rotation_x.cross(rotation_x) + + self.matrix = np.matmul( + self.matrix, + matrix44.create_from_quaternion(rotation_z.cross(self.rotation_x)) + ) + self.inverse_matrix = np.linalg.inv(self.matrix) + + def get(self): + return self.matrix + + def get_inverse(self): + return self.inverse_matrix + + + class View: def __init__(self, gas, decorations, windows): self.gas = gas @@ -116,12 +184,17 @@ class View: self.particle_shader = Shader(*particle_shader) self.decoration_shader = Shader(*decoration_shader) + self.projection3d = Projection(distance = 7) + self.rotation3d = Rotation([-1/2, -1/2, -1/2], 5*np.pi/4, np.pi/4) + def reshape(self, width, height): glViewport(0,0,width,height) world_height = 1.4 world_width = world_height / height * width + projection = Projection(10) + 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]) @@ -143,10 +216,15 @@ class View: window.display(self.texture_shader.uniform) self.particle_shader.use() - glUniformMatrix4fv(self.particle_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection)) + #glUniformMatrix4fv(self.particle_shader.uniform['projection'], 1, False, np.asfortranarray(self.projection)) + + glUniformMatrix4fv(self.particle_shader.uniform['projection'], 1, False, np.ascontiguousarray(self.projection3d.get())) + glUniformMatrix4fv(self.particle_shader.uniform['rotation'], 1, False, np.ascontiguousarray(self.rotation3d.get())) + 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'], -1) + glEnable(GL_POINT_SPRITE) glPointSize(2*self.gas.radius*self.pixels_per_unit) self.gas.gl_draw_particles() diff --git a/random_velocities.py b/random_velocities.py index cf7a503..2d60f3b 100644 --- a/random_velocities.py +++ b/random_velocities.py @@ -4,20 +4,20 @@ from boltzgas import HardSphereSetup, HardSphereSimulation from boltzgas.initial_condition import grid_of_random_velocity_particles from boltzgas.visual import VelocityHistogram, ColoredBox, Tracer -grid_width = 50 -radius = 0.002 +grid_width = 10 +radius = 0.005 char_u = 1120 position, velocity = grid_of_random_velocity_particles(grid_width, radius, char_u) config = HardSphereSetup(radius, char_u, position, velocity) -gas = HardSphereSimulation(config, opengl = True) +gas = HardSphereSimulation(config, opengl = True, t_scale=1) -tracer = Tracer(gas, int((grid_width**2)/2+grid_width/2)) +#tracer = Tracer(gas, int((grid_width**2)/2+grid_width/2)) histogram = VelocityHistogram(gas, [1.1,0], [1,1]) -decorations = [ ColoredBox([0,0], [1,1], (0.2,0.2,0.2)), tracer ] -instruments = [ histogram, tracer ] +decorations = [ ] +instruments = [ histogram ] windows = [ histogram ] boltzgas.visualizer.simulate(config, gas, instruments, decorations, windows) |