WEEK01: introduction, mathematical software, parallel programming, benchmarking

[demon image borrowed from SkiFree, a silly Windows game from the early 90s]

This week we got an introduction to the course. Neil went through a number of open-source software tools for symbolic and numerical math work, along with programming, visualization, and 3D graphics. We also touched on benchmarking using an iterative pi calculation program; this is a useful tool for quantitatively comparing the performance of various computing platforms across a range of scales.

Our assignment was to model Maxwell's Demon using tools of our choice. Maxwell's Demon is the creation of James Clerk Maxwell; in the thought experiment, a sealed container is divided into two zones separated by a door that can be opened or closed without expending any energy. The container is filled with gas molecules that have different energy levels (speeds). A demon controls the door, only opening it when fast molecules come near and closing it to prevent slow molecules from passing through. By doing so, the container eventually sorts the gas into hot and cold sections, thereby decreasing the overall system entropy and violating the second law of thermodynamics.

files

• Maxwell's Demon Jupyter Notebook

model

I decided to use Python along with matplotlib and numpy to run the simulation. I did the development within a Jupyter Notebook for ease of rapid iteration and because it can output static HTML documents for posting. The project was a bit of a learning curve -- I'm used to semicolons rather than tabs (i.e., new to Python)-- but eventually I got the simulation to work.

I created a structured array called particles that keeps track of each particle's x/y position, previous x/y position, x/y velocity, and color. These were initialized and filled with random values using the numpy random generator, then color transparency values were set to the velocity vector magnitude. I added a line to represent the demon's door, then created an animation loop that updates particle position and tests for wall and door collisions. Here is that animation playing from initialization through sorting; it's a bit scratchy but you can see the particles sorting themselves out across the door (animation loops every ~25 seconds or so):

The code is fairly simple; figuring out how to test and modify the structured array's values took a bit of time but it worked out:

%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# fixed random seed
np.random.seed(8675309)

# create an empty box for particles
box = plt.figure(figsize=(6, 6))
ax = box.add_axes([0,0,1,1], frameon=False)
ax.set_xlim(0,1), ax.set_xticks([])
ax.set_ylim(0,1), ax.set_yticks([])

# create structured array of particles
n_particles = 500
particles = np.zeros(n_particles, dtype=[('position', float, 2),
                                         ('position_prev', float, 2),
                                         ('speed', float, 2),
                                         ('color', float, 4)])

# initialize particles: random speed, random position; color = sqrt(xspd^2 + yspd^2)
particles['position'] = np.random.uniform(0,1,(n_particles, 2))
particles['speed'] = np.random.uniform(-0.005,0.005, (n_particles, 2))
particles['color'][:, 3] = np.sqrt(np.add(np.square(particles['speed'][:, 0]),
                                np.square(particles['speed'][:, 0]))) * 100

# make scatter plot which will be updated
scat = ax.scatter(particles['position'][:,0], particles['position'][:,1],
                                facecolors=particles['color'])

# add demon gate
gate = ax.plot([.5, .5],[0,1])


def update(frame_number):
    # update particle position based on speed and direction
    particles['position'][:, 0] += particles['speed'][:, 0]
    particles['position'][:, 1] += particles['speed'][:, 1]

    for p in particles:

        # check for wall collisions and bounce
        if p['position'][0] > 1 or p['position'][0] < 0:
            p['speed'][0] *= -1
        if p['position'][1] > 1 or p['position'][1] < 0:
            p['speed'][1] *= -1

        # check for speed at gate and bounce
        if p['position'][0] > 0.5 and p['position_prev'][0] < 0.5 and p['color'][3] < 0.35:
            p['speed'][0] *= -1
        if p['position'][0] < 0.5 and p['position_prev'][0] > 0.5 and p['color'][3] >= 0.35:
            p['speed'][0] *= -1

        # update previous particle positions for gate check
        p['position_prev'] = p['position']

    # update scatter plot
    scat.set_offsets(particles['position'])


# make the animation using the update function as the director
animation = FuncAnimation(box, update, interval=20)
plt.show