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.     

Progress in development of HEDP capabilities in FLASH’s Unsplit Staggered Mesh MHD solver

FLASH is a publicly available astrophysical community code designed to solve highly compressible multi-physics reactive flows. We are adding capabilities to FLASH that will make it an open science code for the academic HEDP community. Among many important numerical requirements, we consider the following features to be important components necessary to meet our goals for FLASH as an HEDP open toolset. First, we are developing computationally efficient time-stepping integration methods that overcome the stiffness that arises in the equations describing a physical problem when there are disparate time scales. To this end, we are adding two different time-stepping schemes to FLASH that relax the time step limit when diffusive effects are present: an explicit super-time-stepping algorithm (Alexiades et al. in Com. Num. Mech. Eng. 12:31–42, 1996) and a Jacobian-Free Newton-Krylov implicit formulation. These two methods will be integrated into a robust, efficient, and high-order accurate Unsplit Staggered Mesh MHD (USM) solver (Lee and Deane in J. Comput. Phys. 227, 2009). Second, we have implemented an anisotropic Spitzer-Braginskii conductivity model to treat thermal heat conduction along magnetic field lines. Finally, we are implementing the Biermann Battery term to account for spontaneous generation of magnetic fields in the presence of non-parallel temperature and density gradients.     

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.     

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.     

Numerical simulation of a permittivity probe for measuring the electric properties of planetary regolith and application to the near‐surface region of asteroids and comets

Abstract— We present a numerical simulation technique for the retrieval of the electric properties relative permittivity and conductivity of planetary, asteroid, and cometary regolith. Our simulation techniques aim at accompanying hardware development and conducting virtual experiments, e.g., to assess the response of arbitrary heterogeneous conductivity and permittivity distributions or to scrutinize possibilities for spatial reconstruction methods using inverse schemes. In a first step, we have developed a finite element simulation code on the basis of unstructured, adaptive triangular grids for arbitrary two‐dimensional axisymmetric distributions of conductivity and permittivity. The code is able to take into account the spatial geometry of the probe and allows for possible inductive effects. In previous studies, the non‐inductive approach has been used to convert potential and phase data into apparent material properties. By our simulations, we have shown that this approach is valid for the frequency range from 10² Hz to 10⁷ Hz and electric conductivities of 10^?8 S/m that are typical for the near‐surface region of asteroids and comets composed of chondritic materials and/or frozen volatiles such as H₂O and CO₂ ice. We prove the accuracy of our code to be better than 10%, using mixed types of boundary conditions and present a simulated vertical log through a horizontally stratified subsurface layer as a representative example of a heterogeneous distribution of the electrical properties. Resolution studies for the given electrode separation reveal that the material parameters of layers having thicknesses of less than about half the electrode spread are not reconstructible if only apparent quantities are considered. Therefore, spatial distributions of the complex sensitivity are presented having in mind a future data inversion concept that will permit the multi‐dimensional reconstruction of material parameters in heterogeneous environments.     

RAPID: A fast,high resolution,flux-conservative algorithm designed for planet–disk interactions

We describe a newly developed hydrodynamic code for studying accretion disk processes. The numerical method uses a finite volume, non-linear, Total Variation Diminishing (TVD) scheme to capture shocks and control spurious oscillations. It is second-order accurate in time and space and makes use of a FARGO-type algorithm to alleviate Courant–Friedrichs–Lewy time step restrictions imposed by the rapidly rotating inner disk region. OpenMP directives are implemented enabling faster computations on shared-memory, multi-processor machines. The resulting code is simple, fast and memory efficient. We discuss the relevant details of the numerical method and provide results of the code’s performance on standard test problems. We also include a detailed examination of the code’s performance on planetary disk–planet interactions. We show that the results produced on the standard problem setup are consistent with a wide variety of other codes.     

NBSymple,a double parallel,symplectic N-body code running on graphic processing units

We present and discuss the characteristics and performance, both in term of computational speed and precision, of a numerical code which integrates the equation of motions of N ‘particles’ interacting via Newtonian gravitation and move in an external galactic smooth field. The force evaluation on every particle is done by mean of direct summation of the contribution of all the other system’s particles, avoiding truncation error. The time integration is done with second-order and sixth-order symplectic schemes. The code, NBSymple, has been parallelized twice, by mean of the Compute Unified Device Architecture (CUDA) to make the all-pair force evaluation as fast as possible on high-performance Graphic Processing Units NVIDIA TESLA C1060, while the O(N) computations are distributed on various CPUs by mean of OpenMP Application Program. The code works both in single-precision floating point arithmetics or in double precision. The use of single-precision allows the use of the GPU performance at best but, of course, limits the precision of simulation in some critical situations. We find a good compromise in using a software reconstruction of double-precision for those variables that are most critical for the overall precision of the code. The code is available on the web site astrowww.phys.uniroma1.it/dolcetta/nbsymple.html.     

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 Code for Automated Tracing of Coronal Loops Approaching Visual Perception

