If Feynman Were Teaching Today… A Simplified Python Simulation of Diffusion
Understanding the real world is not always easy. A Python simulation can help. And let's find ways of making it efficient, too (Part 1)
How often has this happened to you? You think you understand a topic well. Then you realise that you hadn't quite understood it as well as you thought.
I've been listening to the Feynman lectures in physics recently. Some of you may be familiar with these lectures. Richard Feynman was a brilliant physicist and an exceptional communicator. In the early 1960s, he was asked to deliver a two-year undergraduate physics course. This would have been "beneath him" as one of the most renowned physicists of his time. But Feynman loved teaching, and he was remarkably good at it. So he accepted on the condition he only delivers the course once!
Luckily for us, the lectures were recorded (audio-only, unfortunately).
What's this got to do with The Python Coding Stack? Some of you may know that I'm a physicist by training, and I worked as a research scientist for many years. As someone with a PhD in physics, surely there's nothing to learn from listening to undergraduate-level lectures, right?
Wrong.
Even in the early lectures when Feynman was covering basic concepts, his descriptions and explanations led me to understand topics I "knew well" at a deeper level. Feynman's talent as a communicator was to help students visualise what's going on intuitively.
If we could transpose Richard Feynman to the present day, he would explain many of the concepts using Python simulations!
Recently, I've been helping physics school teachers communicate ideas in physics and demonstrate key concepts using Python. Today's step-by-step tutorial is inspired by this work, and by one of Feynman's lectures.
But this is a Python publication, not a physics one. So, this article will focus on the coding, including looking at clever ways to optimise the simulation. Don't worry if you're not into physics. This article is mostly about Python!
I'll also simplify the physics concepts in this simulation. This is also a technique Feynman often used. When the priority is understanding the main concepts, cutting some corners is acceptable.
Here's the animation you'll create in this article. Notice how the particles collide with each other and with the sides:
I'll use one of my favourite modules in this tutorial, the turtle module. However, I plan to follow this tutorial with Part 2 in a few weeks, which will look at a different solution using NumPy and Matplotlib.
Particles Moving Around, Crashing Into Each Other
Let's start with a gas in a closed container. Let's simplify by saying that there are particles moving around the container, without worrying about the physical chemistry. And yes, let's make them perfect spheres, too.
These particles are moving around in random directions, and each one has a particular speed. In physics, we use the term velocity to represent the speed and direction of an object.
Each particle keeps moving with the same velocity—the same speed and direction—until it hits another particle or the sides of the container. We'll also assume all the particles are identical.
And one more simplification: the gas and its container exist in two dimensions. Who needs a third dimension?
Right, enough of the physics for now–it's time to start coding.
I'll assume you're familiar with the basics of object-oriented programming, but if you're not, you can read The Python Coding Book, which has a chapter introducing OOP, or the seven-part Harry Potter-themed series on this topic here on The Python Coding Stack: A Magical Tour Through Object-Oriented Programming in Python • Hogwarts School of Codecraft and Algorithmancy
A few more seconds of your time before you carry on with the tutorial. I'm making a number of free passes available for readers of The Stack to access the learning platform at The Python Coding Place. The passes give access to all the video courses for one week and a live session, similar to the ones I have weekly with members of The Python Coding Place. Here's the link to get your one week's pass:
The Container • The Simulation Area
Let's represent the container as a 2D simulation area. You can create a class to represent this area and create the screen as an attribute in this class. Let's define the class in particles_simulation.py:
When you create an instance of this class, you'll need to pass the width and height in pixels, and you can optionally pass a colour if you don't want a white background.
The width and height are used to calculate the x- and y-coordinates representing the four sides of this 2D container. You'll need these to be integers later as you'll use them as arguments in random.randint(), so you can use the floor division operator // to ensure the values are integers.
The last four lines in the .__init__() special method create the screen using the turtle module. They set the size of this screen using the width and height. And .bgcolor() is the method in the turtle module that sets the background colour.
The .tracer(0) method ensures that changes to the state of turtle objects are not drawn on the screen instantly. Instead, you'll need to call .update(), which is a screen method in the turtle module. You can improve the interface by making this a SimulationArea method so that someone using this class doesn't need to access the turtle method directly:
You also create a .set_title() method to provide access to the turtle module's method to update the window's title, and .run_mainloop() to call turtle.mainloop(), which is needed to run the event loop in the turtle module. If you're familiar with some of my other turtle tutorials, I often use the alias turtle.done() instead of turtle.mainloop(). They're the same function. The .run_mainloop() method doesn't need access to the data within the instance. So, it's a static method instead of an instance method.
Now, let's start working on a second script to run the demo: run_demo.py:
When you run this script, Python creates the window that represents the 2D container:
Now, we need to put some particles in the container and get them to move around.
The Particles
We're making some assumptions to simplify this demo. One of these is to assume that the molecules in the gas are spheres. Let's just call them particles.
You can create a Particle class. Let's start by listing some of the attributes a particle should have:
Position: This can be a list containing the x- and y-coordinates.
Velocity: This can be another list containing the x- and y-components of the particle's velocity. Each component can be a negative or positive number to account for direction.
Colour and size: These attributes enable you to customise the way each particle is displayed on the screen.
You'll also need to link the particle with the container so that the particle can access information about the simulation area. Here's the class with the first part of the .__init__() method. We'll work on the position and velocity later:
The class's .__init__() method includes several arguments. I used a type hint for simulation_area to clarify this is an instance of the SimulationArea class. Another advantage of using type hints for .simulation_area is that tools such as IDEs can assist with autocomplete and other aids.
The Turtle() object is the primary tool in the turtle module to create graphics. You could define Particle to inherit from turtle.Turtle so that a Particle object is also a turtle.Turtle object. However, I chose to create a .particle attribute containing a turtle.Turtle() object instead. This is a design choice you can make in situations like this.
The rest of the .__init__() method so far changes the shape and colour of the sprite representing the Turtle object, "raises the pen up" so that the Turtle object doesn't draw a line when it moves, and changes the size of the sprite. The Turtle method .shapesize() accepts a multiplicative factor as an argument. Therefore, the default value of size_factor = 0.5 displays the sprite at half its default size. By default, a sprite in the turtle module occupies a 20 x 20 pixel area.
You can check that this code works in run_demo.py. Note that the import line now includes Particle:
You also need to add simulation_area.update() to ensure the screen's display is updated. There's a particle in the middle of the container now:
Setting the particles' positions and velocities
You can pass arguments for position and velocity when you create a Particle instance. However, these arguments are optional as they have a default value of None. If you need to create lots of particles, you want to create each particle in a random position in the container and assign a random velocity within a desired range.
Let's fix this by adding a few more lines to the .__init__() method. You also need to import the random module:
Let's review the changes in this version:
You import the
randommodule, which you use in the.__init__()method.You define a class attribute,
.MAX_VELOCITY. This attribute is not specific to each instance of the class but is shared by all instances.If
positionis the default valueNone, you assign a list containing two random values. The first value represents the x-position and is a random number between the left and right edges of the screen. The second value is the y-position, and it's chosen using the top and bottom edges as limits. The arguments inrandom.randint()need to be integers. This is why you use the floor division//when defining the edges of the container from the width and height values inSimulationArea.If
velocityis the default valueNone, you assign another list containing two random values. These represent the x- and y-components of the particle's velocity. Sincerandom.randint()only deals with integers, you can create values with one digit after the decimal point using a bit of arithmetic. You multiply the minimum and maximum values you use as arguments inrandom.randint()by 10, but then you divide the answer by 10 again.
Let's create 100 particles in run_demo.py:
This code creates a list with 100 Particle objects. A Turtle object is created in the centre of the screen. Even though all the Particle objects have a .position attribute with x- and y-coordinates, you're not moving the Turtle objects to their new positions yet. We'll deal with this in the next section.
Moving the particles
There are two key steps to start the animation and move all the particles. One of them is to create an animation loop in run_demo.py. But first, let's define Particle.move() in particles_simulation.py:
A particle's speed tells us how much the position changes in a specific time unit. Often, the unit is m/s in scientific applications. However, speed can be defined using any measure of distance and time. To keep this animation simple, you'll use pixels per frame as a measure of a particle's speed.
Therefore, each time the particle moves in the animation, its x-position changes by an amount equal to its velocity's x-component. And the y-position changes by the velocity's y-component.
The final line in .move() calls .setposition(), which is a method in the turtle module to place the Turtle object in its new position. Recall that .particle is the data attribute that represents the Turtle object.
Let's add the animation loop in run_demo.py:
You call .move() for each particle in the list of 100 particles. The display is only updated after the for loop to ensure all particle positions are updated on the screen at the same time. Recall that you set .tracer(0) when setting up SimulationArea to prevent the sprites from being updated before you call .update().
Here's the animation so far. Note that the speed of the animation may vary depending on your operating system and setup. You can adjust the class attribute .MAX_VELOCITY in the Particle class to control the speed of your animation.
You can see particles starting at random positions on the screen, moving in different directions at different speeds.
This is great, but we still need to do more work to get a basic animation running:
Account for particles hitting the edges of the container/screen
Account for particles hitting each other
Let's work on the edges of the container first.
Adding hard sides to the container
Let's get back to the .move() method in Particle. You can check whether the particle has hit any of the four edges. If the particle hits the left or right edge, you can change the sign of the velocity's x-component. If the velocity's x-component is negative, the particle is moving towards the left. After it hits the side, it's moving towards the right, so the x-component is positive. The opposite happens on the other side.
And when the ball hits the top or bottom, you can swap the sign of the velocity's y-component. Let's add this to .move():
Note that you shouldn't check for equality using == when checking whether the particle has hit any of the four sides since the particle's position may never be exactly equal to that side's coordinate.
For example, if the x-velocity is 3 pixels per frame, the particle's x-position changes by 3 pixels every time it moves in the animation. The particle may "jump over" the container side's x-coordinate when it moves. Instead, use < for the negative values (left and bottom) and > for the positive values (right and top) to check whether a particle has hit one of the container's sides.
And the nice thing about creating classes is that you don't need to make any changes to run_demo.py to verify whether this code works since you changed the definition of Particle.move(). Here's the animation now:
The particles stay in the container by bouncing off the sides.
Dealing with collisions between particles
But particles don't just collide with the container's sides. They also collide with each other. We're assuming the particles are spheres of the same mass. Here's another simplification: when particles collide, they always collide head-on. This is not the case in real life, but it's a reasonable simplification for a demonstration designed to help understand how gas particles move in a container.
Let's start with a method to check whether two particles have collided. You need to check whether the distance between the centres of two particles is less than a certain threshold. You can't check whether this distance is exactly zero, as was the case when checking whether a particle has hit the sides. To make the demo visually realistic, you can set the threshold to represent the centre-to-centre distance when the surface of one particle touches the surface of another.
By default, a sprite in the turtle module occupies a 20 x 20 pixel footprint. However, you increase or decrease this size when creating the Particle instance. You can access this size factor using the Turtle class's .shapesize() method. You used this method to set the size of the sprite in Particle.__init__(). However, you can call this method with no arguments since it returns a value representing the shape's size. The .shapesize() method returns a tuple containing the stretch factors along the x- and y-directions. The tuple also includes a third element, which I'll ignore here as it's not relevant.
Since you called .shapesize() with one argument when you created the Turtle object, you applied the same stretch along both x- and y-directions. Therefore, you can use either the first or the second element in the tuple returned by .shapesize() since they're equal:
The hard-coded value of 20 represents the default sprite size in the turtle module, and it's multiplied by the scaling factor you chose when creating the particle (or 0.5 if using the default, as we are doing in this example). This method calls .distance(), which is a method in the turtle library that returns the distance between two turtles. Recall that a Turtle object's position is the centre of the dot displayed. So, .distance() returns the centre-to-centre distance between two particles. The .check_collision() method returns a Boolean to indicate whether the particles have collided.
Now, you can define .collide(), which determines what happens when two particles collide. The assumption that particles all collide head-on will simplify this method. Those of you who want to create a more precise simluation can account for the correct geometry in this method. However, this won't make a huge difference to this demonstration:
When two particles of equal mass collide head-on (assuming an elastic collision), they simply swap their velocities. Note: in future code snippets in this tutorial, I'll break this line into multiple lines to keep the overall width of the code narrower, but I left it as a single line on this first occasion as it's more readable.
Next, you need to check all the particles to find out which ones are colliding during each frame. In this part of the tutorial, we'll use the brute force approach. However, I'll look at a more efficient option later in this article. And Part 2, which I'll publish at a later date, will also examine different solutions.
Let's update run_demo.py:
For each particle, you first move it, and then you check all the other particles in the container to look for collisions. This works. There's a "but" coming—several of them. Before we deal with the potential pitfalls, let's look at the animation now:
If you look closely at the animation, you'll notice some particles clumping together. There are a few issues we need to consider about the algorithm dealing with the collisions. Let's look at this in the next section.
Refining the collision code • Particle collides with itself
First of all, let's consider a particle hitting itself. Eh?! That can't happen, right?
Right, but it currently happens in the code. Let's say you just moved particle A. The code in the while loop in run_demo.py loops through all the particles, including A itself. The distance between A and itself is zero, of course!
Now, this doesn't really matter since the particle is swapping velocities with itself when it collides with itself—this sentence is so weird out of context! So, no harm done. But it's tidier and a bit more efficient if we avoid this.
One option is to loop through an iterable that contains all the particles except particle A. In the inner for loop in the animation loop, instead of iterating through particles again, you could create a generator expression (part for part in particles if part is not particle). It's best not to create a copy of particles with one particle removed, as that option uses up more memory.
As much as I love to use generators, in this case, I'll opt for a simpler and more readable option, which you can implement directly in Particle.collide():
When .collide() is called with the instance as its own argument so that self and other_particle are the same object, the method returns early.
Refining the collision code • Particles colliding twice per frame
The main issue with our algorithm is that if A collides with B, say, the animation loop may later find that B collides with A again, and the collision happens twice in the same frame. If this happens, you swap the velocities between A and B and then swap them back again in the same iteration of the game loop.
This double collision doesn't occur in most collisions in the animation since when A hits B, their positions change. So, it's possible that when it's B's turn, it's no longer close enough to A. However, this double collision does occur in some cases. This artefact gives rise to the clumping of particles you see in the animation.
One way of dealing with this is to keep track of which particles are in collision with each other within a frame. You can add a new attribute to Particles to keep a record of these particles:
Next, you can update .collide() to ensure particles haven't already collided in the frame. You also need to add self to other_particle.in_collision_with:
Note: the highlights in the code blocks show new or modified lines compared to the previous step, but they don't show lines if only their indentation has changed.
Since you iterate through the particles in the animation loop, you'll always process one of the colliding particles before the other. Therefore, you only need to keep one record of the collision. If we consider particles A and B again, and let's say A comes first in the animation loop because it's at a lower index in the list of particles. You only need to store A in B's .in_collision_with list. There's no need to also store B in A's.
Let's go through the process step-by-step:
Each particle has a list
.in_collision_with.When you're processing A and find that it collides with B, you execute the collision code, and you add A to B's list of collisions.
When it's B's turn to be processed, you check A's list of collisions. B will be in that list. Therefore, you skip the collision in this case.
There's one problem. At some point, you need to remove a particle from the collision list of another particle to make sure you can deal with future collisions of the same pair of particles. Let's start by defining a method in Particle to provide an interface for removing particles:
And you can call this method in an elif clause in the animation loop:
If the particles aren't colliding in a frame, but one is on the collision list of the other, then you can remove the particle from the collision list. The collision must have occurred in the previous frame.
Here's the animation you get at this stage. Pick one particle and follow its path:
This is good enough for now. You may still get the odd glitch due to some edge cases, but this happens much less frequently now, if at all!
We'll explore another issue with this version later in this article.
Diffusion
Remember the motivation for writing this simulation? I've been assisting physics school teachers with Python-based exercises to visualise physics concepts. Now that we have a simple model of gas molecules moving around in a 2D container, let's assume that we start the animation with one type of gas in the left half of the container and a different type of gas in the right half. I'll stick to the assumption that all particles are of equal mass, which is unlikely in the real-case scenario, but it doesn't affect the demo.
You can assign different colours to the particles depending on whether they're in the left or right half of the container. First, you can define two helper methods in Particle:
The first of these methods, .assign_colour(), provides an interface to change the colour directly through the Particle class. Therefore, anyone using this class doesn't need to use the turtle module directly. The second method you define, .get_position(), is what some languages call a getter method. It provides access to the values of a data attribute without the need to access the data attribute directly. This step discourages a user from accessing .position and accidentally changing its value. There are other steps you can take to prevent .position from being modified by mistake, but I won't go down that path in this tutorial.
Now, you can replace the list comprehension that creates the 100 particles in run_demo.py to assign different colours to particles in the left and right halves of the container:
The middle of the screen corresponds to x = 0. Therefore, you colour those particles with an x-position less than zero in one colour and the rest in another colour. The hex colours are the main brand colours for The Python Coding Stack (and The Python Coding Place)!
Here's the animation. Notice how all the green particles are in the left half and all the orange ones on the right at the start of the animation. However, the gasses gradually diffuse so that after some time, there's an equal number of green and orange particles in both halves:
The random motion and collisions lead to an even distribution of the two gasses within the container.
Here's the full code so far for reference. First, particles_simulation.py, which contains the two classes:
And next, run_demo.py:
Breathe! This is a long article, and there's plenty going on. If you're looking for a good place to stop for the day and return to this article another day, this is it.
Or you can proceed directly to the second half of the article, which looks at a clever optimisation technique to deal with the slowing down you'll see as you add more particles.
…and if you're pausing here, don't forget to come back later!
More Particles?
Great. Let's add more particles! I just noticed that the number of particles, currently 100, is hard-coded in the current version. Let's create a variable and assign a larger number to it:
Here's the problem. This video is recorded using the same settings as all the previous videos:
It's slow. Painfully slow. Life's too short to wait for this animation to reach a state when sufficient diffusion has occurred.
The problem is that our nested for loops are not efficient. You need to loop through all the particles to move each one. You can't quite avoid that! But the current algorithm loops through each particle again for every particle. Let's go back to particles labelled with letters, and let's start with particle A. After you move particle A, you iterate through all the particles A, B, C, D, E, and so on to check whether A has collided with them.
Then, it's B's turn to move. And once more, you iterate through the same full list of particles A, B, C, D, E, and so on to check for collisions. And you repeat this process for every particle. The computations needed increase rapidly as you add more particles, slowing the animation. Those familiar with Big-O notation may recognise this as O(n^2) but don't worry if you're not.
However, a lot of this work is unnecessary. Take particle A, wherever it is in the container. At any point in time, it's too far away from most of the other particles to be able to collide with them. So, it's wasteful to check all the particles for particle A. Let's find a better way.
Checking for Collisions Only in a Particle's Neighbourhood
What if we divide the whole simulation area into a grid? Then, we figure out which cell each particle is in. Finally, we only check for collisions within each cell separately instead of over the whole area. This avoids checking whether a particle in one corner of the container is about to collide with a particle at the other end!
This approach splits the simulation area into many smaller ones. Since there are fewer particles in each cell, checking for collisions is quicker using this approach.
But don't take my word for it. Let's make changes to SimulationArea and Particles to implement this technique.
Dividing the simulation area into a grid
Let's start with SimulationArea, and add a new parameter to its .__init__() method to define the size of the grid. The efficiency gains will depend on how many cells you have in the grid. You can also work out the width and height of each cell and create a dictionary to keep track of which particles are in which cell:
The data attribute .cells is a dictionary. The keys are tuples of indices representing each cell in the grid. We'll use i and j to represent the indices or coordinates of the cells and x and y to represent the pixel coordinates of the particles. In this animation, the tuple (0, 0) represents the cell in the bottom left corner, and so on. The value associated with each key is an empty list. You'll add particles to these lists later.
The default grid size defined in the signature is 1 x 1, which is equivalent to having no grid at all. But you'll be testing this code with grids with more than one cell. So, in run_demo.py, you can update the call to SimulationArea() to set the grid size:
Since the simulation area is 1000 x 700 pixels, you can choose a grid of 10 x 7 for now to keep things simple. This means that each cell is 100 x 100 pixels.
Assigning particles to cells
Next, you need to assign each particle to one of these cells. You can start by defining a .cell data attribute in Particle.__init__(). Initially, you can assign None to this attribute:
The code needs to associate each particle with a cell based on the particle's position. You need to perform this operation at the start of the animation when you create each particle in a random position on the screen. However, you also need to reassign the particles to a new cell each time they move from one cell to the next as they're moving around the container.
Let's focus on the horizontal dimension first. You can find the difference between a particle's x-position and the left edge of the container and divide by the width of a cell. For example, let's consider a particle with an x-coordinate of -424. The simulation area is 1000 pixels wide. Therefore, the left edge is at -500 since the centre of the container is 0. The difference between these values is 76, which means that the particle is 76 pixels from the left edge.
If the simulation area is split into a 10 x 7 grid, the width of each cell is 100 pixels. Therefore, the particle is in the first column of cells from the left since these cells include any particles that are within 100 pixels of the left edge. The index i is equal to 0 for all cells in the first column. Since 76 divided by 100 is 0.76, you can find the index by taking the floor of this value, 76 // 100.
A new method, .assign_to_cell(), can deal with this calculation:
The local variables cell_i and cell_j contain the i and j indices that represent the cell. If you're using a 10 x 7 grid, there are ten columns of cells with indices 0 to 9. However, particles can briefly have an x-coordinate that's further left from the left edge. For example, if a particle is at x = -499 and it's moving at 2 pixels per frame, it will move to x = -501 in the next frame. You'll recall that you reverse the horizontal velocity of these particles so that they bounce back into the container on the following iteration. However, the calculation in .assign_to_cell() assigns this particle to a non-existant cell with index -1. You need to correct this and for a similar issue on the remaining three sides of the container:
Note that the right and top side corrections are inclusive, whereas the left and bottom ones are exclusive. A particle with x = 0 goes in the cell with i = 0, but a particle with x = 500 goes in the cell with i = 9 even though 1000 / 100 = 10.
You determined the i and j indices corresponding to the particle's cell. If this is the first time you're assigning a particle to a cell, which occurs when you create the particle at the beginning of the animation, you need to perform two steps:
Store the i and j indices in the particle's
.cellattribute. You set this attribute toNonewhen you created it.Place the particle in the correct list in
simulation_area.cells. This is the dictionary containing all the particles in each cell. The tuple with the cell indices is the key in this dictionary.
Here are these two steps added to .assign_to_cell():
However, during the rest of the animation, you need to check whether the particle has drifted into a different cell. Here are the scenarios:
If the particle is still in the same cell, do nothing.
If the particle is no longer in the same cell:
Remove it from the current cell in
simulation_area.cellsAssign the new indices to the particles
.cellattributePlace the particle in the new cell in
simulation_area.cells
Let's implement this logic by adding an elif clause to the if statement from the previous step:
You need to call this method in two places within the class:
When you create each particle
When you move each particle
Let's conclude this section by adding these two calls at the end of .__init__() and .move():
Verifying that particles are assigned correctly
Let's add a few temporary lines in run_demo.py. You can pick one particle as a sample. Let's pick the first one in the list. Each time this particle moves, you can display the indices of the cell it belongs to on the screen.
Here are the changes you need. You can reduce the number of particles to 50 in this section to ensure you can spot the odd one out:
You need to clear the writing from the previous iteration and write the new indices in each iteration. Here's the animation showing the cell indices for one of the particles. Look for the particle with a label next to it:
Even though the grid is not shown, you can see the particle moving from one cell to another when the indices change. Recall that you created a 10 x 7 grid. Therefore, the i indices range from 0 to 9 and the j indices from 0 to 6.
Before moving on to improving the main collision algorithm, let's work on another visualisation. You can pick a cell and change the colour of any particle in that cell. This will show us particles entering and exiting cells by changing their colours.
For this exercise, you can use more particles and a smaller grid size. The smaller grid size leads to larger cells, and more particles ensure that there are always several particles in the cells.
You can add a few more temporary lines in the for loop within the while loop:
The particles in cell (1, 1) will light up in a different colour to the rest. The grid is 4 x 3 in this example. Here's the animation:
The particles that enter cell (1, 1) turn purple, and then they turn grey again when they exit the cell.
You can now remove or comment out the temporary code you wrote in this section.
Running the Animation Within the Individual Cells
The final step is to change the iteration in the while loop to check for collisions only within a cell in the grid instead of across the whole simulation area.
You can loop through the values of the dictionary simulation_area.cells, and then you can loop through each particle in each cell. Before making the change in the animation loop, let's create a helper function in SimulationArea to facilitate accessing the values in simulation_area.cells:
This method ensures that anyone using the SimulationArea doesn't need to dive too deeply into the class's implementation. Note that I'm not including docstrings in this example since it's part of a detailed tutorial. However, these methods normally have docstrings to assist users.
Finally, let's update the animation loop in run_demo.py. You can also go back to using 500 particles, which led to a very slow animation earlier in this tutorial. You should increase the number of cells by returning to a grid size of 10 x 7:
A reminder that the highlighted lines are the ones with changes in the code. However, the rest of the lines in the while loop have also changed as their level of indentation increased.
You divide the animation into several smaller units within each cell. Particles can still move from one cell to another, of course. However, the time-consuming collision checks are significantly reduced.
Here's the final animation:
Compare this with the attempt to run the basic version of the simulation with 500 particles! This simulation is much smoother. In the appendix, I'll compare the performance of the two methods and the effect of different grid sizes.
You can see the gradual diffusion of the two gasses leading to a homogenous distribution after some time. The random movement of particles and the collisions between them are all we need to explain diffusion. Feynman explains many physical processes using this type of visualisation—but using diagrams and imagination instead of Python demonstrations!
The final versions of particles_simulation.py and run_demo.py are available below, after the conclusion.
Final Words
A reminder: I'm simplifying physics concepts and cutting some corners since this is a demonstration rather than a more rigorous real-world simulation. For example, if two particles are sufficiently close to each other to collide but are on opposite sides of the boundary between two cells, they won't collide in the code. You could place particles close enough to a boundary in two cells, but you'll need to tidy up other parts of the code for this since the code currently assumes that a particle can only be in one cell at any one time.
This article is already long—very long. So, I'll conclude the main article here. However, I can't resist adding a bit more in two appendices. Firstly, I'll compare the efficiency gains of this method by comparing different grid sizes—and none at all! Secondly, I'll add a metric to monitor the diffusion more quantitatively.
And if that's not enough, I'll follow this article with Part 2 in a few weeks, exploring yet another technique using Python tools other than turtle.
Code in this article uses Python 3.12
Final Version of Code
# particles_simulation.py
import random
import turtle
class SimulationArea:
def __init__(self, width, height, colour="white", cell_grid=(1, 1)):
self.width = width
self.height = height
# Define the edges of the simulation area as integers
# (we'll need them as integers later)
self.left_edge = -width // 2
self.right_edge = width // 2
self.top_edge = height // 2
self.bottom_edge = -height // 2
# Create a turtle screen object
self.window = turtle.Screen()
self.window.setup(width, height)
self.window.tracer(0)
self.window.bgcolor(colour)
# Define the cell grid
self.cell_grid = cell_grid
self.cell_width = width / cell_grid[0]
self.cell_height = height / cell_grid[1]
self.cells = {
(i, j): []
for i in range(cell_grid[0])
for j in range(cell_grid[1])
}
def update(self):
self.window.update()
def set_title(self, title):
self.window.title(title)
@staticmethod
def run_mainloop():
turtle.mainloop()
def get_sequence_of_cells(self):
return self.cells.values()
class Particle:
MAX_VELOCITY = 4
def __init__(
self,
simulation_area: SimulationArea,
position=None,
velocity=None,
colour="black",
size_factor=0.5,
):
self.simulation_area = simulation_area
self.particle = turtle.Turtle()
self.particle.shape("circle")
self.particle.color(colour)
self.particle.penup()
self.particle.shapesize(size_factor)
if position is None:
position = [
random.randint(
self.simulation_area.left_edge,
self.simulation_area.right_edge,
),
random.randint(
self.simulation_area.bottom_edge,
self.simulation_area.top_edge,
),
]
self.position = position
if velocity is None:
velocity = [
random.randint(
-self.MAX_VELOCITY * 10, self.MAX_VELOCITY * 10
) / 10,
random.randint(
-self.MAX_VELOCITY * 10, self.MAX_VELOCITY * 10
) / 10,
]
self.velocity = velocity
self.in_collision_with = []
self.cell = None
self.assign_to_cell()
def assign_to_cell(self):
# Convert the x and y coordinates to cell indices i and j
cell_i = int(
(self.position[0] - self.simulation_area.left_edge)
// self.simulation_area.cell_width
)
cell_j = int(
(self.position[1] - self.simulation_area.bottom_edge)
// self.simulation_area.cell_height
)
# Ensure the cell indices are within the grid
if cell_i < 0:
cell_i = 0
elif cell_i >= self.simulation_area.cell_grid[0]:
cell_i = self.simulation_area.cell_grid[0] - 1
if cell_j < 0:
cell_j = 0
elif cell_j >= self.simulation_area.cell_grid[1]:
cell_j = self.simulation_area.cell_grid[1] - 1
# Add the particle if it's not already in the
# cell (first iteration)
if not self.cell:
self.cell = (cell_i, cell_j)
self.simulation_area.cells[self.cell].append(self)
elif self.cell != (cell_i, cell_j):
# Remove the particle from the current cell
self.simulation_area.cells[self.cell].remove(self)
# Update the cell assignment
self.cell = (cell_i, cell_j)
# Add the particle to the new cell
self.simulation_area.cells[self.cell].append(self)
def move(self):
self.position[0] += self.velocity[0]
self.position[1] += self.velocity[1]
# Bounce off the edges
if (
self.position[0] < self.simulation_area.left_edge
or self.position[0] > self.simulation_area.right_edge
):
self.velocity[0] *= -1
if (
self.position[1] < self.simulation_area.bottom_edge
or self.position[1] > self.simulation_area.top_edge
):
self.velocity[1] *= -1
self.particle.setposition(self.position)
self.assign_to_cell()
def check_collision(self, other_particle):
return (
self.particle.distance(other_particle.particle)
< 20 * self.particle.shapesize()[0]
)
def collide(self, other_particle):
"""
This assumes that the particles are of equal mass and
that the collision is a head-on elastic collision.
"""
if other_particle is self:
return
if self not in other_particle.in_collision_with:
self.velocity, other_particle.velocity = (
other_particle.velocity,
self.velocity,
)
# Keep a record of collision in other_particle.
# No need to keep it in self
other_particle.in_collision_with.append(self)
def remove_collision(self, other_particle):
self.in_collision_with.remove(other_particle)
def assign_colour(self, colour):
self.particle.color(colour)
def get_position(self):
return self.position # run_demo.py
from collections import deque
from particles_simulation import SimulationArea, Particle
n_particles = 500
simulation_area = SimulationArea(1000, 700, cell_grid=(10, 7))
simulation_area.set_title("Diffusion Demo")
particles = []
for _ in range(n_particles):
particle = Particle(simulation_area)
if particle.get_position()[0] < 0:
particle.assign_colour("#1A6B72")
else:
particle.assign_colour("#FDB33B")
particles.append(particle)
# Animation loop
while True:
for cell in simulation_area.get_sequence_of_cells():
for particle in cell:
particle.move()
for other_particle in cell:
if particle.check_collision(other_particle):
particle.collide(other_particle)
elif other_particle in particle.in_collision_with:
particle.remove_collision(other_particle)
simulation_area.update()Appendix 1 • Quantify Speed Improvement When Using Grids
You can see the difference between the animation speeds when you use the grid approach. But let's time the execution under different grid conditions. To make the comparisons, you can set a seed in the random module to ensure all the starting conditions are identical:
You can use any integer as an argument in random.seed().
You can also comment the .update() call out in the animation loop to prevent the frames from being displayed. This allows you to time the operations in the while loop without including the time it takes to update the screen. Note that as a result of this change, you won't see any output when you run this code.
Finally, you can run the animation for a fixed number of iterations and time the execution time:
The cell_grid argument in SimulationArea is (1, 1), which means that you perform the animation using a single simulation area. This is equivalent to the first version of the code you created earlier in this tutorial before you implemented the grid approach.
You also replace the while loop with a for loop to time how long it takes the simulation to run for 500 iterations. The execution time will vary depending on your computing setup. I'll present the values I got on my system here:
Time elapsed with grid (1, 1): 87.53 sNext, you can use a 4 x 3 grid, which has twelve cells:
Here's the output:
Time elapsed with grid (4, 3): 8.61 sAnd finally, a 10 x 7 grid, which has seventy cells:
The 500 iterations run much faster now:
Time elapsed with grid (10, 7): 2.87 sHere's a table comparing the output for several grid sizes:
Recall that these times represent the code when it doesn't display each frame on the screen. However, when the loop execution time is small, the time it takes to display the animation becomes a bottleneck in the full animation. For example, the 10 x 7 grid animation takes over 12 seconds to run when you include .update() in the animation loop to draw each frame, even though the same animation clocked under 3 seconds when you don't display the frames.
Appendix 2 • Quantify Diffusion
At the start of the animation, all the green particles are in the left half, and all the orange particles are in the right half. At the end of the animation, both types of particles are spread across the whole container. You can quantify this using the mean x-position of all the green and orange particles. Let's make some changes to the code to implement this:
You run this simulation with 200 particles and a grid size of 10 x 7. You also add a data attribute .colour directly in the for loop when creating the particles to keep a record of the type of particle you create. Another variable keeps a record of how many particles there are of each colour.
You store the mean x-coordinates in each frame in a deque instead of a list. You set the deque's maximum length to 500 to keep only the values from the most recent 500 frames. The display on the window's title bar shows the average of these last 500 frames. This moving average minimises noise in the value displayed.
Here's the animation showing the average x-positions approach 0 as the simulation progresses. The particles’ .MAX_VELOCITY has been changed to 4 for this animation:
Stop Stack
#66
I'm making a number of free passes available for readers of The Stack to access the learning platform at The Python Coding Place. The passes give access to all the video courses for one week and a live session, similar to the ones I have weekly with members of The Python Coding Place. Here's the link to get your one week's pass:
Thank you to all those who supported me with a one-off donation recently. This means a lot and helps me focus on writing more articles and keeping more of these articles free for everyone. Here's the link again for anyone who wants to make a one-off donation to support The Python Coding Stack
The Python Coding Book is available (Ebook and paperback). This is the First Edition, which follows from the "Zeroth" Edition that has been available online for a while—Just ask Google for "python book"!
And if you read the book already, I'd appreciate a review on Amazon. These things matter so much for individual authors!
I'm also releasing The NumPy Mindset at the moment. Currently, this is available as an Early Release—I'm publishing chapters as and when they're ready. Members of The Python Coding Place already have access to this Early Release. Everyone else can get it here—You'll get the final ebook version too once it's ready if you get the Early Release version, of course!
And for those who want to join The Python Coding Place to access all of my video courses—past and future—join regular live sessions, and interact with me and other learners on the members-only forum, here's the link:
Any questions? Just ask…
Appendix: Code Blocks
Code Block #1
# particles_simulation.py
import turtle
class SimulationArea:
def __init__(self, width, height, colour="white"):
self.width = width
self.height = height
# Define the edges of the simulation area as integers
# (we'll need them as integers later)
self.left_edge = -width // 2
self.right_edge = width // 2
self.top_edge = height // 2
self.bottom_edge = -height // 2
# Create a turtle screen object
self.window = turtle.Screen()
self.window.setup(width, height)
self.window.tracer(0)
self.window.bgcolor(colour)
Code Block #2
# particles_simulation.py
import turtle
class SimulationArea:
def __init__(self, width, height, colour="white"):
# ...
def update(self):
self.window.update()
def set_title(self, title):
self.window.title(title)
@staticmethod
def run_mainloop():
turtle.mainloop()
Code Block #3
# run_demo.py
from particles_simulation import SimulationArea
simulation_area = SimulationArea(1000, 700)
simulation_area.set_title("Diffusion Demo")
simulation_area.run_mainloop()
Code Block #4
# particles_simulation.py
import turtle
class SimulationArea:
# ...
class Particle:
def __init__(
self,
simulation_area: SimulationArea,
position=None,
velocity=None,
colour="black",
size_factor=0.5,
):
self.simulation_area = simulation_area
self.particle = turtle.Turtle()
self.particle.shape("circle")
self.particle.color(colour)
self.particle.penup()
self.particle.shapesize(size_factor)
Code Block #5
# run_demo.py
from particles_simulation import SimulationArea, Particle
simulation_area = SimulationArea(1000, 700)
simulation_area.set_title("Diffusion Demo")
test_particle = Particle(simulation_area)
simulation_area.update()
simulation_area.run_mainloop()
Code Block #6
# particles_simulation.py
import random
import turtle
class SimulationArea:
# ...
class Particle:
MAX_VELOCITY = 2
def __init__(
self,
simulation_area: SimulationArea,
position=None,
velocity=None,
colour="black",
size_factor=0.5,
):
# ... <Rest of `__init__()` not shown>
# Set a random position if not provided
if position is None:
position = [
random.randint(
self.simulation_area.left_edge,
self.simulation_area.right_edge,
),
random.randint(
self.simulation_area.bottom_edge,
self.simulation_area.top_edge,
),
]
self.position = position
# Set a random velocity if not provided
if velocity is None:
velocity = [
random.randint(
-self.MAX_VELOCITY * 10, self.MAX_VELOCITY * 10
) / 10,
random.randint(
-self.MAX_VELOCITY * 10, self.MAX_VELOCITY * 10
) / 10,
]
self.velocity = velocity
Code Block #7
# run_demo.py
from particles_simulation import SimulationArea, Particle
simulation_area = SimulationArea(1000, 700)
simulation_area.set_title("Diffusion Demo")
particles = [Particle(simulation_area) for _ in range(100)]
simulation_area.update()
simulation_area.run_mainloop()
Code Block #8
# particles_simulation.py
import random
import turtle
class SimulationArea:
# ...
class Particle:
# ...
def move(self):
self.position[0] += self.velocity[0]
self.position[1] += self.velocity[1]
self.particle.setposition(self.position)
Code Block #9
# run_demo.py
from particles_simulation import SimulationArea, Particle
simulation_area = SimulationArea(1000, 700)
simulation_area.set_title("Diffusion Demo")
particles = [Particle(simulation_area) for _ in range(100)]
## You can now remove these two lines since there's a `while` loop
# simulation_area.update()
# simulation_area.run_mainloop()
# Animation loop
while True:
for particle in particles:
particle.move()
simulation_area.update()
Code Block #10
# particles_simulation.py
import random
import turtle
class SimulationArea:
# ...
class Particle:
# ...
def move(self):
self.position[0] += self.velocity[0]
self.position[1] += self.velocity[1]
# Bounce off the edges
if (
self.position[0] < self.simulation_area.left_edge
or self.position[0] > self.simulation_area.right_edge
):
self.velocity[0] *= -1
if (
self.position[1] < self.simulation_area.bottom_edge
or self.position[1] > self.simulation_area.top_edge
):
self.velocity[1] *= -1
self.particle.setposition(self.position)
Code Block #11
# particles_simulation.py
import random
import turtle
class SimulationArea:
# ...
class Particle:
# ...
def check_collision(self, other_particle):
return (
self.particle.distance(other_particle.particle)
< 20 * self.particle.shapesize()[0]
)
Code Block #12
# particles_simulation.py
import random
import turtle
class SimulationArea:
# ...
class Particle:
# ...
def collide(self, other_particle):
"""
This assumes that the particles are of equal mass and
that the collision is a head-on elastic collision.
"""
self.velocity, other_particle.velocity = other_particle.velocity, self.velocity
Code Block #13
# run_demo.py
from particles_simulation import SimulationArea, Particle
simulation_area = SimulationArea(1000, 700)
simulation_area.set_title("Diffusion Demo")
particles = [Particle(simulation_area) for _ in range(100)]
# Animation loop
while True:
for particle in particles:
particle.move()
for other_particle in particles:
if particle.check_collision(other_particle):
particle.collide(other_particle)
simulation_area.update()
Code Block #14
# particles_simulation.py
import random
import turtle
class SimulationArea:
# ...
class Particle:
# ...
def collide(self, other_particle):
"""
This assumes that the particles are of equal mass and
that the collision is a head-on elastic collision.
"""
if other_particle is self:
return
self.velocity, other_particle.velocity = (
other_particle.velocity,
self.velocity,
)
Code Block #15
# particles_simulation.py
# ...
class Particle:
MAX_VELOCITY = 2
def __init__(...): # Parameters not displayed
# ...
self.in_collision_with = []
# ...
Code Block #16
# particles_simulation.py
# ...
class Particle:
# ...
def collide(self, other_particle):
"""
This assumes that the particles are of equal mass and
that the collision is a head-on elastic collision.
"""
if other_particle is self:
return
if self not in other_particle.in_collision_with:
self.velocity, other_particle.velocity = (
other_particle.velocity,
self.velocity,
)
# Keep a record of collision in other_particle.
# No need to keep it in self
other_particle.in_collision_with.append(self)
Code Block #17
# particles_simulation.py
# ...
class Particle:
# ...
def remove_collision(self, other_particle):
self.in_collision_with.remove(other_particle)
Code Block #18
# run_demo.py
# ...
# Animation loop
while True:
for particle in particles:
particle.move()
for other_particle in particles:
if particle.check_collision(other_particle):
particle.collide(other_particle)
elif other_particle in particle.in_collision_with:
particle.remove_collision(other_particle)
simulation_area.update()
Code Block #19
# particles_simulation.py
# ...
class Particle:
# ...
def assign_colour(self, colour):
self.particle.color(colour)
def get_position(self):
return self.position
Code Block #20
# run_demo.py
from particles_simulation import SimulationArea, Particle
simulation_area = SimulationArea(1000, 700)
simulation_area.set_title("Diffusion Demo")
particles = []
for _ in range(100):
particle = Particle(simulation_area)
if particle.get_position()[0] < 0:
particle.assign_colour("#1A6B72")
else:
particle.assign_colour("#FDB33B")
particles.append(particle)
# Animation loop
# ...
Code Block #21
# particles_simulation.py
import random
import turtle
class SimulationArea:
def __init__(self, width, height, colour="white"):
self.width = width
self.height = height
# Define the edges of the simulation area as integers
# (we'll need them as integers later)
self.left_edge = -width // 2
self.right_edge = width // 2
self.top_edge = height // 2
self.bottom_edge = -height // 2
# Create a turtle screen object
self.window = turtle.Screen()
self.window.setup(width, height)
self.window.tracer(0)
self.window.bgcolor(colour)
def update(self):
self.window.update()
def set_title(self, title):
self.window.title(title)
@staticmethod
def run_mainloop():
turtle.mainloop()
class Particle:
MAX_VELOCITY = 2
def __init__(
self,
simulation_area: SimulationArea,
position=None,
velocity=None,
colour="black",
size_factor=0.5,
):
self.simulation_area = simulation_area
self.particle = turtle.Turtle()
self.particle.shape("circle")
self.particle.color(colour)
self.particle.penup()
self.particle.shapesize(size_factor)
if position is None:
position = [
random.randint(
self.simulation_area.left_edge,
self.simulation_area.right_edge,
),
random.randint(
self.simulation_area.bottom_edge,
self.simulation_area.top_edge,
),
]
self.position = position
if velocity is None:
velocity = [
random.randint(
-self.MAX_VELOCITY * 10, self.MAX_VELOCITY * 10
) / 10,
random.randint(
-self.MAX_VELOCITY * 10, self.MAX_VELOCITY * 10
) / 10,
]
self.velocity = velocity
self.in_collision_with = []
def move(self):
self.position[0] += self.velocity[0]
self.position[1] += self.velocity[1]
# Bounce off the edges
if (
self.position[0] < self.simulation_area.left_edge
or self.position[0] > self.simulation_area.right_edge
):
self.velocity[0] *= -1
if (
self.position[1] < self.simulation_area.bottom_edge
or self.position[1] > self.simulation_area.top_edge
):
self.velocity[1] *= -1
self.particle.setposition(self.position)
def check_collision(self, other_particle):
return (
self.particle.distance(other_particle.particle)
< 20 * self.particle.shapesize()[0]
)
def collide(self, other_particle):
"""
This assumes that the particles are of equal mass and
that the collision is a head-on elastic collision.
"""
if other_particle is self:
return
if self not in other_particle.in_collision_with:
self.velocity, other_particle.velocity = (
other_particle.velocity,
self.velocity,
)
# Keep a record of collision in other_particle.
# No need to keep it in self
other_particle.in_collision_with.append(self)
def remove_collision(self, other_particle):
self.in_collision_with.remove(other_particle)
def assign_colour(self, colour):
self.particle.color(colour)
def get_position(self):
return self.position
Code Block #22
# run_demo.py
from particles_simulation import SimulationArea, Particle
simulation_area = SimulationArea(1000, 700)
simulation_area.set_title("Diffusion Demo")
particles = []
for _ in range(100):
particle = Particle(simulation_area)
if particle.get_position()[0] < 0:
particle.assign_colour("#1A6B72")
else:
particle.assign_colour("#FDB33B")
particles.append(particle)
# Animation loop
while True:
for particle in particles:
particle.move()
for other_particle in particles:
if particle.check_collision(other_particle):
particle.collide(other_particle)
elif other_particle in particle.in_collision_with:
particle.remove_collision(other_particle)
simulation_area.update()
Code Block #23
# run_demo.py
from particles_simulation import SimulationArea, Particle
n_particles = 500
simulation_area = SimulationArea(1000, 700)
simulation_area.set_title("Diffusion Demo")
particles = []
for _ in range(n_particles):
particle = Particle(simulation_area)
if particle.get_position()[0] < 0:
particle.assign_colour("#1A6B72")
else:
particle.assign_colour("#FDB33B")
particles.append(particle)
# Animation loop
# ...
Code Block #24
# particles_simulation.py
import random
import turtle
class SimulationArea:
def __init__(self, width, height, colour="white", cell_grid=(1, 1)):
self.width = width
self.height = height
# Define the edges of the simulation area as integers
# (we'll need them as integers later)
self.left_edge = -width // 2
self.right_edge = width // 2
self.top_edge = height // 2
self.bottom_edge = -height // 2
# Create a turtle screen object
self.window = turtle.Screen()
self.window.setup(width, height)
self.window.tracer(0)
self.window.bgcolor(colour)
# Define the cell grid
self.cell_grid = cell_grid
self.cell_width = width / cell_grid[0]
self.cell_height = height / cell_grid[1]
self.cells = {
(i, j): []
for i in range(cell_grid[0])
for j in range(cell_grid[1])
}
# ...
Code Block #25
# run_demo.py
from particles_simulation import SimulationArea, Particle
n_particles = 500
simulation_area = SimulationArea(1000, 700, cell_grid=(10, 7))
simulation_area.set_title("Diffusion Demo")
# ...
Code Block #26
# particles_simulation.py
# ...
class Particle:
MAX_VELOCITY = 2
def __init__(...):
# ...
self.cell = None
Code Block #27
# particles_simulation.py
# ...
class Particle:
# ...
def assign_to_cell(self):
# Convert the x and y coordinates to cell indices i and j
cell_i = int(
(self.position[0] - self.simulation_area.left_edge)
// self.simulation_area.cell_width
)
cell_j = int(
(self.position[1] - self.simulation_area.bottom_edge)
// self.simulation_area.cell_height
)
# ...
Code Block #28
# particles_simulation.py
# ...
class Particle:
# ...
def assign_to_cell(self):
# Convert the x and y coordinates to cell indices i and j
cell_i = int(
(self.position[0] - self.simulation_area.left_edge)
// self.simulation_area.cell_width
)
cell_j = int(
(self.position[1] - self.simulation_area.bottom_edge)
// self.simulation_area.cell_height
)
# Ensure the cell indices are within the grid
if cell_i < 0:
cell_i = 0
elif cell_i >= self.simulation_area.cell_grid[0]:
cell_i = self.simulation_area.cell_grid[0] - 1
if cell_j < 0:
cell_j = 0
elif cell_j >= self.simulation_area.cell_grid[1]:
cell_j = self.simulation_area.cell_grid[1] - 1
# ...
Code Block #29
# particles_simulation.py
# ...
class Particle:
# ...
def assign_to_cell(self):
# Convert the x and y coordinates to cell indices i and j
cell_i = int(
(self.position[0] - self.simulation_area.left_edge)
// self.simulation_area.cell_width
)
cell_j = int(
(self.position[1] - self.simulation_area.bottom_edge)
// self.simulation_area.cell_height
)
# Ensure the cell indices are within the grid
if cell_i < 0:
cell_i = 0
elif cell_i >= self.simulation_area.cell_grid[0]:
cell_i = self.simulation_area.cell_grid[0] - 1
if cell_j < 0:
cell_j = 0
elif cell_j >= self.simulation_area.cell_grid[1]:
cell_j = self.simulation_area.cell_grid[1] - 1
# Add the particle if it's not already in the
# cell (first iteration)
if not self.cell:
self.cell = (cell_i, cell_j)
self.simulation_area.cells[self.cell].append(self)
# ...
Code Block #30
# particles_simulation.py
# ...
class Particle:
# ...
def assign_to_cell(self):
# Convert the x and y coordinates to cell indices i and j
cell_i = int(
(self.position[0] - self.simulation_area.left_edge)
// self.simulation_area.cell_width
)
cell_j = int(
(self.position[1] - self.simulation_area.bottom_edge)
// self.simulation_area.cell_height
)
# Ensure the cell indices are within the grid
if cell_i < 0:
cell_i = 0
elif cell_i >= self.simulation_area.cell_grid[0]:
cell_i = self.simulation_area.cell_grid[0] - 1
if cell_j < 0:
cell_j = 0
elif cell_j >= self.simulation_area.cell_grid[1]:
cell_j = self.simulation_area.cell_grid[1] - 1
# If the particle is not already in the cell (first iteration), add it
if not self.cell:
self.cell = (cell_i, cell_j)
self.simulation_area.cells[self.cell].append(self)
elif self.cell != (cell_i, cell_j):
# Remove the particle from the current cell
self.simulation_area.cells[self.cell].remove(self)
# Update the cell assignment
self.cell = (cell_i, cell_j)
# Add the particle to the new cell
self.simulation_area.cells[self.cell].append(self)
# ...
Code Block #31
# particles_simulation.py
# ...
class Particle:
MAX_VELOCITY = 2
def __init__(...)
# ...
self.assign_to_cell()
# ...
def move(self):
# ...
self.assign_to_cell()
# ...
Code Block #32
# run_demo.py
from particles_simulation import SimulationArea, Particle
n_particles = 50
simulation_area = SimulationArea(1000, 700, cell_grid=(10, 7))
# ...
# Temporary checks
sample_particle = particles[0]
# Animation loop
while True:
# ...
# Temporary checks
sample_particle.particle.clear()
sample_particle.particle.write(
f"{sample_particle.cell}", font=("Courier", 20, "bold")
)
simulation_area.update()
Code Block #33
# run_demo.py
from particles_simulation import SimulationArea, Particle
n_particles = 200
simulation_area = SimulationArea(1000, 700, cell_grid=(4, 3))
simulation_area.set_title("Diffusion Demo")
# ...
# Animation loop
while True:
for particle in particles:
# ...
# Temporary code: Visualise one of the cells, (1, 1)
if particle in simulation_area.cells[(1, 1)]:
particle.assign_colour("#6D227A")
else:
particle.assign_colour("#AABFBF")
# Temporary checks
sample_particle.particle.clear()
sample_particle.particle.write(
f"{sample_particle.cell}", font=("Courier", 20, "bold")
)
simulation_area.update()
Code Block #34
# particles_simulation.py
import random
import turtle
class SimulationArea:
# ...
def get_sequence_of_cells(self):
return self.cells.values()
# ...
Code Block #35
# run_demo.py
from particles_simulation import SimulationArea, Particle
n_particles = 500
simulation_area = SimulationArea(1000, 700, cell_grid=(10, 7))
simulation_area.set_title("Diffusion Demo")
particles = []
for _ in range(n_particles):
particle = Particle(simulation_area)
if particle.get_position()[0] < 0:
particle.assign_colour("#1A6B72")
else:
particle.assign_colour("#FDB33B")
particles.append(particle)
# Animation loop
while True:
for cell in simulation_area.get_sequence_of_cells():
for particle in cell:
particle.move()
for other_particle in cell:
if particle.check_collision(other_particle):
particle.collide(other_particle)
elif other_particle in particle.in_collision_with:
particle.remove_collision(other_particle)
simulation_area.update()
Code Block #36
# particles_simulation.py
import random
import turtle
random.seed(0)
class SimulationArea:
# ...
Code Block #37
# run_demo.py
import time
from particles_simulation import SimulationArea, Particle
n_particles = 500
grid = (1, 1)
simulation_area = SimulationArea(1000, 700, cell_grid=grid)
# ...
start_time = time.time()
# Animation loop
for _ in range(500):
for cell in simulation_area.get_sequence_of_cells():
for particle in cell:
particle.move()
for other_particle in cell:
if particle.check_collision(other_particle):
particle.collide(other_particle)
elif other_particle in particle.in_collision_with:
particle.remove_collision(other_particle)
# simulation_area.update()
print(f"Time elapsed with grid {grid}: {time.time() - start_time:.2f} s")
Code Block #38
# run_demo.py
# ...
grid = (4, 3)
# ...
Code Block #39
# run_demo.py
# ...
grid = (10, 7)
# ...
Code Block #40
# run_demo.py
from collections import deque
from particles_simulation import SimulationArea, Particle
n_particles = 200
simulation_area = SimulationArea(1000, 700, cell_grid=(10, 7))
simulation_area.set_title("Diffusion Demo")
particles = []
total_green_particles = 0
total_orange_particles = 0
for _ in range(n_particles):
particle = Particle(simulation_area)
if particle.get_position()[0] < 0:
particle.assign_colour("#1A6B72")
particle.colour = "green"
total_green_particles += 1
else:
particle.assign_colour("#FDB33B")
particle.colour = "orange"
total_orange_particles += 1
particles.append(particle)
# Animation loop
green_mean_x = deque(maxlen=500)
orange_mean_x = deque(maxlen=500)
while True:
sum_green_x_coords, sum_orange_x_coords = 0, 0
for cell in simulation_area.get_sequence_of_cells():
for particle in cell:
particle.move()
if particle.colour == "green":
sum_green_x_coords += particle.position[0]
elif particle.colour == "orange":
sum_orange_x_coords += particle.position[0]
for other_particle in cell:
if particle.check_collision(other_particle):
particle.collide(other_particle)
elif other_particle in particle.in_collision_with:
particle.remove_collision(other_particle)
green_mean_x.append(sum_green_x_coords / total_green_particles)
orange_mean_x.append(sum_orange_x_coords / total_orange_particles)
simulation_area.set_title(
f"Green: {sum(green_mean_x)/len(green_mean_x):0.0f}, "
f"Orange: {sum(orange_mean_x)/len(orange_mean_x):0.0f}"
)
simulation_area.update()
All this has absolutely nothing to do with my day to day work, yet I find your post fascinating!
It gives me the urge to try 😁
Thank you for sharing!!