Computational Astrophysics and Cosmology最新文献

We present a new version of conservative ADER-WENO finite volume schemes, in which both the high order spatial reconstruction as well as the time evolution of the reconstruction polynomials in the local space-time predictor stage are performed in primitive variables, rather than in conserved ones. To obtain a conservative method, the underlying finite volume scheme is still written in terms of the cell averages of the conserved quantities. Therefore, our new approach performs the spatial WENO reconstruction twice: the first WENO reconstruction is carried out on the known cell averages of the conservative variables. The WENO polynomials are then used at the cell centers to compute point values of the conserved variables, which are subsequently converted into point values of the primitive variables. This is the only place where the conversion from conservative to primitive variables is needed in the new scheme. Then, a second WENO reconstruction is performed on the point values of the primitive variables to obtain piecewise high order reconstruction polynomials of the primitive variables. The reconstruction polynomials are subsequently evolved in time with a novel space-time finite element predictor that is directly applied to the governing PDE written in primitive form. The resulting space-time polynomials of the primitive variables can then be directly used as input for the numerical fluxes at the cell boundaries in the underlying conservative finite volume scheme. Hence, the number of necessary conversions from the conserved to the primitive variables is reduced to just one single conversion at each cell center. We have verified the validity of the new approach over a wide range of hyperbolic systems, including the classical Euler equations of gas dynamics, the special relativistic hydrodynamics (RHD) and ideal magnetohydrodynamics (RMHD) equations, as well as the Baer-Nunziato model for compressible two-phase flows. In all cases we have noticed that the new ADER schemes provide less oscillatory solutions when compared to ADER finite volume schemes based on the reconstruction in conserved variables, especially for the RMHD and the Baer-Nunziato equations. For the RHD and RMHD equations, the overall accuracy is improved and the CPU time is reduced by about 25?%. Because of its increased accuracy and due to the reduced computational cost, we recommend to use this version of ADER as the standard one in the relativistic framework. At the end of the paper, the new approach has also been extended to ADER-DG schemes on space-time adaptive grids (AMR).

In this paper we describe a simple numerical approach which allows to study the structure of steady-state axisymmetric relativistic jets using one-dimensional time-dependent simulations. It is based on the fact that for narrow jets with (v_{z}approx c) the steady-state equations of relativistic magnetohydrodynamics can be accurately approximated by the one-dimensional time-dependent equations after the substitution (z=ct). Since only the time-dependent codes are now publicly available this is a valuable and efficient alternative to the development of a high-specialised code for the time-independent equations. The approach is also much cheaper and more robust compared to the relaxation method. We tested this technique against numerical and analytical solutions found in literature as well as solutions we obtained using the relaxation method and found it sufficiently accurate. In the process, we discovered the reason for the failure of the self-similar analytical model of the jet reconfinement in relatively flat atmospheres and elucidated the nature of radial oscillations of steady-state jets.

Astrophysical direct N-body methods have been one of the first production algorithms to be implemented using NVIDIA’s CUDA architecture. Now, almost seven years later, the GPU is the most used accelerator device in astronomy for simulating stellar systems. In this paper we present the implementation of the Sapporo2 N-body library, which allows researchers to use the GPU for N-body simulations with little to no effort. The first version, released five years ago, is actively used, but lacks advanced features and versatility in numerical precision and support for higher order integrators. In this updated version we have rebuilt the code from scratch and added support for OpenCL, multi-precision and higher order integrators. We show how to tune these codes for different GPU architectures and present how to continue utilizing the GPU optimal even when only a small number of particles ((N < 100)) is integrated. This careful tuning allows Sapporo2 to be faster than Sapporo1 even with the added options and double precision data loads. The code runs on a range of NVIDIA and AMD GPUs in single and double precision accuracy. With the addition of OpenCL support the library is also able to run on CPUs and other accelerators that support OpenCL.

