Daniel R. Reynolds

ENZO: AN ADAPTIVE MESH REFINEMENT CODE FOR ASTROPHYSICS

This paper describes the open-source code Enzo, which uses block-structured adaptive mesh refinement to provide high spatial and temporal resolution for modeling astrophysical fluid flows. The code is Cartesian, can be run in one, two, and three dimensions, and supports a wide variety of physics including hydrodynamics, ideal and non-ideal magnetohydrodynamics, N-body dynamics (and, more broadly, self-gravity of fluids and particles), primordial gas chemistry, optically thin radiative cooling of primordial and metal-enriched plasmas (as well as some optically-thick cooling models), radiation transport, cosmological expansion, and models for star formation and feedback in a cosmological context. In addition to explaining the algorithms implemented, we present solutions for a wide range of test problems, demonstrate the codes parallel performance, and discuss the Enzo collaborations code development methodology.

Multiphysics simulations: Challenges and opportunities

We consider multiphysics applications from algorithmic and architectural perspectives, where “algorithmic” includes both mathematical analysis and computational complexity, and “architectural” includes both software and hardware environments. Many diverse multiphysics applications can be reduced, en route to their computational simulation, to a common algebraic coupling paradigm. Mathematical analysis of multiphysics coupling in this form is not always practical for realistic applications, but model problems representative of applications discussed herein can provide insight. A variety of software frameworks for multiphysics applications have been constructed and refined within disciplinary communities and executed on leading-edge computer systems. We examine several of these, expose some commonalities among them, and attempt to extrapolate best practices to future systems. From our study, we summarize challenges and forecast opportunities.

Cosmological radiative transfer comparison project - II. The radiation-hydrodynamic tests

The development of radiation hydrodynamical methods that are able to follow gas dynamics and radiative transfer (RT) self-consistently is key to the solution of many problems in numerical astrophysics. Such fluid flows are highly complex, rarely allowing even for approximate analytical solutions against which numerical codes can be tested. An alternative validation procedure is to compare different methods against each other on common problems, in order to assess the robustness of the results and establish a range of validity for the methods. Previously, we presented such a comparison for a set of pure RT tests (i.e. for fixed, non-evolving density fields). This is the second paper of the Cosmological Radiative Transfer Comparison Project, in which we compare nine independent RT codes directly coupled to gas dynamics on three relatively simple astrophysical hydrodynamics problems: (i) the expansion of an H ii region in a uniform medium, (ii) an ionization front in a 1/r2 density profile with a flat core and (iii) the photoevaporation of a uniform dense clump. Results show a broad agreement between the different methods and no big failures, indicating that the participating codes have reached a certain level of maturity and reliability. However, many details still do differ, and virtually every code has showed some shortcomings and has disagreed, in one respect or another, with the majority of the results. This underscores the fact that no method is universal and all require careful testing of the particular features which are most relevant to the specific problem at hand.

A fully implicit numerical method for single-fluid resistive magnetohydrodynamics

We present a nonlinearly implicit, conservative numerical method for integration of the single-fluid resistive MHD equations. The method uses a high-order spatial discretization that preserves the solenoidal property of the magnetic field. The fully coupled PDE system is solved implicitly in time, providing for increased interaction between physical processes as well as additional stability over explicit-time methods. A high-order adaptive time integration is employed, which in many cases enables time steps ranging from one to two orders of magnitude larger than those constrained by the explicit CFL condition. We apply the solution method to illustrative examples relevant to stiff magnetic fusion processes which challenge the efficiency of explicit methods. We provide computational evidence showing that for such problems the method is comparably accurate with explicit-time simulations, while providing a significant runtime improvement due to its increased temporal stability.

FULLY COUPLED SIMULATION OF COSMIC REIONIZATION. II. RECOMBINATIONS, CLUMPING FACTORS, AND THE PHOTON BUDGET FOR REIONIZATION