We develop a new numerical code with automated feature extraction, customized for tracing of coronal loops, a method we call Oriented Coronal CUrved Loop Tracing (OCCULT), which for the first time breaks even with the results of visual tracing. The method used is based on oriented-directivity tracing of curvi-linear features, but in contrast to other general feature-extraction algorithms, it is customized for solar EUV and SXR images by taking advantage of the specific property that coronal loops have large curvature radii compared with their widths. We evaluate the performance of this new code by comparing the cumulative distribution of loop lengths, the median and maximum loop lengths, the completeness of detection, and the congruency of the detected features with other numerical codes and visual tracings. We find that the new code closely approaches the results of visual perception and outperforms the other existing numerical codes. This algorithm is useful for the 3D reconstruction of the geometry, motion, and oscillations of coronal loops, with single or stereoscopic spacecraft, as well as for modeling of the loop hydrodynamics and the coronal magnetic field.     

Non-linear hydrodynamical evolution of rotating relativistic stars: numerical methods and code tests

We present numerical hydrodynamical evolutions of rapidly rotating relativistic stars, using an axisymmetric, non-linear relativistic hydrodynamics code. We use four different high-resolution shock-capturing (HRSC) finite-difference schemes (based on approximate Riemann solvers) and compare their accuracy in preserving uniformly rotating stationary initial configurations in long-term evolutions. Among these four schemes, we find that the third-order piecewise parabolic method scheme is superior in maintaining the initial rotation law in long-term evolutions, especially near the surface of the star. It is further shown that HRSC schemes are suitable for the evolution of perturbed neutron stars and for the accurate identification (via Fourier transforms) of normal modes of oscillation. This is demonstrated for radial and quadrupolar pulsations in the non-rotating limit, where we find good agreement with frequencies obtained with a linear perturbation code. The code can be used for studying small-amplitude or non-linear pulsations of differentially rotating neutron stars, while our present results serve as testbed computations for three-dimensional general-relativistic evolution codes.     

Star formation and chemical evolution in smoothed particle hydrodynamics simulations: a statistical approach

In smoothed particle hydrodynamics (SPH) codes with a large number of particles, star formation as well as gas and metal restitution from dying stars can be treated statistically. This approach allows one to include detailed chemical evolution and gas re-ejection with minor computational effort. Here we report on a new statistical algorithm for star formation and chemical evolution, especially conceived for SPH simulations with large numbers of particles, and for parallel SPH codes.
For the sake of illustration, we also present two astrophysical simulations obtained with this algorithm, implemented into the Tree-SPH code by Lia & Carraro .
In the first simulation, we follow the formation of an individual disc-like galaxy, predict the final structure and metallicity evolution, and test resolution effects. In the second simulation we simulate the formation and evolution of a cluster of galaxies, to demonstrate the capabilities of the algorithm in investigating the chemo-dynamical evolution of galaxies and of the intergalactic medium in a cosmological context.     

Nonlinear Force-Free Magnetic Field Extrapolations: Comparison of the Grad Rubin and Wheatland Sturrock Roumeliotis Algorithm

We compare the performance of two alternative algorithms which aim to construct a force-free magnetic field given suitable boundary conditions. For this comparison, we have implemented both algorithms on the same finite element grid which uses Whitney forms to describe the fields within the grid cells. The additional use of conjugate gradient and multigrid iterations result in quite effective codes. The Grad Rubin and Wheatland Sturrock Roumeliotis algorithms both perform well for the reconstruction of a known analytic force-free field. For more arbitrary boundary conditions the Wheatland Sturrock Roumeliotis approach has some difficulties because it requires overdetermined boundary information which may include inconsistencies. The Grad Rubin code on the other hand loses convergence for strong current densities. For the example we have investigated, however, the maximum possible current density seems to be not far from the limit beyond which a force-free field cannot exist anymore for a given normal magnetic field intensity on the boundary.     

Reconstruction of RHESSI Solar Flare Images with a Forward Fitting Method

