Chapter 3
Particle-Particle method

1. Introduction

The Particle-Particle method is the oldest and most basic particle simulation technique. It is a very straight forward way to describe the behaviour of a system of particles that attract or repulse each other through pairwise interactions.

The pairwise force from particle i on particle j is given by :

Ž
F

_{iŽ j} = -

Ž
F

_{jŽ i} = c·

q₁·q₂

|r_{iŽ j}|³

Ž
r

_{iŽ j}

(1)

In this formula c is a constant that determines the direction and scale of the force : negative values result in attracting forces, positive values in repulsive forces. The factor q is some property of the interacting particle (i.e. mass or charge) and |r_{iŽ j}| the length of the vector connecting both particles. The last term calculates the different components of the force in different directions ( X , Y and possibly Z ).

The total force on one particle in a system of N particles is the vector sum of all separate contributions as given by 1 , resulting in :

Ž
F

_{i_tot} =

N
ĺ
i š j

F_{iŽ j}

(2)

If we know the total force acting on a particle at a certain time, we can use this to calculate the acceleration working on that particle :

Ž
a

_{i_tot} =

Ž
F

_{i_tot}

m_i

(3)

Note that this m is always the mass of the particle, no matter if we calculate gravitational or Coulomb forces ! It may sound trivial, but I saw people using the charge instead : this happens easily if we calculate the acceleration directly, without using the force as an intermediate. This is correct for gravitational forces, but not for Coulomb forces !!

After that we can determine the new speed and position using one of the integration schemes as discussed in . Roughly the complete algorithm looks as in table :

- calculate the total energy and momentum
- set the timecounter to 0
- while ( timecounter < endtime )
  - for all particles do
    - set running total for acceleration a to 0
    - for all particles except the current one do
      - calculate the distances in all directions ( r_x , r_y and r_z )
      - calculate the length of the vector r = Ö{r_x²+r_y²+(r_z²)}
      - calculate the acceleration contribution
      and add to the running total for a
    - endfor
    - calculate new speed and position using a suitable scheme
  - endfor
  - write or plot the new particle positions
  - increase timecounter with the timestep
- endwhile
- check the energy and momentum for conservation

Table 1: Flowchart Particle-Particle code

2. Implementation aspects

During development several minor improvements were added :

The code uses a signal handler to force the calculation of the end momentum and energy and shows this together with the starting energy and momentum and some other information as soon as the simulation is ended prematurely (i.e. by pressing Ctrl-C).
To be able to use arbitrary testfiles, all necessary arrays are being dynamically allocated. This required a few extra lines, but it made the code much more flexible.
Particles with zero mass are skipped in the total force calculation and in the energy and momentum calculation. This is done to provide a mechanism through which we may perform a non-elastic collision, without the need to remove particles physically from our particle list.

2.1 Tunable parameters

There are many variables that can be used to alter the specific behaviour of the simulation. In the datafile the stepsize, the number of iterations and the number of iterations between a screen update can be specified. For the exact meaning of the datafile layout, see the appropriate appendix in .

Commandline options :

00.00.0000

[: -i < nr > ]it is possible to ask for a specific integration method with the -i option. This needs an additional number, that is currently given by the following enumerator list : DUMMY (0), EULER (1), LEAPFROG (2), WOBBLE (3), RUNGA1 (4), RUNGA2 (5), RKWOBBLE (6)
It is by default set to the Leapfrog integration scheme (which is number 2).

[: -a]This causes the graphics display to be cleared before every new update, so the result looks like an animation instead of a cumulative position plot.

[: -c < nr > ]The type of collisionmanagement to be used. Can be either an elastic collision (1), a non-elastic collision (2) or the use of an extra repulsive component (3). Makes it possible to choose another collision management algorithm for particles that come within the COLLISIONDISTANCE. (See section ) The default is currently to perform a complete elastic collision.

[: -x]This switch suppresses all graphics output. It is useful for i.e. batchprocessing, if only the total running time and total error in energy and momentum are of interest. Without this switch the program can only run if an X-windows capable terminal is used and the DISPLAY variable is set.

[: -f]Skip the energy and momentum computation at the start and the end. Omitting the energy computation isn't of much use with Particle-Particle, but it can save an enormous amount of time in case of the Particle-Mesh and Tree-Code programs for the larger data sets, as the computation of the energy is of order O(N²) !