We use a fully self-consistent cosmological simulation including dark matter dynamics, multispecies hydrodynamics, chemical ionization, flux limited diffusion radiation transport, and a parameterized model of star formation and feedback (thermal and radiative) to investigate the epoch of hydrogen reionization in detail. Our numerical method is scalable with respect to the number of radiation sources, size of the mesh, and the number of computer processors employed, and is described in Paper I of this series. In this the first of several application papers, we investigate the mechanics of reionization from stellar sources forming in high-z galaxies, the utility of various formulations for the gas clumping factor on accurately estimating the effective recombination time in the IGM, and the photon budget required to achieve reionization. We also test the accuracy of the static and time-dependent models of Madau et al. as predictors of reionization completion/maintenance. We simulate a WMAP7 CDM cosmological model in a 20 Mpc comoving cube, resolved with 800 3 uniform fluid cells and dark matter particles. By tuning our star formation recipe to approximately match the observed high redshift star formation rate density and galaxy luminosity function, we have created a fully coupled radiation hydrodynamical realization of hydrogen reionization which begins to ionize at z 10 and completes at z 5:8 without further tuning. The complicated events during reionization that lead to this number can be generally described as inside-out, but in reality the narrative depends on the level of ionization of the gas one attributes to as ionized. We find that roughly 2 ionizing photons per H atom are required to convert the neutral IGM to a highly ionized state, which supports the “photon starved” reionization scenario discussed by Bolton & Haehnelt. We find that the formula for the ionizing photon production rate _ Nion(z) needed to maintain the IGM in an ionized state derived by Madau et al. should not be used to predict the epoch of reionization completion because it ignores history-dependent terms in the global ionization balance which are not ignorable. We find that the time-dependent model for the ionized volume fraction QHII is more predictive, but overestimates the redshift of reionization completion zreion by z 1. We propose a revised formulation of the time-dependent model which agrees with our simulation to high accuracy. Finally, we use our simulation data to estimate a globally averaged ionizing escape fraction due to circumgalactic gas resolved on our mesh fesc(CGM) 0:7. Subject headings: cosmology: theory ‐ intergalactic medium ‐ reionization ‐ large-scale structure of universe ‐ methods: numerical ‐ radiative transfer

Fully coupled simulation of cosmic reionization. I. numerical methods and tests