In the implicit large eddy simulation (ILES) paradigm, the dissipative nature of high-resolution shock-capturing schemes is exploited to provide an implicit model of turbulence. The ILES approach has been applied to different contexts, with varying degrees of success. It is the de-facto standard in many astrophysical simulations and in particular in studies of core-collapse supernovae (CCSN). Recent 3D simulations suggest that turbulence might play a crucial role in core-collapse supernova explosions, however the fidelity with which turbulence is simulated in these studies is unclear. Especially considering that the accuracy of ILES for the regime of interest in CCSN, weakly compressible and strongly anisotropic, has not been systematically assessed before. Anisotropy, in particular, could impact the dissipative properties of the flow and enhance the turbulent pressure in the radial direction, favouring the explosion. In this paper we assess the accuracy of ILES using numerical methods most commonly employed in computational astrophysics by means of a number of local simulations of driven, weakly compressible, anisotropic turbulence. Our simulations employ several different methods and span a wide range of resolutions. We report a detailed analysis of the way in which the turbulent cascade is influenced by the numerics. Our results suggest that anisotropy and compressibility in CCSN turbulence have little effect on the turbulent kinetic energy spectrum and a Kolmogorov (k^{-5/3}) scaling is obtained in the inertial range. We find that, on the one hand, the kinetic energy dissipation rate at large scales is correctly captured even at low resolutions, suggesting that very high “effective Reynolds number” can be achieved at the largest scales of the simulation. On the other hand, the dynamics at intermediate scales appears to be completely dominated by the so-called bottleneck effect, i.e., the pile up of kinetic energy close to the dissipation range due to the partial suppression of the energy cascade by numerical viscosity. An inertial range is not recovered until the point where high resolution ～512³, which would be difficult to realize in global simulations, is reached. We discuss the consequences for CCSN simulations.

We describe the implementation and performance of the (mathrm {P}^{3}mathrm{T}) (Particle-Particle Particle-Tree) scheme for simulating dense stellar systems. In (mathrm{P}^{3}mathrm{T}), the force experienced by a particle is split into short-range and long-range contributions. Short-range forces are evaluated by direct summation and integrated with the fourth order Hermite predictor-corrector method with the block timesteps. For long-range forces, we use a combination of the Barnes-Hut tree code and the leapfrog integrator. The tree part of our simulation environment is accelerated using graphical processing units (GPU), whereas the direct summation is carried out on the host CPU. Our code gives excellent performance and accuracy for star cluster simulations with a large number of particles even when the core size of the star cluster is small.

We introduce phew (Parallel HiErarchical Watershed), a new segmentation algorithm to detect structures in astrophysical fluid simulations, and its implementation into the adaptive mesh refinement (AMR) code ramses. phew works on the density field defined on the adaptive mesh, and can thus be used on the gas density or the dark matter density after a projection of the particles onto the grid. The algorithm is based on a ‘watershed’ segmentation of the computational volume into dense regions, followed by a merging of the segmented patches based on the saddle point topology of the density field. phew is capable of automatically detecting connected regions above the adopted density threshold, as well as the entire set of substructures within. Our algorithm is fully parallel and uses the MPI library. We describe in great detail the parallel algorithm and perform a scaling experiment which proves the capability of phew to run efficiently on massively parallel systems. Future work will add a particle unbinding procedure and the calculation of halo properties onto our segmentation algorithm, thus expanding the scope of phew to genuine halo finding.

The Roche lobe geometry is important to understand and study the properties of the mass-losing component in a semi-detached binary system. However it is not easy to calculate accurately, and existing tables usually do not include the parameters of the binary system under study, nor do they allow for non-synchronous rotation.

A calculator for properties of the Roche lobe is presented in two formats. An easy-to-use Java version has a graphic interface, and a Fortran 90 version has a command line interface. The Fortran version allows for easy modifications by the user. Both versions have two basic output options: one provides values of a set of various quantities (such as the Lagrange points along the binary axis, and area and volume of the Roche lobe); the second provides (R(theta,phi)), the distance from the stellar center to the stellar surface for any specified polar angle. A single set of input parameters can be entered directly or a large set of input parameters can be specified in a text file. The calculator includes the options to have non-synchronous rotation of the star, or to have the star underfill its Roche lobe. It can be used to calculate Roche lobe properties for the case of elliptical orbits, with some restrictions.

