A glimpse into a world of 2D particles

Modeling Forces

Earlier we briefly hinted at a table to represent pairwise forces. here we make it more explicit by considering it as a matrix called the force matrix. following is an example of this matrix for a system with three particle types

Let's call this matrix :

Each entry represents the force exerted by a particle of color associated with row towards a particle of color associated with column . Diagonal elements (e.g. ) define self-interaction between particles of the same color. For instance in the above example

• 
  means red particles attract each other by a force with a magnitude of
• 
  means red particles are attracted towards green particles by a force with a magnitude of . (a little weaker than the force from before)
• 
  means blue particles repel green ones by a force with a magnitude of

So far for a pair of particles we know how much force is determined by their types, but if this is the only source of force and assuming we have only 2 types that attract each other, after a while all particles collapse into each other. basically we need to model a global repulsive force as well so when particles are too close to one another they repel each other regardless of their type. there are so many possibilities but here we use a very simple one. we model a multiplier which is a function of the distance between particles. let's call it . we use a simple piecewise linear function with a range of where for particles closer than a short distance it gives a negative value regardless of their type (indicating repulsion). In modeling we have 3 key distances

• Min Cutoff distance
  For particles any closer we clamp the repulsion force to -1
• Zero Force distance
  A distance at which particles do not apply any force on one another, but any closer they repel and any farther they attract each other
• Max Cutoff distance
  For particles any farther we assume the force between them is negligible i.e. 0

Here's an illustration of this multiplier function

Using , , and for the three distances introduced above, the multiplier in explicit form is:

With and as the row and column of for the colors of particles and , the force strength from particle towards particle is:

Here is the distance between particles and at time , and the sign of determines attraction () or repulsion ().

We can then define force vector , force from particle towards particle as the following

The net force on a particle is the sum of all such vectors i.e. forces coming from all the other particles towards

To understand the vector math above a little better imagine we have only 4 particles , , , and . we want to compute the overal force on particle at a any given time depending on where all the particles are placed and what kinds of attraction and repulsion forces they apply to each other. let's further assume at a given time particle is somewhere in the middle of the others as illustrated in the figure below

Here

• Arrows with solid lines represent (difference) vectors between each pair of particles that include particle . (these are the only pairs that affect )
• The arrow with dotted line represents the negated vector between particles and , or more precisely . as you can see in the formula to compute this points towards the direction when or when . since we are are talking about forces applied to particle , this represents a repel force from particle
• Dashed arrows represent vectors and (to reduce visual clutter, arrow notation denoting vector quantities has been omitted from the graphics)

Let's further assume green and red both attract orange, but blue repels it (each can have a different constant factor). with these assumptions we can compute , , and and their sum looks like the orange vector below. Plugging in this force into our modeling of motion we can compute in which direction and by how much particle moves and will be positioned at time

Implementation

To better understand how the theoretical model translates into a working simulation, we go one step closer towards the final visualisation by mapping out a flowchart of the simulation logic. this simulation is essentially governed by the modeling we already discussed. so for instance "Compute net force" step is a short way of saying we compute

Thoughts on Performance

So far we saw the modeling and a demo, but didn't really discuss implementation details. the flowchart from earlier does give us an overal picture however. we have a long running loop and we want each iteration not to take any longer than 16 milliseconds or else it means our framerate drops to below 60fps. part of this has to do with the complexity of the problem itself i.e. the number of particles involved. no matter what we do there is always a large enough number of particles our implementation cannot handle efficiently. so let's focus on our 3rd objective (1k particles at 60fps on a modern Mac Mini) and see what we can do to support simulating that many particles.

Why 60 FPS and scale matter?

Our goal isn’t just to visualize particles, it’s to create a smooth 60fps simulation where you can watch complex behaviors emerge in real time. Think of it like a video game. If the framerate stutters, interactivity almost disappears. Every frame requires calculating forces between every pair of particles, which can quickly tank performance. So basically with each iteration of our implementation we are interested in knowing how many particles can we simulate before framerate drops below 60fps.

Let's briefly discuss the performance of a naive implementation of our modeling. this part requires some familiarity with typescript and p5js. (I am using p5js in instance mode so calls to p5 specific API look like p.createVector instead of createVector

In the spirit of OOP let's start with defining a couple of interfaces that we are going to use in our main code

We can then define our update logic like the following

The three methods used within updateParticles are as follows

• computeForceBetweenParticles
  This is where we compute force on particle from particle
• updatePhysics
  This is where we compute new velocity and position of a particle given a net force vector according to our modeling of motion
• handleBoundaryConditions
  This is where we constrain particle motion based on reflecting boundary conditions

Here's how they can be implemented

Then we have

There are several options how to define boundary conditions and each has its own consequences and meanings. In this post for simplicity we only consider reflecting boundary conditions (without any energy loss) which can be implemented simply as below

In part 2 we will also explore implementation details of periodic boundary conditions.

Now that we have familiarised ourselves with the code, first point that comes to mind is that within updateParticles we are looking at every pair twice. once computing force from particle towards and vice versa. these forces are not necessarily equal as you'd expect because the colorForceMatrix is not necessarily a symmetric one. However, one major (heavy) computation here is computing the distance between the two particles and in computing both and this distance is the same. So we better only compute it once per pair. Let's refactor updateParticles based on this observation (notice the highlighted lines)

With this simple change we reduce redundant force calculations by roughly half. This is just the start.

Now let's look at something a little more subtle. within updateParticles we use p.createVector and a couple of related methods such as .sub(), .mag(), and .setMag(). This essentialy incurs some object creation overhead and adds pressure to javascript engine's garbage collection. We could simply get rid of these and use direct numerical operations and primitive types. let's see how that looks

Even though in the first step we directly targetted redundant computations, performance effects from this new update were actually more noticable. This shows how much time savings within nested loops can impact performance.

The next small but useful step (V3) is to hoist simulation constants (cutoff distances, the force matrix) out of the nested loop into local const variables before the loop starts. Accessing a local variable is cheaper than a property lookup inside a hot loop that runs millions of times per frame.

V4 then goes one step further and looks at how we store particle data itself. OOP is nice and all but to speed up further we can look for more memory lookup efficiency. This can be achieved by storing particle data in a flat array instead of an array of objects. Here's how it looks

Other methods of course need to be adapted as we no longer have objects of type Particle.

With these optimisations we can easily render slightly over 2000 particles at 60fps on a Mac Mini M4 Pro which was our initial objective.

Examples

In this section I want to share a couple of examples of how this simulation can be used to create interesting visual effects. in all of these examples the simulation parameters are kept the same except for the color force matrix. Initially particles are distrubuted around a circle and the boundary conditions make it so they are constrained within a white rectangular box.