The Astrophysical Journal Supplement Series, 216:16 (24pp), 2015 January C 2015. doi:10.1088/0067-0049/216/1/16 The American Astronomical Society. All rights reserved. FULLY COUPLED SIMULATION OF COSMIC REIONIZATION. I. NUMERICAL METHODS AND TESTS Michael L. Norman 1,2 , Daniel R. Reynolds 3 , Geoffrey C. So 1 , Robert P. Harkness 2,4 , and John H. Wise 5 CASS, University of California, San Diego, 9500 Gilman Drive La Jolla, CA 92093-0424, USA SDSC, University of California, San Diego, 9500 Gilman Drive La Jolla, CA 92093-0505, USA 3 Southern Methodist University, 6425 Boaz Lane, Dallas, TX 75205, USA 4 NICS, Oak Ridge National Laboratory, 1 Bethel Valley Road, Oak Ridge, TN 37831, USA 5 Center for Relativistic Astrophysics, Georgia Institute of Technology, 837 State Street, Atlanta, GA 30332, USA Received 2013 June 8; accepted 2014 November 23; published 2015 January 9 ABSTRACT We describe an extension of the Enzo code to enable fully coupled radiation hydrodynamical simulation of inhomogeneous reionization in large ∼(100 Mpc) 3 cosmological volumes with thousands to millions of point sources. We solve all dynamical, radiative transfer, thermal, and ionization processes self-consistently on the same mesh, as opposed to a postprocessing approach which coarse-grains the radiative transfer. We do, however, employ a simple subgrid model for star formation which we calibrate to observations. The numerical method presented is a modification of an earlier method presented in Reynolds et al. differing principally in the operator splitting algorithm we use to advance the system of equations. Radiation transport is done in the gray flux-limited diffusion (FLD) approximation, which is solved by implicit time integration split off from the gas energy and ionization equations, which are solved separately. This results in a faster and more robust scheme for cosmological applications compared to the earlier method. The FLD equation is solved using the hypre optimally scalable geometric multigrid solver from LLNL. By treating the ionizing radiation as a grid field as opposed to rays, our method is scalable with respect to the number of ionizing sources, limited only by the parallel scaling properties of the radiation solver. We test the speed and accuracy of our approach on a number of standard verification and validation tests. We show by direct comparison with Enzo’s adaptive ray tracing method Moray that the well-known inability of FLD to cast a shadow behind opaque clouds has a minor effect on the evolution of ionized volume and mass fractions in a reionization simulation validation test. We illustrate an application of our method to the problem of inhomogeneous reionization in a 80 Mpc comoving box resolved with 3200 3 Eulerian grid cells and dark matter particles. Key words: cosmology: theory – methods: numerical – radiative transfer redshift interval 6 z 10 using the Hubble Space Telescope support the galaxy reionizer hypothesis, with the caveat that the faint end of the luminosity function which contributes substantially to the number of ionizing photons has not yet been measured (Robertson et al. 2010; Bouwens et al. 2012). Given the paucity of observational information about the pro- cess of cosmic reionization, researchers have resorted to theory and numerical simulation to fill in the blanks. As reviewed by Trac & Gnedin (2011), progress in this area has been dramatic, driven by a synergistic interplay between semi-analytic ap- proaches and numerical simulations. The combination of these two approaches have converged on a qualitative picture of how H reionization proceeds assuming the primary ionizing sources are young, star-forming galaxies. The physics of the reioniza- tion process is determined by the physics of the sources and sinks of ionizing radiation in an expanding universe. Adopting the ΛCDM model of structure formation, galaxies form hierar- chically through the merger of dark matter halos. The structure and evolution of the dark matter density field is now well un- derstood through ultra-high-resolution numerical N-body sim- ulations (Springel et al. 2005; Klypin et al. 2011) and through analytic models based on these simulations (Cooray & Sheth 2002). By making certain assumptions about how ionizing light traces mass and the dynamics of H ii regions, a basic picture of the reionization process has emerged (Furlanetto et al. 2004, 2006; Iliev et al. 2006b; Zahn et al. 2007) that is confirmed by detailed numerical simulations; e.g., Zahn et al. (2011). The basic picture is that galaxies form in the peaks of the dark matter density field and drive expanding H ii regions into their surroundings by virtue of the UV radiation emitted from young, massive stars. These H ii regions are initially isolated, 1. INTRODUCTION The epoch of reionization (EoR) is a current frontier of cos- mological research both observationally and theoretically. Ob- servations constrain the transition from a largely neutral inter- galactic medium (IGM) of primordial gas to a largely ionized one (singly ionized H and He) to the redshift interval z ∼ 11–6, which is a span of roughly 500 Myr. The completion of H reion- ization by z ≈ 6 is firmly established through quasar absorption line studies to luminous, high-redshift quasars which exhibit Lyα Gunn–Peterson absorption troughs (Fan et al. 2006). The precise onset of H reionization (presumably tied to the formation of the first luminous ionizing sources) is presently unknown ob- servationally, however, cosmic microwave background (CMB) measurements of the Thomson optical depth to the surface of last scattering by the Wilkinson Microwave Anisotropy Probe (WMAP) and Planck satellites indicates that the IGM was sub- stantially ionized by z ∼ 10 (Spergel et al. 2003; Komatsu et al. 2009; Jarosik et al. 2011; Planck Collaboration et al. 2014). Since the optical depth measurement is redshift-integrated and averaged over the sky, the CMB observations provide no infor- mation about how reionization proceeded or the nature of the radiation sources that caused it. It is generally believed that reionization begins with the formation of Population III stars at z ∼ 20–30 (Abel et al. 2002; Yoshida et al. 2003; Bromm & Larson 2004; Sokasian et al. 2004), but that soon the ionizing photon budget becomes dominated by young, star forming galaxies (see e.g., Wise et al. 2012; Xu et al. 2013), and to a lesser extent by the first quasars (Madau et al. 1999; Bolton & Haehnelt 2007; Haardt & Madau 2012; Becker & Bolton 2013). Observations of galaxies in the