We present a convenient software tool for quickly and accurately calculating Roche lobe properties for mass ratio in the range 0.01 to 100, for Roche lobe fill-out factor in the range 0.1 to 1.0, and for dimensionless rotation rate of the star in the range 0.1 to 2.0. This will allow anyone working with a binary star system to obtain the Roche lobe or stellar surface geometry for their system.

ChaNGa is an N-body cosmology simulation application implemented using Charm++. In this paper, we present the parallel design of ChaNGa and address many challenges arising due to the high dynamic ranges of clustered datasets. We propose optimizations based on adaptive techniques. We evaluate the performance of ChaNGa on highly clustered datasets: a (z sim0) snapshot of a 2 billion particle realization of a 25 Mpc volume, and a 52 million particle multi-resolution realization of a dwarf galaxy. For the 25 Mpc volume, we show strong scaling on up to 128K cores of Blue Waters. We also demonstrate scaling up to 128K cores of a multi-stepping run of the 2 billion particle simulation. While the scaling of the multi-stepping run is not as good as single stepping, the throughput at 128K cores is greater by a factor of 2. We also demonstrate strong scaling on up to 512K cores of Blue Waters for two large, uniform datasets with 12 and 24 billion particles.

The first stars are believed to have formed a few hundred million years after the big bang in so-called dark matter minihalos with masses (sim 10^{6}mbox{ M}_{odot}). Their radiation lit up the Universe for the first time, and the supernova explosions that ended their brief lives enriched the intergalactic medium with the first heavy elements. Influenced by their feedback, the first galaxies assembled in halos with masses (sim10^{8}mbox{ M}_{odot}), and hosted the first metal-enriched stellar populations. In this review, I summarize the theoretical progress made in the field of high-redshift star and galaxy formation since the turn of the millennium, with an emphasis on numerical simulations. These have become the method of choice to understand the multi-scale, multi-physics problem posed by structure formation in the early Universe. In the first part of the review, I focus on the formation of the first stars in minihalos - in particular the post-collapse phase, where disk fragmentation, protostellar evolution, and radiative feedback become important. I?also discuss the influence of additional physical processes, such as magnetic fields and streaming velocities. In the second part of the review, I summarize the various feedback mechanisms exerted by the first stars, followed by a discussion of the first galaxies and the various physical processes that operate in them.

The general consensus in the N-body community is that statistical results of an ensemble of collisional N-body simulations are accurate, even though individual simulations are not. A way to test this hypothesis is to make a direct comparison of an ensemble of solutions obtained by conventional methods with an ensemble of true solutions. In order to make this possible, we wrote an N-body code called Brutus, that uses arbitrary-precision arithmetic. In combination with the Bulirsch-Stoer method, Brutus is able to obtain converged solutions, which are true up to a specified number of digits.

We perform simulations of democratic 3-body systems, where after a sequence of resonances and ejections, a final configuration is reached consisting of a permanent binary and an escaping star. We do this with conventional double-precision methods, and with Brutus; both have the same set of initial conditions and initial realisations. The ensemble of solutions from the conventional simulations is compared directly to that of the converged simulations, both as an ensemble and on an individual basis to determine the distribution of the errors.

We find that on average at least half of the conventional simulations diverge from the converged solution, such that the two solutions are microscopically incomparable. For the solutions which have not diverged significantly, we observe that if the integrator has a bias in energy and angular momentum, this propagates to a bias in the statistical properties of the binaries. In the case when the conventional solution has diverged onto an entirely different trajectory in phase-space, we find that the errors are centred around zero and symmetric; the error due to divergence is unbiased, as long as the time-step parameter, (etale2^{-5}) and when simulations which violate energy conservation by more than 10% are excluded. For resonant 3-body interactions, we conclude that the statistical results of an ensemble of conventional solutions are indeed accurate.