Methods for computing giant planet formation and evolution

We present a numerical code for computing all stages of the formation and evolution of giant planets in the framework of the core instability mechanism. This code is a non-trivial adaption of the stellar binary evolution code and is based on a standard Henyey technique. To investigate the performance of this code we applied it to the computation of the formation and evolution of a Jupiter mass object from a half Earth core mass to ages in excess of the age of the Universe.
We also present a new smoothed linear interpolation algorithm devised especially for the purpose of circumventing some problems found when some physical data (e.g. opacities, equation of state, etc.) are introduced into an implicit algorithm like the one employed in this work.     

A New Code for Nonradial Stellar Pulsations and its Application to Low-Mass, Helium White Dwarfs

We present a finite difference code intended for computing linear, adiabatic, non radial pulsations of spherical stars. This code is based on a slight modification of the general Newton-Raphson technique in order to handle the relaxation of the eigenvalue(square of the eigenfrequency) of the modes and their corresponding eigenfunctions. This code has been tested computing the pulsation spectra of polytropic spheres finding a good agreement with previous work. Then, we have coupled this code to our evolutionary code and applied it to the computation of the pulsation spectrum of a low mass, pure-helium white dwarf of 0.3 M _⊙ for a wide range of effective temperatures. In making this calculation we have taken an evolutionary time step short enough such that eigenmodes corresponding to a given model are used as initial approximation to those of the next one. Specifically, we have computed periods, period spacing, eigenfunctions, weight functions, kinetic energies and variational periods for a wide range of modes. To our notice this is the first effort in studying the pulsation properties of helium white dwarfs. The solution we have found working with these realistic white dwarf models are in good accord with the predictions of the asymptotic theory of Tassoul (1980) for high order modes. This indicates that the code presented here is able to work adequately also with realistic stellar models. This revised version was published online in July 2006 with corrections to the Cover Date.     

A numerical scheme and benchmark tests for non-isothermal two-fluid ambipolar diffusion

Stars form in magnetized molecular clouds composed primarily of neutral gas with a trace amount of ions. We present a semi-implicit strategy for incorporating the equations that describe the coupled ion and neutral two-fluid equations, with a full energy equation, into the RIEMANN code that uses a TR-BDF2 algorithm to stably handle the stiffness of the source terms. We demonstrate that the numerical implementation works through the use of a suite of test problems that we catalog here. We show that reproducing the analytic dispersion analysis for the propagation of waves in a two-fluid plasma is an especially strong code test. We also present a two-fluid analogue of the Noh wall-shock problem and demonstrate the performance of the code on the Wardle instability. We also present a novel blast wave test, showing that the results reduce to the single fluid results under strong coupling, yet differing considerably when the coupling is weak. These test problems demonstrate that the numerical implementation can accurately capture the dissipation rate of waves and reproduce the structure of a C-shock.     

N-body integrators with individual time steps from Hierarchical splitting

We review the implementation of individual particle time-stepping for N-body dynamics. We present a class of integrators derived from second order Hamiltonian splitting. In contrast to the usual implementation of individual time-stepping, these integrators are momentum conserving and show excellent energy conservation in conjunction with a symmetrized time step criterion. We use an explicit but approximate formula for the time symmetrization that is compatible with the use of individual time steps. No iterative scheme is necessary. We implement these ideas in the HUAYNO¹ code and present tests of the integrators and show that the presented integration schemes shows good energy conservation, with little or no systematic drift, while conserving momentum and angular momentum to machine precision for long term integrations.     

The Adaptive TreePM: an adaptive resolution code for cosmological N-body simulations

Cosmological N -body simulations are used for a variety of applications. Indeed progress in the study of large-scale structures and galaxy formation would have been very limited without this tool. For nearly 20 yr the limitations imposed by computing power forced simulators to ignore some of the basic requirements for modelling gravitational instability. One of the limitations of most cosmological codes has been the use of a force softening length that is much smaller than the typical interparticle separation. This leads to departures from collisionless evolution that is desired in these simulations. We propose a particle-based method with an adaptive resolution where the force softening length is reduced in high-density regions while ensuring that it remains well above the local interparticle separation. The method, called the Adaptive TreePM (ATreePM), is based on the TreePM code. We present the mathematical model and an implementation of this code, and demonstrate that the results converge over a range of options for parameters introduced in generalizing the code from the TreePM code. We explicitly demonstrate collisionless evolution in collapse of an oblique plane wave. We compare the code with the fixed resolution TreePM code and also an implementation that mimics adaptive mesh refinement methods and comment on the agreement and disagreements in the results. We find that in most respects the ATreePM code performs at least as well as the fixed resolution TreePM in highly overdense regions, from clustering and number density of haloes to internal dynamics of haloes. We also show that the adaptive code is faster than the corresponding high-resolution TreePM code.     