Implicit solvers for large-scale nonlinear problems

Computational scientists are grappling with increasingly complex, multi-rate applications that couple such physical phenomena as fluid dynamics, electromagnetics, radiation transport, chemical and nuclear reactions, and wave and material propagation in inhomogeneous media. Parallel computers with large storage capacities are paving the way for high-resolution simulations of coupled problems; however, hardware improvements alone will not prove enough to enable simulations based on brute-force algorithmic approaches. To accurately capture nonlinear couplings between dynamically relevant phenomena, often while stepping over rapid adjustments to quasi-equilibria, simulation scientists are increasingly turning to implicit formulations that require a discrete nonlinear system to be solved for each time step or steady state solution. Recent advances in iterative methods have made fully implicit formulations a viable option for solution of these large-scale problems. In this paper, we overview one of the most effective iterative methods, Newton-Krylov, for nonlinear systems and point to software packages with its implementation. We illustrate the method with an example from magnetically confined plasma fusion and briefly survey other areas in which implicit methods have bestowed important advantages, such as allowing high-order temporal integration and providing a pathway to sensitivity analyses and optimization. Lastly, we overview algorithm extensions under development motivated by current SciDAC applications.

Self-consistent solution of cosmological radiation-hydrodynamics and chemical ionization

We consider a PDE system comprising compressible hydrodynamics, flux-limited diffusion radiation transport and chemical ionization kinetics in a cosmologically-expanding universe. Under an operator-split framework, the cosmological hydrodynamics equations are solved through the piecewise parabolic method, as implemented in the Enzo community hydrodynamics code. The remainder of the model, including radiation transport, chemical ionization kinetics, and gas energy feedback, form a stiff coupled PDE system, which we solve using a fully-implicit inexact Newton approach, and which forms the crux of this paper. The inner linear Newton systems are solved using a Schur complement formulation, and employ a multigrid-preconditioned conjugate gradient solver for the inner Schur systems. We describe this approach and provide results on a suite of test problems, demonstrating its accuracy, robustness, and scalability to very large problems.

Operator-Based Preconditioning of Stiff Hyperbolic Systems

We introduce an operator-based scheme for preconditioning stiff components encountered in implicit methods for hyperbolic systems of PDEs posed on regular grids. The method is based on a directional splitting of the implicit operator, followed by a characteristic decomposition of the resulting directional parts. This approach allows for the solution of any number of characteristic components, from the entire system to only the fastest, stiffness-inducing waves. We apply the preconditioning method to stiff hyperbolic systems arising in magnetohydrodynamics and gas dynamics. We then present numerical results showing that this preconditioning scheme works well on problems where the underlying stiffness results from the interaction of fast transient waves with slowly-evolving dynamics, scales well to large problem sizes and numbers of processors, and allows for additional customization based on the specific problems under study.

Efficient and automatic implementation of the adjoint state method

Combination of object-oriented programming with automatic differentiation techniques facilitates the solution of data fitting, control, and design problems driven by explicit time stepping schemes for initial-boundary value problems. The C++ class fdtd takes a complete specification of a single step, along with some associated code, and assembles from it a complete simulator, along with the linearized and adjoint simulations. The result is a (nonlinear) operator in the sense of the Hilbert Class Library (HCL), a C++ software package for optimization. The HCL operator so produced links directly with any of the HCL optimization algorithms. Moreover the performance of simulators constructed in this way is equivalent to that of optimized Fortran implementations.