Implement a (virtual) Maxwell’s Demon.
For a prior research project, I implemented a general particle simulation framework. It’s still under active development today, and I will use Maxwell’s Demon as a new test case.
The simulation is written in C++ and allows particle state data and particle interaction laws to be specified via templates (i.e. at compile time). I’ve used Lennard-Jones potentials before; this was my starting point, but I removed the attractive term in order to get a simple “billiard ball” gas model. The particles automatically have essential physical properties like mass, velocity, etc. I ended up adding a globally unique id to each particle to make it easier to implement the demon.
bool GasPotential::interaction(
Scalar& magnitude,
Vector2<Scalar>& force,
Particle<GasParticle>& p0,
Particle<GasParticle>& p1,
PhysicalConstants<GasConstants> const& constants
) {
force = p1.location - p0.location;
Scalar distance = force.squaredNorm();
if (distance == 0) {
return false;
}
if (distance >= constants.cutoff_distance_squared) {
return false;
}
// Compute the magnitude of the interaction force.
distance = std::sqrt(distance);
magnitude = -constants.strength * std::pow(distance, -13);
magnitude = std::min(magnitude, constants.maximum_force);
magnitude = std::max(magnitude, -constants.maximum_force);
// Apply the interaction force.
force *= magnitude / distance;
p0.force += force;
p1.force -= force;
return true;
}
struct GasParticle {
uint32_t uid;
};
With these changes we’re ready to instantiate a model with the appropriate physics.
FunSimModel<GasLaw> model;
This model can be plugged into one of two visualization applications I’ve written: one using Vulkan, and one using GLUT. I’ll use the latter since we’ll be making a few modifications to the rendering and it’s a lot less cumbersome than Vulkan.
So without further ado, let’s build some walls. I already have a pipeline in place for rendering FReps with particles, so I’ll just use this.
// Create a box
auto const outer_box = rectangle(
Vector2<Scalar>(walls_x_low, walls_y_low),
Vector2<Scalar>(walls_x_high, walls_y_high)
);
auto const left_compartment = rectangle(
Vector2<Scalar>(walls_x_low + thickness, walls_y_low + thickness),
Vector2<Scalar>(half_width - 0.9 * thickness, walls_y_high - thickness)
);
auto const right_compartment = rectangle(
Vector2<Scalar>(half_width + 0.9 * thickness, walls_y_low + thickness),
Vector2<Scalar>(walls_x_high - thickness, walls_y_high - thickness)
);
auto const opening = rectangle(
Vector2<Scalar>(half_width - thickness, opening_y_lo),
Vector2<Scalar>(half_width + thickness, opening_y_hi)
);
auto const box_frep = frep_difference(outer_box, left_compartment, right_compartment, opening);
FrepModel box_model;
// These particles have infinite mass (i.e. zero inverse mass) so they can't move.
box_model.add_component(box_frep, 0);
auto box_particles = render_frep<GasLaw>(box_model, lattice_bounds, constants.typical_distance);
for (uint32_t i = 0; i < box_particles.size(); ++i) {
result.particles.add_particle(box_particles.particles()[i]);
result.particles.back().uid = puid++;
}
At this point we have our walls. Note that the walls are made out of particles. I gave them infinite mass, so no force will ever be able to move them — a convenient hack to create fixed structures. The walls aren’t smooth, but since particle interactions already have all the physics we need, this saves me some coding.
Now let’s add some particles.
for (uint32_t i = 0; i < n_particles; ++i) {
Vector2<Scalar> potential_position;
do {
potential_position = {x_distribution(gen), y_distribution(gen)};
} while (!is_valid(potential_position));
particle_positions.emplace_back(potential_position);
particle_velocities.emplace_back(v_distribution(gen), v_distribution(gen));
//particle_velocities.emplace_back(7, 0);
}
for (uint32_t i = 0; i < n_particles; ++i) {
model.particles.add_particle(particle_positions[i], particle_velocities[i]);
model.particles.back().uid = puid++;
}
The particles are spatially distributed according to a uniform distribution, and their x and y
velocity components are drawn from a normal distribution with mean zero. I wrote a helper method
is_valid that makes sure the particles don’t start inside the divider, and are an acceptable
distance away from each other. The latter is necessary since if two particles start very close
together, they could have arbitrarily high potential energy in the initial configuration. Here the
box is filled with 500 particles.
Now we’re ready to program the demon. To do this I’ll use one of the hooks built into my simulation:
the frame_update function. This is an optional user supplied method that is evaluated globally
after each step of the simulation. (In multithreaded simulations it’s more efficient to let each
thread run a more specific function on its own, but since the demon may want to analyze all the
particles, and this isn’t a demanding simulation, I’m fine with the easy slow way.)
I’m too lazy to actually program a moving door, so I’m just going to implement a Maxwell membrane. It will exert a force on particles it wants to bounce back into their respective compartments, and will let others pass through. The code is quick and dirty. The demon pays attention to a limited area of the simulation. When a particle enters this area, it decides if it wants to let it pass through the opening or not. If not, it adds the particle to a running list of blocked particles. (We can’t just re-evaluate from scratch each frame, since the demon will end up slowing down particles before it ends up reversing their direction — we don’t want the demon sucking so much energy from the system!). All blocked particles then have a force applied, calculated as a one dimensional analog to the normal particle interaction force. Any particle that the demon was interacting with that leaves the demon’s domain is removed from the list.
model.frame_update = [=](Domain<GasLaw>& domain) mutable {
// Reset all particle visibility flags and energy variables.
for (auto& pair : interaction_flags) {
pair.second = false;
}
average_energy = (total_left_energy + total_right_energy) / (n_left + n_right);
total_left_energy = 0;
total_right_energy = 0;
n_left = 0;
n_right = 0;
// Iterate through all particles.
for (uint32_t j = 0; j < domain.n_regions(); ++j) {
auto& chunk = domain.region(j).chunk();
for (uint32_t i = 0; i < chunk.size(); ++i) {
auto const& particle = chunk.particle(i);
Scalar const energy = particle.velocity.squaredNorm();
// Update energy averages.
if (chunk.location(i)[0] < half_width) {
++n_left;
total_left_energy += energy;
} else {
++n_right;
total_right_energy += energy;
}
// Check if we should block it.
if (
opening_y_lo <= chunk.location(i)[1] && chunk.location(i)[1] <= opening_y_hi &&
x_maxwell_lo <= chunk.location(i)[0] && chunk.location(i)[0] <= x_maxwell_hi
) {
// Figure out if we need to block this particle.
uint32_t const puid = particle.uid;
bool interact = false;
if (interaction_flags.find(puid) != interaction_flags.end()) {
interaction_flags[puid] = true;
interact = true;
} else {
if (energy >= average_energy && particle.velocity[0] >= 0) {
interaction_flags[puid] = true;
interact = true;
}
if (energy < average_energy && particle.velocity[0] <= 0) {
interaction_flags[puid] = true;
interact = true;
}
}
// Block if necessary.
if (interact) {
Scalar const distance = particle.location[0] - half_width;
Scalar magnitude = force_strength * std::pow(distance, -13);
magnitude = std::min(magnitude, max_force);
magnitude = std::max(magnitude, -max_force);
chunk.force(i)[0] += magnitude;
}
}
}
}
// If we didn't see a particle this frame, remove it from the list.
dead_particles.clear();
for (auto& pair : interaction_flags) {
if (!pair.second) {
dead_particles.push_back(pair.first);
}
}
for (uint32_t puid : dead_particles) {
interaction_flags.erase(puid);
}
// Print total energies occasionally
if (domain.region(0).frame() % 1000 == 0) {
std::cout << "left and right energy sums: " << total_left_energy << ", " << total_right_energy << '\n';
std::cout << "average: " << average_energy << '\n';
}
// Always keep the simulation running.
return true;
};
One final modification: let’s color the particles according to their energies. This is a quick tweak to the GL code in my GLUT app.
if (chunk.inverse_mass(i) != 0) {
Scalar const energy = chunk.velocity(i).squaredNorm();
Scalar const scale = std::min(energy / 30, Scalar(1.0));
glColor3f(scale, 0, (1.0f - scale));
} else {
glColor3f(0, 0, 0);
}
Let’s let the demon do it’s thing! Here it’s sorting faster particles to the left.
In hindsight it seems obvious, but my demon sorts the high energy particles faster than the low energy ones, because the fast ones hit the boundary more frequently. As a result the left side fills up with more particles than the right. In the video above the left side ends up with more than twice as much energy as the right.