An implementation of radiative transfer in the cosmological simulation code gadget

We present a novel numerical implementation of radiative transfer in the cosmological smoothed particle hydrodynamics (SPH) simulation code gadget . It is based on a fast, robust and photon-conserving integration scheme where the radiation transport problem is approximated in terms of moments of the transfer equation and by using a variable Eddington tensor as a closure relation, following the Optically Thin Variable Eddington Tensor suggestion of Gnedin & Abel. We derive a suitable anisotropic diffusion operator for use in the SPH discretization of the local photon transport, and we combine this with an implicit solver that guarantees robustness and photon conservation. This entails a matrix inversion problem of a huge, sparsely populated matrix that is distributed in memory in our parallel code. We solve this task iteratively with a conjugate gradient scheme. Finally, to model photon sink processes we consider ionization and recombination processes of hydrogen, which is represented with a chemical network that is evolved with an implicit time integration scheme. We present several tests of our implementation, including single and multiple sources in static uniform density fields with and without temperature evolution, shadowing by a dense clump and multiple sources in a static cosmological density field. All tests agree quite well with analytical computations or with predictions from other radiative transfer codes, except for shadowing. However, unlike most other radiative transfer codes presently in use for studying re-ionization, our new method can be used on-the-fly during dynamical cosmological simulation, allowing simultaneous treatments of galaxy formation and the re-ionization process of the Universe.     

On the dynamical evolution of the brown dwarf populationin open clusters

It is by now well established that open clusters contain a considerable fraction of brown dwarfs (BDs). This paper investigates the dynamical evolution of this substellar population by using simulations with Aarseth's (1994) NBODY5 code. A noticeable preferential escape of BDs is found, which may influence the determination of the IMF of substellar objects in dynamically evolved open clusters. This small dynamical-in-origin depletion may not explain, however, the scarcity of BDs observed in some evolved clusters, as the Hyades. On the other hand, BD cooling processes are able to reduce our ability to detect BDs in old clusters in a very significant way. Our results confirm that the probability of observing BDs in open clusters is almost the same over the whole cluster area because they are distributed quite uniformly even at late stages of the evolution of the cluster. This is expected to be a general feature as observed for low-mass stars in well studied open clusters (Pleiades, Praesepe). Our present calculations show that clusters as old as the Pleiades may have lost about 10% of their initial BD population but the number ratio of BDs to normal (not substellar) stars must remain almost unchanged. However, the long-term behavior of the relative percentage of BDs depends strongly on the initial mass function (IMF) assumed in the calculations. Clusters with a Salpeterian IMF evolve to reach relative percentages of BDs as low as 40% for a starting value around 70%. Our results suggest that BDs in clusters escape preferentially by evaporation. This revised version was published online in July 2006 with corrections to the Cover Date.     

crash2: coloured packets and other updates

In this paper, we report on the improvements implemented in the cosmological radiative transfer code crash . In particular, we present a new multi-frequency algorithm for spectra sampling which makes use of coloured photon packets: we discuss the need for the multi-frequency approach, describe its implementation and present the improved crash performance in reproducing the effects of ionizing radiation with an arbitrary spectrum. We further discuss minor changes in the code implementation which allow for more efficient performance and an increased precision.     

A GPU accelerated Barnes–Hut tree code for FLASH4

