An investigation, implementation and comparison of 3 important particle simulation techniques:

PP : Particle-Particle

PM : Particle-Mesh and

TC : Tree-code

Kees Lemmens
Faculty of Applied Physics,
University of Technology,
Delft.

Feb 26, 1998

Version 1.1

Members of committee :

Prof.dr.ir. L. Dekker.
Prof.dr.ir. S.W. de Leeuw.
Dr.ir. J.J. Derksen.
Dr.ir. E.E.E. Frietman.
Dr.ir. H.X. Lin.

This report was completely realised without using any commercial software: all programs were written using the Public Domain JED editor and compiled with one of the best C-compilers available: the GNU-C Compiler.
All drawings were made using the Xfig vector drawing package, the 2D and 3D graphs were generated through GNUplot and the typesetting was done using LATEX 2E with the LYX ``wysiwyg'' frontend. As operating system the excellent Public-Domain Unix system called Linux was used, installed on a Pentium 75 with 32 MB internal memory.

Contents

1  Introduction
    1.1  Overview of discussed multibody simulation techniques
        1.1.1  Particle-Particle method
        1.1.2  Particle-Mesh method
        1.1.3  Tree-Code method
2  General aspects
    2.1  Integration methods
        2.1.1  Brief discussion of several methods
        2.1.2  Using only one equation of motion
        2.1.3  Discretisation errors
        2.1.4  Round-off errors
            2.1.4.1  Euler error-propagation
            2.1.4.2  Leapfrog error-propagation
    2.2  How accurate is our simulation ?
        2.2.1  Momentum conservation
        2.2.2  Energy conservation
    2.3  Long range versus short range forces
        2.3.1  2D and 3D aspects of uniform density
    2.4  How to choose the timestep
        2.4.1  Particle-Mesh considerations
        2.4.2  Variable timesteps
    2.5  Collision management
        2.5.1  Collisionless systems: Particle-Mesh
        2.5.2  Elastic collisions
        2.5.3  Non-Elastic collisions
        2.5.4  Extra repulsive components
3  Particle-Particle method
    3.1  Introduction
    3.2  Implementation aspects
        3.2.1  Tunable parameters
    3.3  Visualisation techniques
    3.4  Restrictions
4  Particle-Mesh method
    4.1  Introduction
    4.2  Mass (Charge) to Mesh assignment
        4.2.1  Nearest-Grid-Point
        4.2.2  Cloud-In-Cell
    4.3  Calculate Potentials on Mesh
    4.4  Fourier Transforms
        4.4.1  Discrete Fourier Transform (DFT)
        4.4.2  Fast Fourier Transform (FFT)
        4.4.3  2D Fourier Transforms
    4.5  Calculate the gravity (or electric) field
        4.5.1  Field interpolation
    4.6  Implementation and practical aspects
        4.6.1  Tunable parameters
    4.7  Visualisation techniques and problems
    4.8  Restrictions
5  Tree-Code method
    5.1  Introduction
    5.2  Algorithm based on Barnes-Hut method
    5.3  Choice of the Q parameter
    5.4  Implementation aspects
        5.4.1  Sorting techniques
        5.4.2  Memory considerations
        5.4.3  Tunable parameters
    5.5  Collision management
    5.6  Visualisation techniques
    5.7  Restrictions
6  Analysis: Which method for what problem ?
    6.1  Introduction
    6.2  Parameters and properties
    6.3  Autocorrelation
        6.3.1  Discretisation
    6.4  Uniform and non-uniform problems
    6.5  Implementation aspects
    6.6  Tunable parameters
    6.7  Restrictions
7  Results and conclusions
    7.1  Results
    7.2  Conclusions

Appendices

A  Parameter files
    A.1  Parameterfile layout
    A.2  Scaling parameter files
        A.2.1  Ekin and Epot both invariant
        A.2.2  Ekin and Epot both scaled with the same factor
B  Some examples of analyser output
    B.1  Example of a 1024 particle mesh problem
    B.2  Example of a 300 particle non-uniform problem
C  Source code
    C.1  Particle-Particle
    C.2  Particle-Mesh
    C.3  Tree-Code
    C.4  Analyser
     Bibliography

Abstract

