aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--gas.py134
1 files changed, 126 insertions, 8 deletions
diff --git a/gas.py b/gas.py
index de047a1..c4c7062 100644
--- a/gas.py
+++ b/gas.py
@@ -6,8 +6,17 @@ 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, ThreadPoolExecutor
+
glutInit()
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH)
glutInitWindowPosition(0, 0)
@@ -80,19 +89,43 @@ decoration_shader = Shader(
uniform = ['projection', 'face_color']
)
+texture_shader = Shader(
+ fragment_src = """
+ #version 130
+
+ uniform sampler2D screenTexture;
+
+ void main() {
+ gl_FragColor = texture(screenTexture, vec2(gl_TexCoord[0].x, 1.0 - gl_TexCoord[0].y));
+ }""",
+ vertex_src = """
+ #version 130
+
+ in vec3 vertex;
+
+ uniform mat4 projection;
+
+ void main() {
+ gl_Position = projection * vec4(vertex, 1.0);
+ gl_TexCoord[0] = gl_MultiTexCoord0;
+ }""",
+ uniform = ['projection']
+)
+
class View:
- def __init__(self, gas, decorations):
+ 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.0
+ 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.0/2, -1.0/2, 0])
+ translation = matrix44.create_from_translation([-1.05, -1.0/2, 0])
self.pixels_per_unit = height / world_height
self.projection = np.matmul(translation, projection)
@@ -106,6 +139,11 @@ class View:
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))
glEnable(GL_POINT_SPRITE)
@@ -147,26 +185,106 @@ class ColoredBox:
glVertex(self.origin[0] + self.extend[1], self.origin[1] + self.extend[1], 0.)
glEnd()
+def get_histogram(velocities):
+ velocities = np.array([np.linalg.norm(v) for v in velocities])
+ maxwellian = stats.maxwell.fit(velocities)
+
+ fig = plt.figure(figsize=(10,10))
+
+ 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()
+ ncols, nrows = fig.canvas.get_width_height()
+ texture = np.frombuffer(buf, dtype=np.uint8).reshape(nrows, ncols, 3)
+ plt.close()
+
+ return (texture, ncols, nrows)
+
+class VelocityHistogram:
+ def __init__(self, gas, origin, extend):
+ self.gas = gas
+ self.origin = origin
+ self.extend = extend
+ self.steps = -1
+
+ self.pool = ThreadPoolExecutor(max_workers=1)
+ self.plotter = None
+
+ def setup(self):
+ self.texture_id = glGenTextures(1)
+ glBindTexture(GL_TEXTURE_2D, self.texture_id)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
+ glEnable(GL_TEXTURE_2D)
+
+ def update(self):
+ self.steps = self.steps + 1
+
+ if self.steps % 100 == 0 and self.plotter == None:
+ self.plotter = self.pool.submit(get_histogram, self.gas.get_velocities())
+
+ else:
+ if self.plotter != None:
+ if self.plotter.done():
+ texture, width, height = self.plotter.result()
+ self.plotter = None
+
+ glEnable(GL_TEXTURE_2D)
+ glBindTexture(GL_TEXTURE_2D, self.texture_id)
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, texture)
+
+
+ def display(self, uniform):
+ glEnable(GL_TEXTURE_2D)
+ glBindTexture(GL_TEXTURE_2D, self.texture_id)
+
+ glBegin(GL_TRIANGLE_STRIP)
+ glTexCoord(0.,0.)
+ glVertex(self.origin[0], self.origin[1] , 0.)
+ glTexCoord(1.,0.)
+ glVertex(self.origin[0] + self.extend[0], self.origin[1] , 0.)
+ glTexCoord(0.,1.)
+ glVertex(self.origin[0] , self.origin[1] + self.extend[1], 0.)
+ glTexCoord(1.,1.)
+ glVertex(self.origin[0] + self.extend[1], self.origin[1] + self.extend[1], 0.)
+ glEnd()
+
grid_width = 30
-radius = 0.002
+radius = 0.003
char_u = 1120
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
+velocity[0,1] = 4*char_u
config = HardSphereSetup(radius, char_u, position, velocity)
gas = GasFlow(config, opengl = True, t_scale = 1.0)
tracer = Tracer(gas, 5)
-view = View(gas, [ColoredBox([0,0], [1,1], (1,1,1)), tracer])
+histogram = VelocityHistogram(gas, [1.1,0], [1,1])
+histogram.setup()
+view = View(gas, [ ColoredBox([0,0], [1,1], (1,1,1)) ], [ histogram ])
def on_display():
- for i in range(0,10):
+ for i in range(0,2):
gas.evolve()
- tracer.update()
+ gas.queue.finish()
+
+ histogram.update()
view.display()