aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--boltzgas/initial_condition.py13
-rw-r--r--boltzgas/kernel.cl120
-rw-r--r--boltzgas/simulation.py4
-rw-r--r--boltzgas/visual/view.py90
-rw-r--r--random_velocities.py12
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)