We present a GPU accelerated CUDA-C implementation of the Barnes Hut (BH) tree code for calculating the gravitational potential on octree adaptive meshes. The tree code algorithm is implemented within the FLASH4 adaptive mesh refinement (AMR) code framework and therefore fully MPI parallel. We describe the algorithm and present test results that demonstrate its accuracy and performance in comparison to the algorithms available in the current FLASH4 version. We use a MacLaurin spheroid to test the accuracy of our new implementation and use spherical, collapsing cloud cores with effective AMR to carry out performance tests also in comparison with previous gravity solvers. Depending on the setup and the GPU/CPU ratio, we find a speedup for the gravity unit of at least a factor of 3 and up to 60 in comparison to the gravity solvers implemented in the FLASH4 code. We find an overall speedup factor for full simulations of at least factor 1.6 up to a factor of 10.     

Computational advances in gravitational microlensing: A comparison of CPU,GPU, and parallel,large data codes

To assess how future progress in gravitational microlensing computation at high optical depth will rely on both hardware and software solutions, we compare a direct inverse ray-shooting code implemented on a graphics processing unit (GPU) with both a widely-used hierarchical tree code on a single-core CPU, and a recent implementation of a parallel tree code suitable for a CPU-based cluster supercomputer. We examine the accuracy of the tree codes through comparison with a direct code over a much wider range of parameter space than has been feasible before. We demonstrate that all three codes present comparable accuracy, and choice of approach depends on considerations relating to the scale and nature of the microlensing problem under investigation. On current hardware, there is little difference in the processing speed of the single-core CPU tree code and the GPU direct code, however the recent plateau in single-core CPU speeds means the existing tree code is no longer able to take advantage of Moore’s law-like increases in processing speed. Instead, we anticipate a rapid increase in GPU capabilities in the next few years, which is advantageous to the direct code. We suggest that progress in other areas of astrophysical computation may benefit from a transition to GPUs through the use of “brute force” algorithms, rather than attempting to port the current best solution directly to a GPU language – for certain classes of problems, the simple implementation on GPUs may already be no worse than an optimised single-core CPU version.     

Adaptive mesh fluid simulations on GPU

We describe an implementation of compressible inviscid fluid solvers with block-structured adaptive mesh refinement on Graphics Processing Units using NVIDIA’s CUDA. We show that a class of high resolution shock capturing schemes can be mapped naturally on this architecture. Using the method of lines approach with the second order total variation diminishing Runge–Kutta time integration scheme, piecewise linear reconstruction, and a Harten–Lax–van Leer Riemann solver, we achieve an overall speedup of approximately 10 times faster execution on one graphics card as compared to a single core on the host computer. We attain this speedup in uniform grid runs as well as in problems with deep AMR hierarchies. Our framework can readily be applied to more general systems of conservation laws and extended to higher order shock capturing schemes. This is shown directly by an implementation of a magneto-hydrodynamic solver and comparing its performance to the pure hydrodynamic case. Finally, we also combined our CUDA parallel scheme with MPI to make the code run on GPU clusters. Close to ideal speedup is observed on up to four GPUs.     

The determination of the gravitational potential of axisymmetric mass distribution in spherical coordinates. II

The numerical computation of gravitational potential of rotational-symmetric mass distribution in spherical coordinates requires a lot of computing time. We discuss direct methods and present a reduction method which is superior to methods used so far. Results for accuracy which is similar to that of other methods and computing times are given in Fig. 2 and 3.     

A new parallel PM code for very large-scale cosmological simulations