[: -o < file > ]Write all particle coordinates to the file specified at every plotstep. Note that this not necessarily has to be at every iteration !

At compile time there are several other parameters that can be set :

MAGX and MAGY: These determine the scaling factor that is being used to plot the particle orbits on the screen. They are for historical reasons both set to 150. It is possible to change the aspect ratio by these changing parameters if necessary. All testfiles are chosen in such a way to fit exactly on the plot window if the default scaling factor is used. (See also the appendix about scaling of parameterfiles)
USE3DFORCES: With normal gravitational or electrical problems we use a 3D force, dropping with the square of the distance ( F ~ [1/(r²)] ). But in a 2D simulation using the Particle-Mesh method, the force only drops with the absolute value of the distance ( F ~ [1/| r| ] ). To be able to compare the results of PP and PM, this macro should be undefined. In that case also PP will use forces proportional to [1/| r| ] .
COLLISIONDISTANCE: This is the minimum distance allowed between two particles, to avoid the acceleration gets to high. One of the collision mechanism routines as described in section is called as soon as this distance is reached.
LAMBDA: This is the constant c as used in the interparticle force calculation (see equation 1). It is currently set to the Cavendish constant, that determines the size of a gravitational force. Note that the code assumes attractive forces by default, so the minus sign is embedded in the code itself.
DEBUGLEVEL: A value of 0 is the default: it shows no debugging information. A value of 1 shows some more information about the simulation, but still shows the graphical output as well. A value of 2 shows even more debug information and also disables the graphical output.

Although it is possible to use a variable timestep in the Particle-Particle simulations (small timesteps for the fast particles and large ones for the slower ones), this is not (currently) implemented as it would be useless for the other two simulation methods. (See discussion of variable timesteps in section ).

3. Visualisation techniques

Just a few remarks about the visualisation:

To show the particle orbits or the new particle positions a simple graphic interface - written by myself using X-windows in combination with Athena widgets - is used. This graphic visualisation process forks from the simulation process and runs as a separate thread that communicates with the simulation program through a Unix pipe.
This construction ensures the simulation program has little overhead for graphics: it only should write the new coordinates to the pipe and then the bulk of the work is done in the graphic process. This is a pleasant feature for timing purposes as the processtime spent in the graphic process doesn't influence the processtime in the simulation process.
Animations are achieved by clearing the whole display before every new screen update. Cumulative plots just leave the old positions on the screen as well as the new ones (which is the default).
Plotting just uses a dot or a small circle to visualise the particle positions: it doesn't connect the previous position with the new one as would be done if we'd draw a straight line between them instead. This makes it possible to get some idea of the speed as well: the larger the distance between 2 positions, the higher the speed of that particle. Drawing a straight line would hide this extra information.
On the parallel (Parsytec) computer only the simulation process is run. The graphic output still uses the same X-windows interface as described earlier, but this time it runs separately on the (Unix) host computer. The communication between them is established by connecting the standard output of the simulation process (through a server process on the host computer called the iserver) to the standard input of the visualisation process.
This is done because the only connection between the Parsytec and the host computer is through the stdin, stdout and stderr channels from the server process, so we must use either stdout or stderr and we cannot open a separate pipe.
Scaling and shifting the picture can be done through some extra buttons in the visualisation window, so it is not strictly necessary to scale the problem until it fits the graphic window or to adapt the scaling parameters in the source code.

4. Restrictions

From all simulation methods discussed in this article, Particle-Particle has the least accuracy restrictions of all. However, things may easily go wrong if particles come to close to each other, so we should be certain that our timestep is small enough to have a reliable simulation also for the shortest possible (or allowable !) distances.

But another much more important restriction rises with Particle-Particle simulation: the acceleration calculation for a system of n particles requires roughly n·(n-1) operations (this could be halved on a sequential system if we assume the force from particle i on j would be minus the force from particle j on i ).

This implies for a system of 10 particles about 50 operations per timestep, but for a system of 100000 particles about 5000000000 operations per timestep ! Of course this makes Particle-Particle simulation completely useless for the simulation of large systems.

The two other methods described in this report both scale with n·log₂n , which is for a system of 10 particles still about 30 operations per timestep, but for a system of 100000 particles now only 1600000 operations per timestep (which is more than 3000 times less as PP) ! Note however that one operation in PM or TC may take up to 10-50 times as much time as one operation under PP.