We describe a forward-fitting method that has been developed to reconstruct hard X-ray images of solar flares from the Ramaty High-Energy Solar Spectroscopic Imager (RHESSI), a Fourier imager with rotation-modulated collimators that was launched on 5 February 2002. The forward-fitting method is based on geometric models that represent a spatial map by a superposition of multiple source structures, which are quantified by circular gaussians (4 parameters per source), elliptical gaussians (6 parameters), or curved ellipticals (7 parameters), designed to characterize real solar flare hard X-ray maps with a minimum number of geometric elements. We describe and demonstrate the use of the forward-fitting algorithm. We perform some 500 simulations of rotation-modulated time profiles of the 9 RHESSI detectors, based on single and multiple source structures, and perform their image reconstruction. We quantify the fidelity of the image reconstruction, as function of photon statistics, and the accuracy of retrieved source positions, widths, and fluxes. We outline applications for which the forward-fitting code is most suitable, such as measurements of the energy-dependent altitude of energy loss near the limb, or footpoint separation during flares. Supplementary material to this paper is available in electronic form at http://dx.doi.org/10.1023/A:1022469811115     

Spotted star mapping by light curve inversion: Tests and application to HD 12545

A code for mapping the surfaces of spotted stars is developed. The concept of the code is to analyze rotational-modulated light curves. We simulate the process of reconstruction for the star surface and the results of simulation are presented. The reconstruction atrifacts caused by the ill-posed nature of the problem are deduced. The surface of the spotted component of system HD 12545 is mapped using the procedure.     

Radiation hydrodynamics using characteristics on adaptive decomposed domains for massively parallel star formation simulations

We present an algorithm for solving the radiative transfer problem on massively parallel computers using adaptive mesh refinement and domain decomposition. The solver is based on the method of characteristics which requires an adaptive raytracer that integrates the equation of radiative transfer. The radiation field is split into local and global components which are handled separately to overcome the non-locality problem. The solver is implemented in the framework of the magneto-hydrodynamics code FLASH and is coupled by an operator splitting step. The goal is the study of radiation in the context of star formation simulations with a focus on early disc formation and evolution. This requires a proper treatment of radiation physics that covers both the optically thin as well as the optically thick regimes and the transition region in particular. We successfully show the accuracy and feasibility of our method in a series of standard radiative transfer problems and two 3D collapse simulations resembling the early stages of protostar and disc formation.     

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.     

Super resolution for astronomical observations

In order to obtain detailed information from multiple telescope observations a general blind super-resolution (SR) reconstruction approach for astronomical images is proposed in this paper. A pixel-reliability-based SR reconstruction algorithm is described and implemented, where the developed process incorporates flat field correction, automatic star searching and centering, iterative star matching, and sub-pixel image registration. Images captured by the 1-m telescope at Yunnan Observatory are used to test the proposed technique. The results of these experiments indicate that, following SR reconstruction, faint stars are more distinct, bright stars have sharper profiles, and the backgrounds have higher details; thus these results benefit from the high-precision star centering and image registration provided by the developed method. Application of the proposed approach not only provides more opportunities for new discoveries from astronomical image sequences, but will also contribute to enhancing the capabilities of most spatial or ground-based telescopes.     

A Time-Transformed Leapfrog Scheme

We present a time-transformed leapfrog scheme combined with the extrapolation method to construct an integrator for orbits in N-body systems with large mass ratios. The basic idea can be used to transform any second-order differential equation into a form which may allow more efficient numerical integration. When applied to gravitating few-body systems this formulation permits extremely close two-body encounters to be considered without significant loss of accuracy. The new scheme has been implemented in a direct N-body code for simulations of super-massive binaries in galactic nuclei. In this context relativistic effects may also be included.     

Automated analysis of eclipsing binary light curves – I. EBAS – a new Eclipsing Binary Automated Solver with ebop

We present a new algorithm, Eclipsing Binary Automated Solver (EBAS), to analyse light curves of eclipsing binaries. The algorithm is designed to analyse large numbers of light curves, and is therefore based on the relatively fast ebop code. To facilitate the search for the best solution, EBAS uses two parameter transformations. Instead of the radii of the two stellar components, EBAS uses the sum of radii and their ratio, while the inclination is transformed into the impact parameter. To replace human visual assessment, we introduce a new 'alarm' goodness-of-fit statistic that takes into account correlation between neighbouring residuals. We perform extensive tests and simulations that show that our algorithm converges well, finds a good set of parameters and provides reasonable error estimation.     

Optimization code with weighting function for the reconstruction of coronal magnetic fields

We developed a code for the reconstruction of nonlinear force-free and non-force-free coronal magnetic fields. The 3D magnetic field is computed numerically with the help of an optimization principle. The force-free and non-force-free codes are compiled in one program. The force-free approach needs photospheric vector magnetograms as input. The non-force-free code additionally requires the line-of-sight integrated coronal density distribution in combination with a tomographic inversion code. Previously the optimization approach has been used to compute magnetic fields using all six boundaries of a computational box. Here we extend this method and show how the coronal magnetic field can be reconstructed only from the bottom boundary, where the boundary conditions are measured with vector magnetographs. The program is planed for use within the Stereo mission.