We have developed a parallel Particle–Particle, Particle–Mesh (P³M) simulation code for the Cray T3E parallel supercomputer that is well suited to studying the time evolution of systems of particles interacting via gravity and gas forces in cosmological contexts. The parallel code is based upon the public-domain serial Adaptive P³M-SPH (http://coho.astro.uwo.ca/pub/hydra/hydra.html) code of Couchman et al. (1995)[ApJ, 452, 797]. The algorithm resolves gravitational forces into a long-range component computed by discretizing the mass distribution and solving Poisson's equation on a grid using an FFT convolution method, and a short-range component computed by direct force summation for sufficiently close particle pairs. The code consists primarily of a particle–particle computation parallelized by domain decomposition over blocks of neighbour-cells, a more regular mesh calculation distributed in planes along one dimension, and several transformations between the two distributions. The load balancing of the P³M code is static, since this greatly aids the ongoing implementation of parallel adaptive refinements of the particle and mesh systems. Great care was taken throughout to make optimal use of the available memory, so that a version of the current implementation has been used to simulate systems of up to 10⁹ particles with a 1024³ mesh for the long-range force computation. These are the largest Cosmological N-body simulations of which we are aware. We discuss these memory optimizations as well as those motivated by computational performance. Performance results are very encouraging, and, even without refinements, the code has been used effectively for simulations in which the particle distribution becomes highly clustered as well as for other non-uniform systems of astrophysical interest.     

Numerical simulations of the Hyades dynamics and the nature of the moving Hyades cluster

We present our numerical simulations of the dynamical evolution of the Hyades open cluster. The simulations were performed usinga modified NBODY6 algorithm that included tidal forces and a realistic orbit of the cluster in a gravitational field described by the Miyamoto-Nagai potential. Our goal was to study the nature of movingclu sters. We show that the stars that were earlier cluster members could be later identified within a sphere of 50 pc in diameter around the Sun. The number of such stars for the chosen initial mass and virial radius of the cluster does not exceed ten. The maximum space velocity of these stars relative to the core of the current cluster does not exceed 3 km s^?1. Our numerical simulations confirm the assumption that some of the moving clusters near the Sun could consist of stars that have escaped from open clusters in the course of their dynamical evolution.     

A tree code for planetesimal dynamics: comparison with a hybrid direct code

We present a tree code for simulations of collisional systems dominated by a central mass. We describe the implementation of the code and the results of some test runs with which the performance of the code was tested. A comparison between the behaviour of the tree code and a direct hybrid integrator is also presented. The main result is that tree codes can be useful in numerical simulations of planetary accretion, especially during intermediate stages, where possible runaway accretion and dynamical friction lead to a population with a few large bodies in low-eccentricity and low-inclination orbits embedded in a large swarm of small planetesimals in rather excited orbits. Some strategies to improve the performance of the code are also discussed.     

Numerical MHD codes for modeling astrophysical flows

We describe a Godunov-type magnetohydrodynamic (MHD) code based on the Miyoshi and Kusano (2005) solver which can be used to solve various astrophysical hydrodynamic and MHD problems. The energy equation is in the form of entropy conservation. The code has been implemented on several different coordinate systems: 2.5D axisymmetric cylindrical coordinates, 2D Cartesian coordinates, 2D plane polar coordinates, and fully 3D cylindrical coordinates. Viscosity and diffusivity are implemented in the code to control the accretion rate in the disk and the rate of penetration of the disk matter through the magnetic field lines. The code has been utilized for the numerical investigations of a number of different astrophysical problems, several examples of which are shown.     

A multidimensional magnetohydrodynamics code based on semidiscrete central WENO and PPM schemes

We present a multidimensional magnetohydrodynamics code based on semidiscrete central Godunov-type schemes and high order weighted essentially non-oscillatory data reconstruction and incorporate elements of the piecewise parabolic method into the reconstruction schemes. This algorithm avoids spectral decomposition into characteristic waves and is easily implemented. We have implemented this algorithm and carried out one and two dimensional tests, which verifies that this code can achieve the accuracy and efficiency to solve the astrophysical problems.     

Using Java for distributed computing in the Gaia satellite data processing

In recent years Java has matured to a stable easy-to-use language with the flexibility of an interpreter (for reflection etc.) but the performance and type checking of a compiled language. When we started using Java for astronomical applications around 1999 they were the first of their kind in astronomy. Now a great deal of astronomy software is written in Java as are many business applications. We discuss the current environment and trends concerning the language and present an actual example of scientific use of Java for high-performance distributed computing: ESA’s mission Gaia. The Gaia scanning satellite will perform a galactic census of about 1,000 million objects in our galaxy. The Gaia community has chosen to write its processing software in Java. We explore the manifold reasons for choosing Java for this large science collaboration. Gaia processing is numerically complex but highly distributable, some parts being embarrassingly parallel. We describe the Gaia processing architecture and its realisation in Java. We delve into the astrometric solution which is the most advanced and most complex part of the processing. The Gaia simulator is also written in Java and is the most mature code in the system. This has been successfully running since about 2005 on the supercomputer “Marenostrum” in Barcelona. We relate experiences of using Java on a large shared machine. Finally we discuss Java, including some of its problems, for scientific computing.