Cactus: Issues for Sustainable Simulation Software

The Cactus Framework is an open-source, modular, portable programming environment for the collaborative development and deployment of scientific applications using high-performance computing. Its roots reach back to 1996 at the National Center for Supercomputer Applications and the Albert Einstein Institute in Germany, where its development jumpstarted. Since then, the Cactus framework has witnessed major changes in hardware infrastructure as well as its own community. This paper describes its endurance through these past changes and, drawing upon lessons from its past, also discusses futureComment: submitted to the Workshop on Sustainable Software for Science: Practice and Experiences 201

Non-symmetric trapped surfaces in the Schwarzschild and Vaidya spacetimes

Marginally trapped surfaces (MTSs) are commonly used in numerical relativity to locate black holes. For dynamical black holes, it is not known generally if this procedure is sufficiently reliable. Even for Schwarzschild black holes, Wald and Iyer constructed foliations which come arbitrarily close to the singularity but do not contain any MTSs. In this paper, we review the Wald-Iyer construction, discuss some implications for numerical relativity, and generalize to the well known Vaidya spacetime describing spherically symmetric collapse of null dust. In the Vaidya spacetime, we numerically locate non-spherically symmetric trapped surfaces which extend outside the standard spherically symmetric trapping horizon. This shows that MTSs are common in this spacetime and that the event horizon is the most likely candidate for the boundary of the trapped region.Comment: 4 pages, 3 figures; v2: minor modifications; v3: clarified conclusion

Numerical relativity with characteristic evolution, using six angular patches

The characteristic approach to numerical relativity is a useful tool in evolving gravitational systems. In the past this has been implemented using two patches of stereographic angular coordinates. In other applications, a six-patch angular coordinate system has proved effective. Here we investigate the use of a six-patch system in characteristic numerical relativity, by comparing an existing two-patch implementation (using second-order finite differencing throughout) with a new six-patch implementation (using either second- or fourth-order finite differencing for the angular derivatives). We compare these different codes by monitoring the Einstein constraint equations, numerically evaluated independently from the evolution. We find that, compared to the (second-order) two-patch code at equivalent resolutions, the errors of the second-order six-patch code are smaller by a factor of about 2, and the errors of the fourth-order six-patch code are smaller by a factor of nearly 50.Comment: 12 pages, 5 figures, submitted to CQG (special NFNR issue

Accurate Evolutions of Orbiting Binary Black Holes

We present a detailed analysis of binary black hole evolutions in the last orbit and demonstrate consistent and convergent results for the trajectories of the individual bodies. The gauge choice can significantly affect the overall accuracy of the evolution. It is possible to reconcile certain gauge-dependent discrepancies by examining the convergence limit. We illustrate these results using an initial data set recently evolved by Brügmann et al. [Phys. Rev. Lett. 92, 211101 (2004)]. For our highest resolution and most accurate gauge, we estimate the duration of this data set's last orbit to be approximately 59MADM

A Massive Data Parallel Computational Framework for Petascale/Exascale Hybrid Computer Systems

Heterogeneous systems are becoming more common on High Performance Computing (HPC) systems. Even using tools like CUDA and OpenCL it is a non-trivial task to obtain optimal performance on the GPU. Approaches to simplifying this task include Merge (a library based framework for heterogeneous multi-core systems), Zippy (a framework for parallel execution of codes on multiple GPUs), BSGP (a new programming language for general purpose computation on the GPU) and CUDA-lite (an enhancement to CUDA that transforms code based on annotations). In addition, efforts are underway to improve compiler tools for automatic parallelization and optimization of affine loop nests for GPUs and for automatic translation of OpenMP parallelized codes to CUDA. In this paper we present an alternative approach: a new computational framework for the development of massively data parallel scientific codes applications suitable for use on such petascale/exascale hybrid systems built upon the highly scalable Cactus framework. As the first non-trivial demonstration of its usefulness, we successfully developed a new 3D CFD code that achieves improved performance.Comment: Parallel Computing 2011 (ParCo2011), 30 August -- 2 September 2011, Ghent, Belgiu

Neutrino-driven Turbulent Convection and Standing Accretion Shock Instability in Three-Dimensional Core-Collapse Supernovae

We conduct a series of numerical experiments into the nature of three-dimensional (3D) hydrodynamics in the postbounce stalled-shock phase of core-collapse supernovae using 3D general-relativistic hydrodynamic simulations of a $27$-$M_\odot$ progenitor star with a neutrino leakage/heating scheme. We vary the strength of neutrino heating and find three cases of 3D dynamics: (1) neutrino-driven convection, (2) initially neutrino-driven convection and subsequent development of the standing accretion shock instability (SASI), (3) SASI dominated evolution. This confirms previous 3D results of Hanke et al. 2013, ApJ 770, 66 and Couch & Connor 2014, ApJ 785, 123. We carry out simulations with resolutions differing by up to a factor of $\sim$4 and demonstrate that low resolution is artificially favorable for explosion in the 3D convection-dominated case, since it decreases the efficiency of energy transport to small scales. Low resolution results in higher radial convective fluxes of energy and enthalpy, more fully buoyant mass, and stronger neutrino heating. In the SASI-dominated case, lower resolution damps SASI oscillations. In the convection-dominated case, a quasi-stationary angular kinetic energy spectrum $E(\ell)$ develops in the heating layer. Like other 3D studies, we find $E(\ell) \propto \ell^{-1}$ in the "inertial range," while theory and local simulations argue for $E(\ell) \propto \ell^{-5/3}$. We argue that current 3D simulations do not resolve the inertial range of turbulence and are affected by numerical viscosity up to the energy containing scale, creating a "bottleneck" that prevents an efficient turbulent cascade.Comment: 24 pages, 15 figures. Accepted for publication in The Astrophysical Journal. Added one figure and made minor modifications to text according to suggestions from the refere

From Physics Model to Results: An Optimizing Framework for Cross-Architecture Code Generation

Starting from a high-level problem description in terms of partial differential equations using abstract tensor notation, the Chemora framework discretizes, optimizes, and generates complete high performance codes for a wide range of compute architectures. Chemora extends the capabilities of Cactus, facilitating the usage of large-scale CPU/GPU systems in an efficient manner for complex applications, without low-level code tuning. Chemora achieves parallelism through MPI and multi-threading, combining OpenMP and CUDA. Optimizations include high-level code transformations, efficient loop traversal strategies, dynamically selected data and instruction cache usage strategies, and JIT compilation of GPU code tailored to the problem characteristics. The discretization is based on higher-order finite differences on multi-block domains. Chemora's capabilities are demonstrated by simulations of black hole collisions. This problem provides an acid test of the framework, as the Einstein equations contain hundreds of variables and thousands of terms.Comment: 18 pages, 4 figures, accepted for publication in Scientific Programmin

Moving black holes via singularity excision

We present a singularity excision algorithm appropriate for numerical simulations of black holes moving throughout the computational domain. The method is an extension of the excision procedure previously used to obtain stable simulations of single, non-moving black holes. The excision procedure also shares elements used in recent work to study the dynamics of a scalarfield in the background of a single, boosted black hole. The robustness of our excision method is tested with single black-hole evolutions using a coordinate system in which the coordinate location of the black hole, and thus the excision boundary, moves throughout the computational domain.Comment: 9 pages and 11 figure