Numerical Methods

The main techniques I use to study astrophysical fluid dynamics are Smoothed Particle Hydrodynamics (SPH) and the publically-available ZEUS grid-based codes for modelling magneto-hydrodynamics (MHD). SPH is useful for very dynamic problems involving high density contrasts (e.g. gravitational collapse and fragmentation of molecular clouds). I use ZEUS for problems that are closer to equilibrium when I am particularly interested in more subtle effects (e.g. the propagation of waves in accretion discs).

Smoothed Particle Hydrodynamics

SPH is a gridless, Lagrangian particle method typically used to simulate self-gravitating astrophysical fluids. It was invented by L. Lucy in 1977 and extensively developed by R.A. Gingold and J.J. Monaghan from 1977 onwards. Although originally invented for astrophysical processes its use is not limited to astrophysics or to fluids; the technique has been applied to problems with equations of state corresponding liquids and solids and to problems such as the modelling of impacts, volcanic eruptions, and tsunami.

The SPH code I use was originally developed by W. Benz in the 1980s. I have extensively altered the original version and it now includes several important features:

• Individual particle timesteps:
  The timesteps required for individual particles to be evolved correctly are binned into multiples of two and evaluated in these groups. This allows calculations with a large range of dynamical times to be evolved efficiently. This is described in: Bate (1995) and Bate, Bonnell & Price (1995) .

• Sink particles:
  Even evolving particles with individual timesteps the time taken to calculate the evolution of systems with very dense objects (e.g. stars embedded in tenous gas clouds) can be prohibitively long. Therefore, Bate, Bonnell & Price (1995) developed a method by which self-gravitating objects consisting of many SPH gas particles can be replaced (during a calculation) by a single, massive, `sink' particle. This particle interacts with surrounding gas particles via gravity and, if necessary, boundary conditions. Any gas particles that pass within an the `accretion radius' of the sink particle and are bound to the sink particle are accreted. Their mass, linear and angular momemta are added to those of the sink particle. In this way, calculations with dense objects can be evolved over long times without having to model the internal dynamics of the dense objects themselves.

• Parallelisation:
  In 1999, I parallelised my SPH code using OpenMP. This code was used on the UKAFF and GRAND SGI Origin and IBM Power5 parallel supercomputers. The efficiency of the code depends on the specific problem and on the number of particles. However, for simulations without sink particles and more than 1 million particles, it achieves close to linear speed-up for up to 64 processors.
  In 2007, I parallelised the SPH code using MPI. The latest version of the code can run as MPI, OpenMP, or hybrid using MPI between compute nodes (or within nodes) and OpenMP within shared-memory nodes.

• Radiation Hydrodynamics:
  During 2003-2006, my PhD student, Stuart Whitehouse, and I developed a new method for modelling two-temperature (matter and radiation) radiation hydrodynamics in the flux-limited diffusion approximation within the SPH formalism. This is now an option in the SPH code. See our papers: Whitehouse & Bate (2004), Whitehouse, Bate & Monaghan (2005), and Whitehouse & Bate (2006).

• Magnetohydrodynamics (MHD):
  During 2005-2006, Daniel Price added his method for modelling MHD with SPH to the code. Initially, we used an Euler potential method to satisfy div(B)=0 (from 2007-2009), but more recently we are using a divergence cleaning technique (2012+) developed by Tricco & Price (2012). See some of our joint papers: Price & Bate (2007), Price & Bate (2008), Price & Bate (2009), Price, Tricco & Bate (2012), Price, Tricco & Bate (2013).

Along with developing my SPH code itself, I have performed several comparisons of the SPH technique with more traditional grid-based fluid dynamics methods (e.g. Bate & Burkert 1997 ; Burkert, Bate & Bodenheimer 1997 ).

In Bate & Burkert 1997 we discovered that in order to model self-gravitating fluids correctly with SPH, a resolution criterion must be obeyed. We found that the Jeans mass must always be resolved by at least twice the number of particles in a smoothing kernel (typically 50). If this criterion is not obeyed then, depending on the specific implementation of an SPH code, gravitationally-unstable objects may be stablised against collapse or stable objects may collapse unphysically. This Jeans mass criterion is similar to the Jeans length criterion for grid-based codes identified by Truelove et al. (1997) but was discovered independently.

ZEUS: A Finite-difference Magnetohydrodynamics Code

ZEUS is a grid-based finite-difference code for simulating fluids. It was originally developed by J. Stone and M. Norman (1992). The original version is a two-dimensional code that includes ideal magnetohydrodynamics and radiative transfer in the flux-limited diffusion approximation. Several three-dimensional versions have been developed. Currently, I make use of the original 2-D code (which I have parallelised using OpenMP) and a 3-D OpenMP parallel version of ZEUS developed by K. Millar and J. Stone for spherical polar coordinates.