In this report 3 basic particle simulation techniques are investigated, implemented and compared with each other: the Particle-Particle method, the Particle-Mesh method and a variant of the Barnes-Hut ([4]) Tree-Code algorithm. All of them were implemented in such a way that it was possible to feed the same datafile with information about particle positions, velocities and masses as well as the number of iterations and the size of the timestep to all 3 different programs. This gives the possibility to compare the 3 methods directly with respect to i.e. speed, speedup, accuracy and/or memory requirements.
Many problems that were encountered during the implementation of the 3 programs and of course possible solutions are discussed as well. One of the intentions of this report is to be able to make someone who is not familiar with particle simulations techniques known with the basics in a relatively short time, without going too deep into the hairy and confusing details of the individual methods.

The Particle-Particle and Particle-Mesh techniques are basically the same as described in i.e. the book ``Computer Simulations using particles'' by Hockney and Eastwood ([9]) except that in this book they are almost exclusively discussed for 1D, while I used 2D instead. The step from 2D to the 3D world - which is the most interesting from a practical point of view - is relatively easy, while 1D is too far of from reality to my opinion.
The Tree-Code algorithm discussed in this report is a variant on the Barnes-Hut algorithm, with the most important differences the building of the tree (where I used a method I decided to call ``sort-and-count'') and the definition of the Q criterion (where I use the diagonal length instead of the length of one rib).

Finally a heuristic is presented, that is supposed to be able to analyse a problem within a reasonable time and then give an indication about which of the 3 discussed methods would be best suited for that problem.
To give this indication the heuristic mainly uses the density auto-correlation and the ratio between the average particle displacement during the simulation (AD) and the average minimum interparticle distance (AMI).

List of Figures

    1.1  Example of Particle-Particle simulation.
    1.2  Example of a mesh for the problem from the Particle-Particle section.
    1.3  Quadtree for the problem from the Particle-Particle section.
    2.1  Comparison between a stable method with a large roundoff error and unstable method with a small roundoff error. The tradeoff is at about 24 iterations in this example.
    2.2  Influence of particles at certain distance r® (r+dr) on the total force.
    2.3  Elastic collisions between attracting particles are often almost central.
    4.1  Here the 2D CIC mass assignment is shown in a graphic way.
    4.2  DC part, right- and leftband in the asymmetrical DFT ( n=0 ... N-1 )
    5.1  Tree example for the 4 body problem from chapter 1.3.
    5.2  Discrimination between faraway and nearby cells using Q = 1 (note that r2 » r1)
    5.3  Fast search of the starting point in the sortlist.
    5.4  Cell array to cell tree mapping.
    6.1  Flowchart for analyser program.
    6.2  A small gridpattern and it's autocorrelation function (256 particles)
    6.3  A larger gridpattern and it's autocorrelation function (4096 particles)
    6.4  A small random generated pattern and it's autocorrelation function (256 particles)
    6.5  A larger random generated pattern and it's autocorrelation function (4096 particles)
    6.6  500 particle random set with a few extra heavy masses causing non-uniformity.
    6.7  300 particle set built by clustering 3 smaller random sets of 100 particles each.
    6.8  Autocorrelation for the 4 particle problem from figure 1.1.
    7.1  Output for the 4 body problem from chapter 1 if run with PP (left) and with PM (right). In both cases a 2D force is used to make them comparable.
    7.2  Comparison of the PP, TC and PM runtimes for the problems from tables 7.1 and 7.2.

List of Tables

    3.1  Flowchart Particle-Particle code
    4.1  The relation between the dimension and the potential, field and force for an electrical problem.
    4.2  Flowchart Particle-Mesh code
    4.3  How to FFT a 2-dimensional mesh using only one 1-dimensional FFT routine
    5.1  Flowchart Tree-Code for a 2D system
    7.1  Comparison between the 3 methods for several non-uniform sets.
    7.2  Comparison between the 3 methods for several smaller uniform sets (meshes).
    7.3  Results for generated sets and Tree-Code.
    7.4  Results for generated sets and Particle-Mesh.
    7.5  Total number of computations for both O(N2) and O(N logN) for the generated problems.

Preface

This report describes the results obtained during a research period from Sept 1996 - August 1997.

The research period was part of the Master Degree program from the Faculty of Applied Physics of the Technical University of Delft, to meet the requirements for the final examination.

I would like to thank the following persons for their support during this research period: