In this chapter we will present examples of so-called open boundary problems. By this we mean that the simulation contains all particles relevant to the problem, and the size of the simulation region is adjusted accordingly at each timestep. This type of problem is by far the easiest one to which one can apply hierarchical data structures. The additional difficulties posed by periodic boundaries will be considered in Chapters 5 and 6.

Gravitational Problems in Astrophysics

Hierarchical tree codes were first developed in the context of astrophysics. This is not surprising because there is a big discrepancy between the number of bodies one would like to study – for example, O (1011) for a galaxy – and the number one can afford to model with a standard N-body code – at present O(105). PIC codes, which employ a grid structure to represent the fields in space, usually cannot handle these problems for two reasons: the complex structure of the investigated object and the large density contrasts such as those found in galaxies. N-body codes are able to avoid the first of these difficulties, because they are gridless and can therefore cope with arbitrarily complicated structures. The second difficulty remains, however, because of the N2 scaling of computation time. This can have two consequences. If the number of simulation particles is too small, the spatial resolution of the simulation might not be good enough to reveal the real dynamic behaviour of the system.

In this chapter we will consider a variety of fields where the tree algorithm in combination with periodic boundaries can be applied and where the speedup in comparison to standard MD and MC codes enables previously inaccessible problems to be investigated. This set of applications is in no way exhaustive, but is intended to indicate the types of problems where the algorithm might best be put to use.

Practically every N-body MD or MC code could incorporate the tree algorithm. However, for systems with short-range potentials, like the Lennard–Jones potentials for the description of solids, this does not bring much of an advantage. Because the potential falls off very rapidly (see Fig. 6.1), a sharp cutoff can be used to exclude interactions of more distant particles whose contribution is negligible. The tree algorithm is mainly suited to systems with long-range forces such as Coulomb, where the summed effect of distant particles is important.

An ideal application is dense, fully-ionized plasmas. Here the particles are so closely packed that an analytical treatment is difficult, and many MD and MC calculations have been carried out to investigate their properties. Limitations in the number of simulation particles make some problems difficult to address due to a combination of poor statistics and small system size. The tree algorithm could be successfully applied here, because the particles interact purely through Coulomb forces.

In previous chapters we occasionally referred to an alternative type of tree code, namely the Fast Multipole Method (FMM). This technique, an elegant refinement of the basic Barnes–Hut algorithm, appears to be best suited to ‘static’ problems, where the particle distribution is more or less uniform. Although it has not been as widely used as the Barnes–Hut (BH) method for dynamic problems – because of either its increased mathematical complexity or the additional computational overhead – it may well become the basis of ‘multimillion’ N-body problems in the near future. We therefore include an introduction to FMM here, based primarily on works by Greengard (1987, 1988, 1990) and Schmidt and Lee (1991). At the same time, we will try to maintain a consistency of notation with the Barnes–Hut algorithm (hereafter referred to as the ‘tree method’ or ‘tree algorithm’), as described in Chapter 2.

Outline of the Fast Multipole Algorithm

The Fast Multipole Method makes use of the fact that a multipole expansion to infinite order contains the total information of a particle distribution. As in the BH algorithm, the interaction between near neighbours is calculated by direct particle–particle force summation, and more distant particles are treated separately. However, the distinction between these two contributions is obtained in a different way. In FMM the distant region is treated as a single ‘far-field’ contribution, which is calculated by a high-order multipole expansion.

The FMM was first formulated by Greengard and Rokhlin (1987).

The difficulty in writing a ‘how-to’ book on numerical methods is to find a form which is accessible to people from various scientific backgrounds. When we started this project, hierarchical N-body techniques were deemed to be ‘too new’ for a book. On the other hand, a few minutes browsing in the References will reveal that the scientific output arising from the original papers of Barnes and Hut (1986) and Greengard and Rohklin (1987) is impressive but largely confined to two or three specialist fields. To us, this suggests that it is about time these techniques became better known in other fields where N-body problems thrive, not least in our own field of computational plasma physics. This book is therefore an attempt to gather everything hierarchical under one roof, and then to indicate how and where tree methods might be used in the reader's own research field. Inevitably, this has resulted in something of a pot-pourri of techniques and applications, but we hope there is enough here to satisfy the beginners and connoisseurs alike.

Tree Construction

We have seen in the preceding chapter that in grid-based codes the particles interact via some averaged density distribution. This enables one to calculate the influence of a number of particles represented by a cell on its neighbouring cells. Problems occur if the density contrast in the simulation becomes very large or the geometry of the problem is very complex.

So why does one bother with a grid at all and not just calculate the interparticle forces? The answer is simply that the computational effort involved quite dramatically limits the number of particles that can be simulated. Particularly with 1/r-type potentials, calculating each particle–particle interaction requires an unnecessary amount of work because the individual contributions of distant particles is small. On the other hand, gridless codes cannot distinguish between near-neighbours and more distant particles; each particle is given the same weighting.

Ideally, the calculation would be performed without a grid in the usual sense, but with some division of the physical space that maintains a relationship between each particle and its neighbours. The force could then be calculated by direct integration while combining increasingly large groups of particles at larger distances. Barnes and Hut (1986) observed that this works in the same way that humans interact with neighbouring individuals, more distant villages, and larger states and countries. A resident of Lower-Wobbleton, Kent, England, is unlikely to undertake a trip to Oberfriedrichsheim, Bavaria, Germany, for a beer and to catch up on the local gossip.

In Chapter 2 we saw the basic workings of the tree algorithm. Now we will discuss some methods that can be used to optimise the performance of this type of code. Although most of these techniques are not specific to tree codes, they are not always straightforward to implement within a hierarchical data structure. It therefore seems worthwhile to reconsider some of the common tricks of the N-body trade, in order to make sure that the tree code is optimised in every sense – not just in its N log N scaling.

There are basically two points of possible improvement:

• Improvement of the accuracy of the particle trajectory calculation by means of higher order integration schemes and individual timesteps. This is especially important for problems involving many close encounters of the particles, that is, ‘collisional’ problems.

• Speedup of the computation time needed to evaluate a single timestep by use of modern software and hardware combinations, such as vectorisation, and special-purpose hierarchical or parallel computer architecture.

Individual Timesteps

For most many-body simulations one would like the total simulated time T = ntΔt (where nt is the number of timesteps) to be as large as possible to approach the hydrodynamic limit. However, the choice of the timestep Δt has to be a compromise between this aim and the fact that as Δt increases, the accuracy of the numerical integration gets rapidly worse.

The hierarchical tree method can not be adapted only for Monte Carlo applications: It can also be modified to perform near-neighbour searches efficiently. This means that the tree algorithm could also have applications for systems with short-range or contact forces. Hernquist and Katz (1989) first showed how the tree structure can be used to find near neighbours through range searching. Following their method, the near-neighbour search is performed the following way.

Consider a system in which only neighbours lying within a distance h will interact with the particle i in question. For the near-neighbour search this sphere is enclosed in a cube whose sides are of length 2h. The tree is built the usual way and the tree search starts at the root. The tree search is performed in a very similar way to the normal force calculation of Section 2.2 by constructing an interaction list. The main difference is that the s/d criterion is substituted by the question: ‘Does the volume of the search cube overlap with the volume of the pseudoparticle presently under consideration?’

If there is no overlap, this branch of the tree contains no near neighbours and is not searched any further. However, if there is an overlap, the cell is subdivided into its daughter cells and the search continues on the next highest level. If the cell is a leaf – meaning there is only one particle in the cell – one has to check whether it actually lies within the radius h of particle i.

Today, the study of The Structure of the Sun is one of the most exciting and rapidly evolving fields in physics. Helioseismology has provided us with a new tool to measure the physical state of the interior of a star, our Sun. This technique is successful to a depth of 0.7 R⊙ (i.e. 0.3 R⊙ from the centre). Deeper than this, observational data has been scarce. However, data are now becoming available from Earth-bound helioseismic networks (GONG, TON, IRIS, BISON,…) and from experiments on board SOHO (GOLF, MDI, VIRGO). These should allow the spectrum of gravity modes for the Sun to be determined, and thus the physical state of the solar core.

This book provides an up-to-date and comprehensive review of our current understanding of the Sun. Each chapter is written by a world expert. They are based on lectures given at the VIth Canary Islands Winter School on Astrophysics. This timely conference brought together leading scientists in the field, postgraduates and recent postdocs students. The aim was to take stock of the new understanding of the Sun and to focus on avenues for fruitful future research. Eight lecturers, around 60 students, and staff from the IAC met in the Hotel Gran Tinerfe in Playa de las Américas (Adeje, Tenerife) from the 5th to the 16th of December, 1994. It was a fortnight of intense and enjoyable scientific work.

INTRODUCTION

This is almost an impossible task, to summarize the subject of Global Changes in the Sun so I must apologize in advance for limiting the scope of these lectures to issues that have been choosen, in part, because of personal interests. I hope that the references provide the reader with footpoints from which to explore a larger set of questions which bear on this subject.

Here we will not discuss the very long, evolutionary, timescales over which the sun changes, nor will we explore the fast changes associated with flares and other transient phenomena. While these discussions depend on some results from MHD models of the solar magnetic cycle, we will not be concerned with the MHD mechanism. These lectures will not address the questions needed to understand local physical models that describe, for example, granulation. On the other hand we will describe some of the physical problems that connect the small-scale behavior of the sun to its global properties. By “global property” I mean an observable that is connected by physically important timescales to the entire sun: limb shape and brightness, largescale magnetic field, oscillation frequencies, solar luminosity, and solar irradiance are all examples of global properties.

Here we are interested in understanding the deviations of the sun from some standard one-dimensional static stellar model. This is a subject we can hardly approach for other stars, and for the sun it is difficult because the physics of magnetic fields and convection are linked over a wide range of spatial and temporal scales.